#include <vector>	// std::vector</vector>

#include <cmath>	 // cos(),sin(),sqrt()</cmath>

#include <limits>	// std::numeric_limits<t></t></limits>

#include <stdexcept> // std::domain_error()</stdexcept>

#include "igs_ifx_common.h"

#include "igs_radial_blur.h"

namespace

{

//------------------------------------------------------------------

template <class t=""></class>

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;</double>

		const double yp = static_cast<double>(yy) + 0.5;</double>

		/* 中心から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_);</double>

		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>

				int yi = static_cast<int>(yps + tt * sinvals);</int>

				/* 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]);</double>

				}

				++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);</double>

			if ((0 <= ref_increase_val) &&

				(in_current_pixel[zz] < accum_val[zz])) {

				/* 増分のみMask! */

				accum_val[zz] =

					static_cast<double>(in_current_pixel[zz]) +</double>

					(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]) +</double>

					(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_);</t>

			} else if (accum_val[zz] < 0) {

				result_pixel[zz] = 0;

			} else {

				result_pixel[zz] = static_cast<t>(accum_val[zz]);</t>

			}

		}

	}

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 class="" it,="" rt=""></class>

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(</it>

		in, hh + margin * 2, ww + margin * 2, cc, xc + margin, yc + margin, 1.0 / sub_div, std::numeric_limits<it>::max() /* サンプリング値 */</it>

		,

		static_cast<double>(std::numeric_limits<it>::max()), intensity, radius);</it></double>

#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();</it>

				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);</double>

					}

				}

			}

		} 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();</rt>

		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();</it>

				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);</double>

					}

				}

			}

		} 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);</int>

		yi = static_cast<int>(y2);</int>

	}

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>

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;</double>

		const double yp = static_cast<double>(yy) + 0.5;</double>

		/* 中心から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_);</double>

		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]);</double>

					}

					++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);</double>

			if ((0 <= ref_increase_val) &&

				(in_current_pixel[zz] < accum_val[zz])) {

				/* 増分のみMask! */

				accum_val[zz] =

					static_cast<double>(in_current_pixel[zz]) +</double>

					(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]) +</double>

					(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_);</t>

			} else if (accum_val[zz] < 0) {

				result_pixel[zz] = 0;

			} else {

				result_pixel[zz] = static_cast<t>(</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 class="" it,="" rt=""></class>

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(</it>

		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);</it></double></it>

#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();</it>

				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);</double>

					}

				}

			}

		} 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();</rt>

		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();</it>

				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);</double>

					}

				}

			}

		} 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) &&</unsigned>

		(std::numeric_limits<unsigned short="">::digits != bits)) {</unsigned>

		throw std::domain_error("Bad in bits,Not uchar/ushort");

	}

	if ((0 != ref) &&

		(std::numeric_limits<unsigned char="">::digits != ref_bits) &&</unsigned>

		(std::numeric_limits<unsigned short="">::digits != ref_bits)) {</unsigned>

		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) {</unsigned>

			igs::image::copy_except_margin(

				in, margin, out, height, width, channels);

		} else if (std::numeric_limits<unsigned short="">::digits == bits) {</unsigned>

			igs::image::copy_except_margin(

				reinterpret_cast<const *="" short="" unsigned="">(in), margin, reinterpret_cast<unsigned *="" short="">(out), height, width, channels);</unsigned></const>

		}

		return;

	}

	if (0.0 == twist_radian) {

		if (std::numeric_limits<unsigned char="">::digits == ref_bits) {</unsigned>

			if (std::numeric_limits<unsigned char="">::digits == bits) {</unsigned>

				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) {</unsigned>

				radial_convert_template_(

					reinterpret_cast<const *="" short="" unsigned="">(in), margin, ref, ref_bits, ref_mode, reinterpret_cast<unsigned *="" short="">(out), height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);</unsigned></const>

			}

		} else { /* ref_bitsがゼロでも(refがなくても)ここにくる */

			if (std::numeric_limits<unsigned char="">::digits == bits) {</unsigned>

				radial_convert_template_(

					in, margin, reinterpret_cast<const *="" short="" unsigned="">(ref), ref_bits, ref_mode, out, height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);</const>

			} else if (std::numeric_limits<unsigned short="">::digits == bits) {</unsigned>

				radial_convert_template_(

					reinterpret_cast<const *="" short="" unsigned="">(in), margin, reinterpret_cast<const *="" short="" unsigned="">(ref), ref_bits, ref_mode, reinterpret_cast<unsigned *="" short="">(out), height, width, channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);</unsigned></const></const>

			}

		}

	} else {

		if (std::numeric_limits<unsigned char="">::digits == ref_bits) {</unsigned>

			if (std::numeric_limits<unsigned char="">::digits == bits) {</unsigned>

				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) {</unsigned>

				twist_convert_template_(

					reinterpret_cast<const *="" short="" unsigned="">(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);</unsigned></const>

			}

		} else { /* ref_bitsがゼロでも(refがなくても)ここにくる */

			if (std::numeric_limits<unsigned char="">::digits == bits) {</unsigned>

				twist_convert_template_(

					in, margin, reinterpret_cast<const *="" short="" unsigned="">(ref), ref_bits, ref_mode, out, height, width, channels, xc, yc, twist_radian, twist_radius, intensity, radius, sub_div, alpha_rendering_sw);</const>

			} else if (std::numeric_limits<unsigned short="">::digits == bits) {</unsigned>

				twist_convert_template_(

					reinterpret_cast<const *="" short="" unsigned="">(in), margin, reinterpret_cast<const *="" short="" unsigned="">(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);</unsigned></const></const>

			}

		}

	}

}

#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) {