#include <vector>	// std::vector
#include <cmath>	 // cos(),sin(),sqrt()
#include <limits>	// std::numeric_limits<T>
#include <stdexcept> // std::domain_error()
#include "igs_ifx_common.h"
#include "igs_radial_blur.h"
namespace
{
//------------------------------------------------------------------
template <class T>
class radial_
{
public:
	radial_(
		const T *in_top, const int height, const int width, const int channels, const double xc, const double yc, const double sub_size, const int imax, const double dmax, const double intensity /* 平均値ぼかし強度 */
		,
		const double radius /* 平均値ぼかしの始まる半径 */
		)
		: in_top_(in_top), hh_(height), ww_(width), cc_(channels), xc_(xc), yc_(yc), sub_size_(sub_size), imax_(imax), dmax_(dmax), intensity_(intensity), radius_(radius)
	{
	}
	void pixel_value(
		const T *in_current_pixel, const int xx, const int yy, const int z1, const int z2, const double ref_increase_val, const double ref_decrease_val, const double each_pixel_blur_ratio, T *result_pixel)
	{
		/* Pixel位置(0.5 1.5 2.5 ...) */
		const double xp = static_cast<double>(xx) + 0.5;
		const double yp = static_cast<double>(yy) + 0.5;
 
		/* 中心からPixel位置へのベクトルと長さ */
		const double xv = xp - this->xc_;
		const double yv = yp - this->yc_;
		const double dist = sqrt(xv * xv + yv * yv);
 
		/* 指定半径の範囲内なら何もしない */
		if (dist <= this->radius_) {
			for (int zz = z1; zz <= z2; ++zz) {
				result_pixel[zz] = in_current_pixel[zz];
			}
			return;
		}
 
		/* 中心からPixel位置への単位ベクトル */
		const double cosval = xv / dist;
		const double sinval = yv / dist;
 
		/* Radial方向のSamplingの開始位置と終了位置 */
		double scale = this->intensity_;
		if (0.0 <= each_pixel_blur_ratio) {
			scale *= each_pixel_blur_ratio;
		}
		const double length = (dist - this->radius_) *
							  scale;
		const double count_half = floor(
			length / 2.0 / this->sub_size_);
		const double sub_sta = -this->sub_size_ * count_half;
		const double sub_end = this->sub_size_ * count_half;
 
		/* 積算値と積算回数 */
		std::vector<double> accum_val(this->cc_);
		int accum_counter = 0;
 
		/* 円周接線方向Samplingの相対位置 */
		for (double ss = this->sub_size_ / 2.0 - 0.5;
			 ss < 0.5; ss += this->sub_size_) {
 
			/* 円周接線方向Sampling位置 */
			const double xps = xp + ss * sinval;
			const double yps = yp + ss * cosval;
 
			/* 中心からSampling位置へのベクトルと長さ */
			const double xvs = xps - this->xc_;
			const double yvs = yps - this->yc_;
			const double dists = sqrt(xvs * xvs + yvs * yvs);
 
			/* 中心からSampling位置への単位ベクトル */
			const double cosvals = xvs / dists;
			const double sinvals = yvs / dists;
 
			/* Radial方向のSampling */
			for (double tt = sub_sta;
				 tt <= sub_end; tt += this->sub_size_) {
				/* Sampling位置からPixel位置を得る */
				int xi = static_cast<int>(xps + tt * cosvals);
				int yi = static_cast<int>(yps + tt * sinvals);
 
				/* clamp */
				xi = (xi < 0) ? 0 : ((this->ww_ <= xi) ? this->ww_ - 1 : xi);
				yi = (yi < 0) ? 0 : ((this->hh_ <= yi) ? this->hh_ - 1 : yi);
 
				/* 画像のPixel位置 */
				const T *in_current = this->in_top_ +
									  this->cc_ * this->ww_ * yi +
									  this->cc_ * xi;
 
				/* 積算 */
				for (int zz = z1; zz <= z2; ++zz) {
					accum_val[zz] +=
						static_cast<double>(in_current[zz]);
				}
				++accum_counter;
			}
		}
		/* 積算しなかったとき(念のためのCheck) */
		if (accum_counter <= 0) {
			for (int zz = z1; zz <= z2; ++zz) {
				result_pixel[zz] = in_current_pixel[zz];
			}
			return;
		}
		/* ここで画像Pixelに保存 */
		for (int zz = z1; zz <= z2; ++zz) {
			accum_val[zz] /= static_cast<double>(accum_counter);
 
			if ((0 <= ref_increase_val) &&
				(in_current_pixel[zz] < accum_val[zz])) {
				/* 増分のみMask! */
				accum_val[zz] =
					static_cast<double>(in_current_pixel[zz]) +
					(accum_val[zz] - in_current_pixel[zz]) *
						ref_increase_val;
			} else if ((0 <= ref_decrease_val) &&
					   (accum_val[zz] < in_current_pixel[zz])) {
				/* 減分のみMask! */
				accum_val[zz] =
					static_cast<double>(in_current_pixel[zz]) +
					(accum_val[zz] - in_current_pixel[zz]) *
						ref_decrease_val;
			}
 
			accum_val[zz] += 0.5; /* 誤差対策 */
			if (this->dmax_ < accum_val[zz]) {
				result_pixel[zz] = static_cast<T>(this->imax_);
			} else if (accum_val[zz] < 0) {
				result_pixel[zz] = 0;
			} else {
				result_pixel[zz] = static_cast<T>(accum_val[zz]);
			}
		}
	}
 
private:
	radial_() {}
 
