#include <iostream>   // std::cout
#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_rotate_blur.h"
namespace {
//------------------------------------------------------------------
template <class T>
class rotate_ {
public:
  rotate_(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 radian, const double blur_radius,
          const double spin_radius)
      : in_top_(in_top)
      , hh_(height)
      , ww_(width)
      , cc_(channels)
      , xc_(xc)
      , yc_(yc)
      , sub_size_(sub_size)
      , imax_(imax)
      , dmax_(dmax)
      , radian_(radian)
      , blur_radius_(blur_radius)
      , spin_radius_(spin_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->blur_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;

    /* 円周接線方向Sampling１つずらした位置への回転角度 */
    const double xw = xv - this->sub_size_ * sinval;
    const double yw = yv + this->sub_size_ * cosval;
    const double sub_radian =
        acos((xv * xw + yv * yw) /
             (sqrt(xv * xv + yv * yv) * sqrt(xw * xw + yw * yw)));

    /* 積算値と積算回数 */
    std::vector<double> accum_val(this->cc_);
    int accum_counter = 0;

    /* 参照画像による強弱付加 */
    double scale = this->radian_;
    if (0.0 <= each_pixel_blur_ratio) {
      scale *= each_pixel_blur_ratio;
    }

    /* Radial方向Sampling */
    for (double ss = this->sub_size_ / 2.0 - 0.5; ss < 0.5;
         ss += this->sub_size_) {
      /* Radial方向位置 */
      const double xp2 = xp + ss * cosval;
      const double yp2 = yp + ss * sinval;

      /* ぼかす回転角度範囲range */
      double radian = scale;

      if (0.0 < this->spin_radius_) { /* 外への強調 */
        const double xv2   = xp2 - this->xc_;
        const double yv2   = yp2 - this->yc_;
        const double dist2 = sqrt(xv2 * xv2 + yv2 * yv2);
        radian *=
            /**(dist2              - this->blur_radius_) /
(this->spin_radius_ - this->blur_radius_);**/
            (dist2 - this->blur_radius_) / this->spin_radius_;
      }

      /* Sampling回数 */
      int sub_count = static_cast<int>(radian / sub_radian);

      /* Sampling開始角度 */
      const double sub_start =
          (radian - sub_count * sub_radian) / 2.0 - radian / 2.0;

      /* 円周方向Sampling */
      for (double rr = sub_start; 0 < sub_count--; rr += sub_radian) {
        /* Radial位置(xp2,yp2)から円周位置(xx,yy)へ */
        const double cosval2 = cos(rr);
        const double sinval2 = sin(rr);

        const double xd = xp2 - this->xc_;
        const double yd = yp2 - this->yc_;

        int xx = static_cast<int>(xd * cosval2 - yd * sinval2 + this->xc_);
        int yy = static_cast<int>(xd * sinval2 + yd * cosval2 + this->yc_);

        /* clamp */
        xx = (xx < 0) ? 0 : ((this->ww_ <= xx) ? this->ww_ - 1 : xx);
        yy = (yy < 0) ? 0 : ((this->hh_ <= yy) ? this->hh_ - 1 : yy);

        /* 画像のPixel位置 */
        const T *in_current =
            this->in_top_ + this->cc_ * yy * this->ww_ + this->cc_ * xx;

        /* 積算 */
        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:
  rotate_() {}

  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 radian_;
  const double blur_radius_;
  const double spin_radius_;

  /* copy constructorを無効化 */
  rotate_(const rotate_ &);

  /* 代入演算子を無効化 */
  rotate_ &operator=(const rotate_ &);
};
//------------------------------------------------------------------
const double pi = 3.14159265358979;
template <class IT, class RT>
void rotate_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 degree,
    const double blur_radius /* ぼかしの始まる半径 */
    ,
    const double spin_radius /* ゼロ以上でspin指定となり、
                            かつぼかし強弱の一定になる半径となる */
    ,
    const int sub_div /* 1ならJaggy、2以上はAntialias */
    ,
    const bool alpha_rendering_sw) {
  ref_bits;  // for warning

  /* 強度のないとき、または、サブ分割がないとき、なにもしない */
  if (degree <= 0.0 || sub_div <= 0
      // || spin_radius<=blur_radius
      ) {
    return;
  }

  rotate_<IT> cl_rota(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()),
                      degree * pi / 180.0, blur_radius, spin_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_rota.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_rota.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_rota.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_rota.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_rota.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_rota.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_rota.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_rota.pixel_value(p_in, xx, yy, 0, cc - 1, -1.0, -1.0, refv, pout);
        }
      }
    }
  }
}
//------------------------------------------------------------------
}

