#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,</double>

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

                              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);</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>(accum_val[zz]);</t>

      }

    }

  }

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(in, hh + margin * 2, ww + margin * 2, cc,</it>

                            xc + margin, yc + margin, 1.0 / sub_div,

                            std::numeric_limits<it>::max(),</it>

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

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

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

                                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) &&</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,</const>

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

    }

    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,</const>

            ref_mode, reinterpret_cast<unsigned *="" short="">(out), height, width,</unsigned>

            channels, xc, yc, intensity, radius, sub_div, alpha_rendering_sw);

      }

    } 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,</const>

            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,</const>

            reinterpret_cast<const *="" short="" unsigned="">(ref), ref_bits, ref_mode,</const>

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

            xc, yc, intensity, radius, sub_div, alpha_rendering_sw);

      }

    }

  } 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,</const>

            ref_mode, reinterpret_cast<unsigned *="" short="">(out), height, width,</unsigned>

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

        twist_convert_template_(

            in, margin, reinterpret_cast<const *="" short="" unsigned="">(ref), ref_bits,</const>

            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,</const>

            reinterpret_cast<const *="" short="" unsigned="">(ref), ref_bits, ref_mode,</const>

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

            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