	const T *in_top_;
	const int hh_;
	const int ww_;
	const int cc_;
	const double xc_;
	const double yc_;
	const double sub_size_;
	const int imax_;
	const double dmax_;
	const double intensity_;
	const double radius_;
 
	/* copy constructorを無効化 */
	radial_(const radial_ &);
 
	/* 代入演算子を無効化 */
	radial_ &operator=(const radial_ &);
};
template <class IT, class RT>
void radial_convert_template_(
	const IT *in, const int margin /* 参照画像(in)がもつ余白 */
 
	,
	const RT *ref /* 求める画像(out)と同じ高さ、幅、チャンネル数 */
	,
	const int ref_bits, const int ref_mode // R,G,B,A,luminance
 
	,
	IT *out, const int hh /* 求める画像(out)の高さ */
	,
	const int ww /* 求める画像(out)の幅 */
	,
	const int cc, const double xc, const double yc, const double intensity /* 強度。ゼロより大きく2以下 */
	/* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */
	,
	const double radius /* 平均値ぼかしの始まる半径 */
	,
	const int sub_div /* 1ならJaggy、2以上はAntialias */
	,
	const bool alpha_rendering_sw)
{
	ref_bits; // for warning
 
	/* 強度のないとき、または、サブ分割がないとき、なにもしない */
	if (intensity <= 0.0 || sub_div <= 0) {
		return;
	}
 
	radial_<IT> cl_ra(
		in, hh + margin * 2, ww + margin * 2, cc, xc + margin, yc + margin, 1.0 / sub_div, std::numeric_limits<IT>::max() /* サンプリング値 */
		,
		static_cast<double>(std::numeric_limits<IT>::max()), intensity, radius);
 
#if defined RGBA_ORDER_OF_TOONZ6
	const int z1 = igs::image::rgba::blu;
	const int z2 = igs::image::rgba::red;
#elif defined RGBA_ORDER_OF_OPENGL
	const int z1 = igs::image::rgba::red;
	const int z2 = igs::image::rgba::blu;
#else
	Must be define / DRGBA_ORDER_OF_TOONZ6 or / DRGBA_ORDER_OF_OPENGL
#endif
	const IT *p_in = in + margin * (ww + margin * 2) * cc + margin * cc;
	IT *pout = out;
 
	if (0 == ref) { /* 参照なし */
		if (igs::image::rgba::siz == cc) {
			using namespace igs::image::rgba;
			if (alpha_rendering_sw) { /* Alphaも処理する */
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc) {
						cl_ra.pixel_value(p_in, xx, yy, alp, alp, -1.0, -1.0, -1.0, pout);
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ra.pixel_value(p_in, xx, yy, z1, z2, -1.0, -1.0, -1.0, pout);
					}
				}
			} else { /* Alpha処理しない、RGB増分をAlphaでMaskする */
				const unsigned int val_max = std::numeric_limits<IT>::max();
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc) {
						pout[alp] = p_in[alp];
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						/* Alpha値あるならRGB増分のみAlphaでMaskする */
						cl_ra.pixel_value(
							p_in, xx, yy, z1, z2, static_cast<double>(p_in[alp]) / val_max, -1.0, -1.0, pout);
					}
				}
			}
		} else { /* Alphaがない, RGB/Grayscale... */
			for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
				for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc) {
					cl_ra.pixel_value(p_in, xx, yy, 0, cc - 1, -1.0, -1.0, -1.0, pout);
				}
			}
		}
	} else { /* 参照あり */
		const RT *refe = ref;
		const int r_max = std::numeric_limits<RT>::max();
		if (igs::image::rgba::siz == cc) {
			using namespace igs::image::rgba;
			if (alpha_rendering_sw) { /* Alpha処理する */
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc, refe += cc) {
						const double refv = igs::color::ref_value(
							refe, cc, r_max, ref_mode);
						if (0 == refv) {
							for (int zz = 0; zz < cc; ++zz) {
								pout[zz] = p_in[zz];
							}
							continue;
						}
 
						cl_ra.pixel_value(p_in, xx, yy, alp, alp, -1.0, -1.0, refv, pout);
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ra.pixel_value(p_in, xx, yy, z1, z2, -1.0, -1.0, refv, pout);
					}
				}
			} else { /* Alpha処理しない、RGB増分をAlphaでMaskする */
				const unsigned int val_max = std::numeric_limits<IT>::max();
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc, refe += cc) {
						const double refv = igs::color::ref_value(
							refe, cc, r_max, ref_mode);
						if (0 == refv) {
							for (int zz = 0; zz < cc; ++zz) {
								pout[zz] = p_in[zz];
							}
							continue;
						}
 
						pout[alp] = p_in[alp];
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ra.pixel_value(p_in, xx, yy, z1, z2, static_cast<double>(p_in[alp]) / val_max, -1.0, refv, pout);
					}
				}
			}
		} else { /* Alphaがない, RGB/Grayscale... */
			for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
				for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc, refe += cc) {
					const double refv = igs::color::ref_value(
						refe, cc, r_max, ref_mode);
					if (0 == refv) {
						for (int zz = 0; zz < cc; ++zz) {
							pout[zz] = p_in[zz];
						}
						continue;
					}
 
					cl_ra.pixel_value(p_in, xx, yy, 0, cc - 1, -1.0, -1.0, refv, pout);
				}
			}
		}
	}
}
}
//--------------------------------------------------------------------
namespace
{
class twist_
{
public:
	twist_(
		const int height, const double xc, const double yc, const double twist_radian, const double twist_radius)
		: xc_(xc), yc_(yc), twist_xp_(0), twist_yp_(0), twist_radian_(twist_radian), twist_radius_(twist_radius), current_xp_(0), current_yp_(0), current_radian_(0), current_cos_(0), current_sin_(0), current_radius_(0)
	{
		/* twist_radiusは、Twist曲線上半径1のPixel画像上半径 */
		/* twist_radiusはゼロのとき、高さの半分を設定する */
		if (this->twist_radius_ <= 0.0) {
			this->twist_radius_ = height / 2.0;
		}
		/* ゼロだと計算不能になる */
		if (this->twist_radius_ <= 0.0) {
			throw std::domain_error(
				"twist_radius is equal less than zero");
		}
	}
	void current_pos(const double xp, const double yp)
	{
		/* 現在位置 */
		this->current_xp_ = xp;
		this->current_yp_ = yp;
 
		/* 中心からの距離 */
		const double xv = xp - this->xc_;
		const double yv = yp - this->yc_;
		this->current_radius_ = sqrt(xv * xv + yv * yv);
 
		/* 基準からの距離比 */
		const double tt = this->current_radius_ / this->twist_radius_;
 
		/* Twist曲線上の位置 */
		const double xx = tt * cos(tt * this->twist_radian_);
		const double yy = tt * sin(tt * this->twist_radian_);
 
		/* Twist曲線の回転角度 */
		this->current_radian_ = atan2(yv, xv) - atan2(yy, xx);
 
		/* Twist曲線上の位置を回転角度に傾けた位置 */
		this->current_cos_ = cos(this->current_radian_);
		this->current_sin_ = sin(this->current_radian_);
		this->twist_xp_ = xx * this->current_cos_ -
						  yy * this->current_sin_;
		this->twist_yp_ = xx * this->current_sin_ +
						  yy * this->current_cos_;
	}
	void twist_pos(
		const double pixel_pos, int &xi, int &yi)
	{
		/* 基準半径からの比 */
		const double tt = (this->current_radius_ + pixel_pos) /
						  this->twist_radius_;
 
		/* Twist曲線上の比による位置 */
		const double xx = tt * cos(tt * this->twist_radian_);
		const double yy = tt * sin(tt * this->twist_radian_);
 
		/* その位置の角度に傾ける */
		double x2 = xx * this->current_cos_ -
					yy * this->current_sin_;
		double y2 = xx * this->current_sin_ +
					yy * this->current_cos_;
 
		/* Twist曲線上の位置の差 */
		x2 = (x2 - this->twist_xp_);
		y2 = (y2 - this->twist_yp_);
 
		/* Pixel画像上の位置の差 */
		x2 *= this->twist_radius_;
		y2 *= this->twist_radius_;
 
		/* Pixel画像上の位置 */
		x2 += this->current_xp_;
		y2 += this->current_yp_;
 
		xi = static_cast<int>(x2);
		yi = static_cast<int>(y2);
	}
 
private:
	twist_();
 