void igs::rotate_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 degree /* ぼかしの回転角度 */
    ,
    const double blur_radius /* ぼかしの始まる半径 */
    ,
    const double spin_radius /* ゼロ以上でspin指定となり、
                            かつぼかし強弱の一定になる半径となる */
    ,
    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 (degree <= 0.0 || sub_div <= 0
      // || spin_radius<=blur_radius
      ) {
    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 (std::numeric_limits<unsigned char>::digits == ref_bits) {
    if (std::numeric_limits<unsigned char>::digits == bits) {
      rotate_convert_template_(in, margin, ref, ref_bits, ref_mode, out, height,
                               width, channels, xc, yc, degree, blur_radius,
                               spin_radius, sub_div, alpha_rendering_sw);
    } else if (std::numeric_limits<unsigned short>::digits == bits) {
      rotate_convert_template_(reinterpret_cast<const unsigned short *>(in),
                               margin, ref, ref_bits, ref_mode,
                               reinterpret_cast<unsigned short *>(out), height,
                               width, channels, xc, yc, degree, blur_radius,
                               spin_radius, sub_div, alpha_rendering_sw);
    }
  } else { /* ref_bitsがゼロでも(refがなくても)ここにくる */
    if (std::numeric_limits<unsigned char>::digits == bits) {
      rotate_convert_template_(
          in, margin, reinterpret_cast<const unsigned short *>(ref), ref_bits,
          ref_mode, out, height, width, channels, xc, yc, degree, blur_radius,
          spin_radius, sub_div, alpha_rendering_sw);
    } else if (std::numeric_limits<unsigned short>::digits == bits) {
      rotate_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, degree, blur_radius, spin_radius, sub_div, alpha_rendering_sw);
    }
  }
}
//--------------------------------------------------------------------
namespace {
double reference_margin_length_(const double xc, const double yc,
                                const double xp, const double yp, double radian,
                                const double blur_radius,
                                const double spin_radius,
                                const double sub_size) {
  sub_size;  // for warning

  const double xv = xp - xc;
  const double yv = yp - yc;

  if (0.0 < spin_radius) { /* 外への強調 */
    const double dist = sqrt(xv * xv + yv * yv);
    // radian *= (dist-blur_radius)/(spin_radius-blur_radius);
    radian *= (dist - blur_radius) / spin_radius;
  }

  double cosval      = cos(radian / 2.0);
  double sinval      = sin(radian / 2.0);
  const double x1    = xv * cosval - yv * sinval + xc;
  const double y1    = xv * sinval + yv * cosval + yc;
  const double xv1   = x1 - xp;
  const double yv1   = y1 - yp;
  const double dist1 = sqrt(xv1 * xv1 + yv1 * yv1); /* 必ずプラス値 */

  cosval             = cos(-radian / 2.0);
  sinval             = sin(-radian / 2.0);
  const double xp2   = xv * cosval - yv * sinval + xc;
  const double yp2   = xv * sinval + yv * cosval + yc;
  const double xv2   = xp2 - xp;
  const double yv2   = yp2 - yp;
  const double dist2 = sqrt(xv2 * xv2 + yv2 * yv2); /* 必ずプラス値 */

  return (dist1 < dist2) ? dist2 : dist1;
}
}
int igs::rotate_blur::reference_margin(
    const int height /* 求める画像(out)の高さ */
    ,
    const int width /* 求める画像(out)の幅 */
    ,
    const double xc, const double yc, const double degree /* ぼかしの回転角度 */
    ,
    const double blur_radius /* ぼかしの始まる半径 */
    ,
    const double spin_radius /* ゼロ以上でspin指定となり、
                            かつぼかし強弱の一定になる半径となる */
    ,
    const int sub_div /* 1ならJaggy、2以上はAntialias */
    ) {
  /* 強度のないとき、または、サブ分割がないとき、なにもしない */
  if (degree <= 0.0 || sub_div <= 0
      // || spin_radius<=blur_radius
      ) {
    return 0;
  }

  double margin1 = 0;
  double margin2 = 0;

  double deg = degree;
  if (180.0 < deg) {
    deg = 180.0;
  }

  margin1 = reference_margin_length_(xc, yc, -width / 2.0, -height / 2.0,
                                     deg * pi / 180.0, blur_radius, spin_radius,
                                     1.0 / sub_div);

  margin2 = reference_margin_length_(xc, yc, -width / 2.0, height / 2.0,
                                     deg * pi / 180.0, blur_radius, spin_radius,
                                     1.0 / sub_div);
  if (margin1 < margin2) {
    margin1 = margin2;
  }

  margin2 = reference_margin_length_(xc, yc, width / 2.0, -height / 2.0,
                                     deg * pi / 180.0, blur_radius, spin_radius,
                                     1.0 / sub_div);
  if (margin1 < margin2) {
    margin1 = margin2;
  }

  margin2 = reference_margin_length_(xc, yc, width / 2.0, height / 2.0,
                                     deg * pi / 180.0, blur_radius, spin_radius,
                                     1.0 / sub_div);
  if (margin1 < margin2) {
    margin1 = margin2;
  }

  return static_cast<int>(ceil(margin1));
}