	const double xc_;
	const double yc_;
 
	double twist_xp_;
	double twist_yp_;
	double twist_radian_;
	double twist_radius_;
 
	double current_xp_;
	double current_yp_;
	double current_radian_;
	double current_cos_;
	double current_sin_;
	double current_radius_;
 
	/* copy constructorを無効化 */
	twist_(const twist_ &);
 
	/* 代入演算子を無効化 */
	twist_ &operator=(const twist_ &);
};
//--------------------------------
template <class T>
class radial_twist_
{
public:
	radial_twist_(
		const T *in_top, const int height /* 求める画像(out)の高さ */
		,
		const int width /* 求める画像(out)の幅 */
		,
		const int channels, const double xc, const double yc, const double sub_size, const int imax, const double dmax, const double intensity /* 平均値ぼかし強度 */
		,
		const double radius /* 平均値ぼかしの始まる半径 */
		,
		const double twist_radian, const double twist_radius)
		: in_top_(in_top), hh_(height), ww_(width), cc_(channels), xc_(xc), yc_(yc), sub_size_(sub_size), imax_(imax), dmax_(dmax), intensity_(intensity), radius_(radius), cl_tw_(
																																												height, xc, yc, twist_radian, twist_radius)
	{
	}
	void pixel_value(
		const T *in_current_pixel, const int xx, const int yy, const int z1, const int z2, const double ref_increase_val, const double ref_decrease_val, const double each_pixel_blur_ratio, T *result_pixel)
	{
		/* Pixel位置(0.5 1.5 2.5 ...) */
		const double xp = static_cast<double>(xx) + 0.5;
		const double yp = static_cast<double>(yy) + 0.5;
 
		/* 中心からPixel位置へのベクトルと長さ */
		const double xv = xp - this->xc_;
		const double yv = yp - this->yc_;
		const double dist = sqrt(xv * xv + yv * yv);
 
		/* 指定半径の範囲内なら何もしない */
		if (dist <= this->radius_) {
			for (int zz = z1; zz <= z2; ++zz) {
				result_pixel[zz] = in_current_pixel[zz];
			}
			return;
		}
 
		/* Radial方向のSamplingの開始位置と終了位置 */
		double scale = this->intensity_;
		if (0.0 <= each_pixel_blur_ratio) {
			scale *= each_pixel_blur_ratio;
		}
		const double length = (dist - this->radius_) *
							  scale;
		const double count_half = floor(length / 2.0 / this->sub_size_);
		const double sub_sta = -this->sub_size_ * count_half;
		const double sub_end = this->sub_size_ * count_half;
 
		/* 積算値と積算回数 */
		std::vector<double> accum_val(this->cc_);
		int accum_counter = 0;
 
		/* SubPixelによるSamplingの相対位置 */
		for (double xsub = this->sub_size_ / 2.0 - 0.5;
			 xsub < 0.5; xsub += this->sub_size_) {
			for (double ysub = this->sub_size_ / 2.0 - 0.5;
				 ysub < 0.5; ysub += this->sub_size_) {
 
				/* Sampling位置をセット */
				this->cl_tw_.current_pos(xp + xsub, yp + ysub);
 
				/* Radial&Twist方向のSampling */
				for (double tt = sub_sta; tt <= sub_end; tt += this->sub_size_) {
					/* Sampling位置からPixel位置を得る */
					int xi = 0;
					int yi = 0;
					this->cl_tw_.twist_pos(tt, xi, yi);
 
					/* clamp */
					xi = (xi < 0) ? 0 : ((this->ww_ <= xi) ? this->ww_ - 1 : xi);
					yi = (yi < 0) ? 0 : ((this->hh_ <= yi) ? this->hh_ - 1 : yi);
 
					/* 画像のPixel位置 */
					const T *in_current = this->in_top_ +
										  this->cc_ * this->ww_ * yi +
										  this->cc_ * xi;
 
					/* 積算 */
					for (int zz = z1; zz <= z2; ++zz) {
						accum_val[zz] +=
							static_cast<double>(in_current[zz]);
					}
					++accum_counter;
				}
			}
		}
		/* 積算しなかったとき(念のためのCheck) */
		if (accum_counter <= 0) {
			for (int zz = z1; zz <= z2; ++zz) {
				result_pixel[zz] = in_current_pixel[zz];
			}
			return;
		}
		/* ここで画像Pixelに保存 */
		for (int zz = z1; zz <= z2; ++zz) {
			accum_val[zz] /= static_cast<double>(accum_counter);
 
			if ((0 <= ref_increase_val) &&
				(in_current_pixel[zz] < accum_val[zz])) {
				/* 増分のみMask! */
				accum_val[zz] =
					static_cast<double>(in_current_pixel[zz]) +
					(accum_val[zz] - in_current_pixel[zz]) *
						ref_increase_val;
			} else if ((0 <= ref_decrease_val) &&
					   (accum_val[zz] < in_current_pixel[zz])) {
				/* 減分のみMask! */
				accum_val[zz] =
					static_cast<double>(in_current_pixel[zz]) +
					(accum_val[zz] - in_current_pixel[zz]) *
						ref_decrease_val;
			}
 
			accum_val[zz] += 0.5; /* 誤差対策 */
			if (this->dmax_ < accum_val[zz]) {
				result_pixel[zz] = static_cast<T>(this->imax_);
			} else if (accum_val[zz] < 0) {
				result_pixel[zz] = 0;
			} else {
				result_pixel[zz] = static_cast<T>(
					accum_val[zz]);
			}
		}
	}
 
private:
	radial_twist_() {}
 
	const T *in_top_;
	const int hh_;
	const int ww_;
	const int cc_;
	const double xc_;
	const double yc_;
	const double sub_size_;
	const int imax_;
	const double dmax_;
	const double intensity_;
	const double radius_;
	twist_ cl_tw_;
 
	/* copy constructorを無効化 */
	radial_twist_(const radial_twist_ &);
 
	/* 代入演算子を無効化 */
	radial_twist_ &operator=(const radial_twist_ &);
};
//--------------------------------
template <class IT, class RT>
void twist_convert_template_(
	const IT *in, const int margin /* 参照画像(in)がもつ余白 */
 
	,
	const RT *ref /* 求める画像(out)と同じ高さ、幅、チャンネル数 */
	,
	const int ref_bits, const int ref_mode // R,G,B,A,luminance
 
	,
	IT *out, const int hh /* 求める画像(out)の高さ */
	,
	const int ww /* 求める画像(out)の幅 */
	,
	const int cc, const double xc, const double yc, const double twist_radian, const double twist_radius, const double intensity /* 強度。ゼロより大きく2以下 */
	/* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */
	,
	const double radius /* 平均値ぼかしの始まる半径 */
	,
	const int sub_div /* 1ならJaggy、2以上はAntialias */
	,
	const bool alpha_rendering_sw)
{
	ref_bits; // for warning
 
	/* 強度のないとき、または、サブ分割がないとき、なにもしない */
	if (intensity <= 0.0 || sub_div <= 0) {
		return;
	}
 
	radial_twist_<IT> cl_ratw(
		in, hh + margin * 2, ww + margin * 2, cc, xc + margin, yc + margin, 1.0 / sub_div, std::numeric_limits<IT>::max(), static_cast<double>(std::numeric_limits<IT>::max()), intensity, radius, twist_radian, twist_radius);
 
#if defined RGBA_ORDER_OF_TOONZ6
	const int z1 = igs::image::rgba::blu;
	const int z2 = igs::image::rgba::red;
#elif defined RGBA_ORDER_OF_OPENGL
	const int z1 = igs::image::rgba::red;
	const int z2 = igs::image::rgba::blu;
#else
	Must be define / DRGBA_ORDER_OF_TOONZ6 or / DRGBA_ORDER_OF_OPENGL
#endif
 
	const IT *p_in = in + margin * (ww + margin * 2) * cc + margin * cc;
	IT *pout = out;
 
	if (0 == ref) { /* 参照なし */
		if (igs::image::rgba::siz == cc) {
			using namespace igs::image::rgba;
			if (alpha_rendering_sw) { /* Alphaも処理する */
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc) {
						cl_ratw.pixel_value(p_in, xx, yy, alp, alp, -1.0, -1.0, -1.0, pout);
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ratw.pixel_value(p_in, xx, yy, z1, z2, -1.0, -1.0, -1.0, pout);
					}
				}
			} else { /* Alpha処理しない、RGB増分をAlphaでMaskする */
				const unsigned int val_max = std::numeric_limits<IT>::max();
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc) {
						pout[alp] = p_in[alp];
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ratw.pixel_value(
							p_in, xx, yy, z1, z2, static_cast<double>(p_in[alp]) / val_max, -1.0, -1.0, pout);
					}
				}
			}
		} else { /* Alphaがない, RGB/Grayscale... */
			for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
				for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc) {
					cl_ratw.pixel_value(p_in, xx, yy, 0, cc - 1, -1.0, -1.0, -1.0, pout);
				}
			}
		}
	} else { /* 参照あり */
		const RT *refe = ref;
		const int r_max = std::numeric_limits<RT>::max();
		if (igs::image::rgba::siz == cc) {
			using namespace igs::image::rgba;
			if (alpha_rendering_sw) { /* Alphaも処理する */
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc, refe += cc) {
						const double refv = igs::color::ref_value(
							refe, cc, r_max, ref_mode);
						if (0 == refv) {
							for (int zz = 0; zz < cc; ++zz) {
								pout[zz] = p_in[zz];
							}
							continue;
						}
 
						cl_ratw.pixel_value(p_in, xx, yy, alp, alp, -1.0, -1.0, refv, pout);
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ratw.pixel_value(p_in, xx, yy, z1, z2, -1.0, -1.0, refv, pout);
					}
				}
			} else { /* Alpha処理しない、RGB増分をAlphaでMaskする */
				const unsigned int val_max = std::numeric_limits<IT>::max();
				for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
					for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc, refe += cc) {
						const double refv = igs::color::ref_value(
							refe, cc, r_max, ref_mode);
						if (0 == refv) {
							for (int zz = 0; zz < cc; ++zz) {
								pout[zz] = p_in[zz];
							}
							continue;
						}
 
						pout[alp] = p_in[alp];
						if (0 == pout[alp]) {
							pout[red] = p_in[red];
							pout[gre] = p_in[gre];
							pout[blu] = p_in[blu];
							continue;
						}
						cl_ratw.pixel_value(p_in, xx, yy, z1, z2, static_cast<double>(p_in[alp]) / val_max, -1.0, refv, pout);
					}
				}
			}
		} else { /* Alphaがない, RGB/Grayscale... */
			for (int yy = margin; yy < hh + margin; ++yy, p_in += 2 * margin * cc) {
				for (int xx = margin; xx < ww + margin; ++xx, p_in += cc, pout += cc, refe += cc) {
					const double refv = igs::color::ref_value(
						refe, cc, r_max, ref_mode);
					if (0 == refv) {
						for (int zz = 0; zz < cc; ++zz) {
							pout[zz] = p_in[zz];
						}
						continue;
					}
 
					cl_ratw.pixel_value(p_in, xx, yy, 0, cc - 1, -1.0, -1.0, refv, pout);
				}
			}
		}
	}
}
}
//--------------------------------------------------------------------
void igs::radial_blur::convert(
	const unsigned char *in, const int margin /* 参照画像(in)がもつ余白 */
 
	,
	const unsigned char *ref /* outと同じ高さ、幅、チャンネル数 */
	,
	const int ref_bits, const int ref_mode // R,G,B,A,luminance
 
	,
	unsigned char *out
 
	,
	const int height /* 求める画像(out)の高さ */
	,
	const int width /* 求める画像(out)の幅 */
	,
	const int channels, const int bits
 
	,
	const double xc, const double yc, const double twist_radian, const double twist_radius
	/*
Twist Radius
	ひねりの基準半径を指定します。
	単位はミリメートルです。
	最小は0。最大1000mm(1m)です。
	初期値は0でひねりはありません。
*/
	,
	const double intensity /* 強度。ゼロより大きく2以下 */
	/* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */
	,
	const double radius /* 平均値ぼかしの始まる半径 */
	,
	const int sub_div /* 1ならJaggy、2以上はAntialias */
	,
	const bool alpha_rendering_sw)
{
	if ((igs::image::rgba::siz != channels) &&
		(igs::image::rgb::siz != channels) &&
		(1 != channels) /* bit(monoBW) */
		) {
		throw std::domain_error("Bad channels,Not rgba/rgb/grayscale");
	}
 
	if (
		(std::numeric_limits<unsigned char>::digits != bits) &&
		(std::numeric_limits<unsigned short>::digits != bits)) {
		throw std::domain_error("Bad in bits,Not uchar/ushort");
	}
	if ((0 != ref) &&
		(std::numeric_limits<unsigned char>::digits != ref_bits) &&
		(std::numeric_limits<unsigned short>::digits != ref_bits)) {
		throw std::domain_error("Bad ref bits,Not uchar/ushort");
	}
 
	/* 強度のないとき、または、サブ分割がないとき */
	if (intensity <= 0.0 || sub_div <= 0) {
		if (std::numeric_limits<unsigned char>::digits == bits) {
			igs::image::copy_except_margin(
				in, margin, out, height, width, channels);
		} else if (std::numeric_limits<unsigned short>::digits == bits) {
			igs::image::copy_except_margin(
				reinterpret_cast<const unsigned short *>(in), margin, reinterpret_cast<unsigned short *>(out), height, width, channels);
		}
		return;
	}
 
	if (0.0 == twist_radian) {
		if (std::numeric_limits<unsigned char>::digits == ref_bits) {
			if (std::numeric_limits<unsigned char>::digits == bits) {
				radial_convert_template_(
					in, margin, ref, ref_bits, ref_mode, out, height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);
			} else if (std::numeric_limits<unsigned short>::digits == bits) {
				radial_convert_template_(
					reinterpret_cast<const unsigned short *>(in), margin, ref, ref_bits, ref_mode, reinterpret_cast<unsigned short *>(out), height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);
			}
		} else { /* ref_bitsがゼロでも(refがなくても)ここにくる */
			if (std::numeric_limits<unsigned char>::digits == bits) {
				radial_convert_template_(
					in, margin, reinterpret_cast<const unsigned short *>(ref), ref_bits, ref_mode, out, height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);
			} else if (std::numeric_limits<unsigned short>::digits == bits) {
				radial_convert_template_(
					reinterpret_cast<const unsigned short *>(in), margin, reinterpret_cast<const unsigned short *>(ref), ref_bits, ref_mode, reinterpret_cast<unsigned short *>(out), height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);
			}
		}
	} else {
		if (std::numeric_limits<unsigned char>::digits == ref_bits) {
			if (std::numeric_limits<unsigned char>::digits == bits) {
				twist_convert_template_(
					in, margin, ref, ref_bits, ref_mode, out, height, width, channels, xc, yc, twist_radian, twist_radius, intensity, radius, sub_div, alpha_rendering_sw);
			} else if (std::numeric_limits<unsigned short>::digits == bits) {
				twist_convert_template_(
					reinterpret_cast<const unsigned short *>(in), margin, ref, ref_bits, ref_mode, reinterpret_cast<unsigned short *>(out), height, width, channels, xc, yc, twist_radian, twist_radius, intensity, radius, sub_div, alpha_rendering_sw);
			}
		} else { /* ref_bitsがゼロでも(refがなくても)ここにくる */
			if (std::numeric_limits<unsigned char>::digits == bits) {
				twist_convert_template_(
					in, margin, reinterpret_cast<const unsigned short *>(ref), ref_bits, ref_mode, out, height, width, channels, xc, yc, twist_radian, twist_radius, intensity, radius, sub_div, alpha_rendering_sw);
			} else if (std::numeric_limits<unsigned short>::digits == bits) {
				twist_convert_template_(
					reinterpret_cast<const unsigned short *>(in), margin, reinterpret_cast<const unsigned short *>(ref), ref_bits, ref_mode, reinterpret_cast<unsigned short *>(out), height, width, channels, xc, yc, twist_radian, twist_radius, intensity, radius, sub_div, alpha_rendering_sw);
			}
		}
	}
}
#if 0  //-------------------- comment out start ------------------------
namespace {
 double enlarge_margin_length_(
	const double xc
	,const double yc
	,const double xp
	,const double yp
	,const double intensity
	,const double radius
	,const double sub_size
 ) {
	const double xv = xp - xc;
	const double yv = yp - yc;
	const double dist = sqrt(xv*xv + yv*yv);
 
	/* "dist = len - (len - radius) * intensity / 2.0;"より */
	const double len = (2.0 * dist - radius * intensity) /
				(2.0-intensity);
 
	/***const double length = (len - radius) * intensity;
	const double count_half = floor(length / sub_size);
	return sub_size * count_half;***/
 
	if (len <= radius) { return 0; }
	return (len - radius) * intensity / 2.0;
 }
}
int igs::radial_blur::enlarge_margin( /* Twist時は正確ではない... */
	const int height
	,const int width
	,const double xc
	,const double yc
	,const double twist_radian
	,const double twist_radius
	,const double intensity/* 強度。ゼロより大きく2以下 */
 /* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */
	,const double radius	/* 平均値ぼかしの始まる半径 */
	,const int sub_div	/* 1ならJaggy、2以上はAntialias */
) {
	/* 強度のないとき、2以上のとき、または、サブ分割がないとき、
	なにもしない */
	if (intensity <= 0.0 || 2.0 <= intensity || sub_div <= 0) {
		return 0;
	}
 
	double margin1 = 0;
	double margin2 = 0;
 
	margin1 = enlarge_margin_length_(
	  xc,yc,-width/2.0,-height/2.0,intensity,radius,1.0/sub_div);
 
	margin2 = enlarge_margin_length_(
	  xc,yc,-width/2.0, height/2.0,intensity,radius,1.0/sub_div);
	if (margin1 < margin2) { margin1 = margin2; }
 
	margin2 = enlarge_margin_length_(
	  xc,yc, width/2.0,-height/2.0,intensity,radius,1.0/sub_div);
	if (margin1 < margin2) { margin1 = margin2; }
 
	margin2 = enlarge_margin_length_(
	  xc,yc, width/2.0, height/2.0,intensity,radius,1.0/sub_div);
	if (margin1 < margin2) { margin1 = margin2; }
 
	return static_cast<int>(ceil(margin1));
}
#endif //-------------------- comment out end -------------------------
namespace
{
double reference_margin_length_(
	const double xc, const double yc, const double xp, const double yp, const double intensity, const double radius, const double sub_size)
{
	const double xv = xp - xc;
	const double yv = yp - yc;
	const double dist = sqrt(xv * xv + yv * yv);
	if (dist <= radius) {
		return 0;
	}
	const double half_length = (dist - radius) * intensity / 2.0;
	const double count_half = floor(half_length / sub_size);
	return sub_size * count_half;
}
}
int igs::radial_blur::reference_margin(																																					 /* Twist時は正確ではない... */
									   const int height, const int width, const double xc, const double yc, const double twist_radian, const double twist_radius, const double intensity /* 強度。ゼロより大きく2以下 */
									   /* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */
									   ,
									   const double radius /* 平均値ぼかしの始まる半径 */
									   ,
									   const int sub_div /* 1ならJaggy、2以上はAntialias */
									   )
{
	twist_radius; // for warning
	twist_radian; // for warning
 
	/* 強度のないとき、2以上のとき、または、サブ分割がないとき、
	なにもしない */
	if (intensity <= 0.0 || 2.0 <= intensity || sub_div <= 0) {
		return 0;
	}
 
	double margin1 = 0;
	double margin2 = 0;
 
	/*
	四隅から参照する外部への最大マージンを計算する
	*/
	margin1 = reference_margin_length_(
		xc, yc, -width / 2.0, -height / 2.0, intensity, radius, 1.0 / sub_div);
 
	margin2 = reference_margin_length_(
		xc, yc, -width / 2.0, height / 2.0, intensity, radius, 1.0 / sub_div);
	if (margin1 < margin2) {
		margin1 = margin2;
	}
 
	margin2 = reference_margin_length_(
		xc, yc, width / 2.0, -height / 2.0, intensity, radius, 1.0 / sub_div);
	if (margin1 < margin2) {
		margin1 = margin2;
	}
 
	margin2 = reference_margin_length_(
		xc, yc, width / 2.0, height / 2.0, intensity, radius, 1.0 / sub_div);
	if (margin1 < margin2) {
		margin1 = margin2;
	}
 
	return static_cast<int>(ceil(margin1));
}