#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"

#include <qvector2d></qvector2d>

#include <qtransform></qtransform>

namespace {

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

enum Type { Accelerator = 0, Uniform_Length };

class radial_ {

  const float* in_top_;

  const int hh_;

  const int ww_;

  const int cc_;

  const TPointD center_;

  const bool antialias_sw_;

  const bool alpha_rendering_sw_;

  const double intensity_;

  const double blur_radius_;

  const double twist_radian_;

  const double pivot_radius_;

  const int type_;

  const double ellipse_aspect_ratio_;

  const double ellipse_angle_;

  const double intensity_correlation_with_ellipse_;

  QTransform tr_, tr_inv_;

public:

  radial_(const float* in_top, const int height, const int width,

          const int channels, const TPointD center, const bool antialias_sw,

          const bool alpha_rendering_sw,

          const double intensity,   /* 平均値ぼかし強度 */

          const double blur_radius, /* 平均値ぼかしの始まる半径 */

          const double twist_radian, const double pivot_radius, const int type,

          const double ellipse_aspect_ratio, const double ellipse_angle,

          const double intensity_correlation_with_ellipse)

      : in_top_(in_top)

      , hh_(height)

      , ww_(width)

      , cc_(channels)

      , center_(center)

      , antialias_sw_(antialias_sw)

      , alpha_rendering_sw_(alpha_rendering_sw)

      , intensity_(intensity)

      , blur_radius_(blur_radius)

      , twist_radian_(twist_radian)

      , pivot_radius_(pivot_radius)

      , type_(type)

      , ellipse_aspect_ratio_(ellipse_aspect_ratio)

      , ellipse_angle_(ellipse_angle)

      , intensity_correlation_with_ellipse_(

            intensity_correlation_with_ellipse) {

    if (ellipse_aspect_ratio_ != 1.0) {

      double axis_x =

          2.0 * ellipse_aspect_ratio_ / (ellipse_aspect_ratio_ + 1.0);

      double axis_y = axis_x / ellipse_aspect_ratio_;

      tr_           = QTransform()

                .rotateRadians(this->ellipse_angle_)

                .scale(axis_x, axis_y);

      tr_inv_ = QTransform(tr_).inverted();

    }

  }

  void pixel_value(const float* in_current_pixel, const int xx, const int yy,

                   const bool isRGB, const double refVal, float* result_pixel) {

    auto in_pixel = [&](int x, int y) {

      /* clamp */

      x = (x < 0) ? 0 : ((this->ww_ <= x) ? this->ww_ - 1 : x);

      y = (y < 0) ? 0 : ((this->hh_ <= y) ? this->hh_ - 1 : y);

      return this->in_top_ + this->cc_ * y * this->ww_ + this->cc_ * x;

    };

    auto interp = [&](float v1, float v2, float r) {

      return v1 * (1.f - r) + v2 * r;

    };

    auto accum_interp_in_values = [&](QPointF pos, std::vector<float>& accumP,</float>

                                      int z1, int z2, float weight) {

      int xId          = (int)std::floor(pos.x());

      float rx         = pos.x() - (float)xId;

      int yId          = (int)std::floor(pos.y());

      float ry         = pos.y() - (float)yId;

      const float* p00 = in_pixel(xId, yId);

      const float* p01 = in_pixel(xId + 1, yId);

      const float* p10 = in_pixel(xId, yId + 1);

      const float* p11 = in_pixel(xId + 1, yId + 1);

      for (int zz = z1; zz <= z2; zz++) {

        accumP[zz] += weight * interp(interp(p00[zz], p01[zz], rx),

                                      interp(p10[zz], p11[zz], rx), ry);

      }

    };

    auto accum_in_values = [&](QPointF pos, std::vector<float>& accumP, int z1,</float>

                               int z2, float weight) {

      int xId        = (int)std::round(pos.x());

      int yId        = (int)std::round(pos.y());

      const float* p = in_pixel(xId, yId);

      for (int zz = z1; zz <= z2; zz++) {

        accumP[zz] += weight * p[zz];

      }

    };

    int c1, c2;

    if (isRGB) {

#if defined RGBA_ORDER_OF_TOONZ6

      c1 = igs::image::rgba::blu;

      c2 = igs::image::rgba::red;

#elif defined RGBA_ORDER_OF_OPENGL

      c1 = igs::image::rgba::red;

      c2 = igs::image::rgba::blu;

#else

      Must be define / DRGBA_ORDER_OF_TOONZ6 or / DRGBA_ORDER_OF_OPENGL

#endif

    } else {

      c1 = igs::image::rgba::alp;

      c2 = igs::image::rgba::alp;

    }

    const QPointF center(this->center_.x, this->center_.y);

    /* Pixel位置 */

    const QPointF p(static_cast<float>(xx), static_cast<float>(yy));</float></float>

    /* 中心からPixel位置へのベクトルと長さ */

    const QVector2D v(p - center);

    const float dist = v.length();

    /* 指定半径の範囲内なら何もしない */

    bool is_in_blur_radius = false;

    if (this->ellipse_aspect_ratio_ == 1.f) {

      is_in_blur_radius = (dist <= this->blur_radius_);

    } else {

      is_in_blur_radius =

          QVector2D(this->tr_inv_.map(v.toPointF())).lengthSquared() <=

          (this->blur_radius_ * this->blur_radius_);

    }

    if (is_in_blur_radius) {

      for (int c = c1; c <= c2; ++c) {

        result_pixel[c] = in_current_pixel[c];

      }

      return;

    }

    /* 積算値と積算回数 */

    std::vector<float> accum_val(this->cc_);</float>

    float accum_counter = 0.f;

    /* Radial方向のSamplingの開始位置と終了位置 */

    float scale = this->intensity_ * refVal;

    float distanceRatio = 1.f;

    float length;

    if (this->ellipse_aspect_ratio_ == 1.) {

      if (this->type_ == Accelerator)

        length = (dist - this->blur_radius_) * scale;

      else  // Uniform_Length

        length = (this->pivot_radius_ - this->blur_radius_) * scale;

    } else {

      float dist_mod = QVector2D(this->tr_inv_.map(v.toPointF())).length();

      distanceRatio  = dist / dist_mod;

      if (this->type_ == Accelerator)

        length = (dist - this->blur_radius_ * distanceRatio) * scale *

                 std::pow(distanceRatio, intensity_correlation_with_ellipse_);

      else  // Uniform_Length

        length = (this->pivot_radius_ - this->blur_radius_) * scale *

                 std::pow(distanceRatio, intensity_correlation_with_ellipse_);

    }

    float blur_length_half = length * 0.5f;

    // sampling in both directions

    float sample_length = 0.f;

    while (sample_length < blur_length_half) {

      // compute weight

      float weight = 1.f;

      if (sample_length > blur_length_half - 1.f) {

        if (antialias_sw_)

          weight = blur_length_half - sample_length;

        else

          break;

      }

      // advance to the next sample

      sample_length += weight;

      float sample_dists[2] = {dist + sample_length, dist - sample_length};

      for (int i = 0; i < 2; i++) {

        if (sample_dists[i] < 0.f) continue;  // 中心を突き抜けた場合は無視する

        QPointF v_smpl;

        float scaleVal = sample_dists[i] / dist;

        if (this->twist_radian_ == 0.0) {

          if (this->ellipse_aspect_ratio_ == 1.f)

            v_smpl = v.toPointF() * scaleVal;

          else

            v_smpl = (this->tr_inv_ *

                      QTransform(this->tr_).scale(scaleVal, scaleVal))

                         .map(v.toPointF());

        } else {  // twisted case

          if (this->ellipse_aspect_ratio_ == 1.f) {

            double tmp_twist_angle =

                twist_radian_ * (sample_dists[i] - dist) / pivot_radius_;

            float cos = std::cos(tmp_twist_angle);

            float sin = std::sin(tmp_twist_angle);

            v_smpl =

                QPointF(cos * v.x() - sin * v.y(), sin * v.x() + cos * v.y()) *

                scaleVal;

          } else {

            double tmp_twist_angle = twist_radian_ * (sample_dists[i] - dist) *

                                     distanceRatio / pivot_radius_;

            v_smpl = (this->tr_inv_ * QTransform(this->tr_)

                                          .rotateRadians(tmp_twist_angle)

                                          .scale(scaleVal, scaleVal))

                         .map(v.toPointF());

          }

        }

        if (antialias_sw_)

          accum_interp_in_values(v_smpl + center, accum_val, c1, c2, weight);

        else

          accum_in_values(v_smpl + center, accum_val, c1, c2, weight);

        accum_counter += weight;

      }

    }

    // sample the original pos

    if (antialias_sw_)

      accum_interp_in_values(p, accum_val, c1, c2, 1.f);

    else

      accum_in_values(p, accum_val, c1, c2, 1.f);

    accum_counter += 1.f;

    /* 積算しなかったとき(念のためのCheck) */

    if (accum_counter <= 0.f) {

      for (int c = c1; c <= c2; ++c) {

        result_pixel[c] = in_current_pixel[c];

      }

      return;

    }

    /* ここで画像Pixelに保存 */

    for (int c = c1; c <= c2; ++c) {

      accum_val[c] /= accum_counter;

      if (isRGB && !this->alpha_rendering_sw_ &&

          (in_current_pixel[c] < accum_val[c]) &&

          result_pixel[igs::image::rgba::alp] < 1.f) {

        result_pixel[c] =

            in_current_pixel[c] + (accum_val[c] - in_current_pixel[c]) *

                                      result_pixel[igs::image::rgba::alp];

      } else {

        result_pixel[c] = accum_val[c];

      }

    }

  }

private:

  radial_();

  /* copy constructorを無効化 */

  radial_(const radial_&);

  /* 代入演算子を無効化 */

  radial_& operator=(const radial_&);

};

void radial_convert(

    const float* in, float* out, const int margin, /* 参照画像(in)がもつ余白 */

    const TDimension out_dim,             /* 求める画像(out)のサイズ */

    const int channels, const float* ref, /* 求める画像(out)と同じ高さ、幅 */

    const TPointD center,

    const double intensity, /* 強度。ゼロより大きく2以下 */

    /* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */

    const double radius, /* 平均値ぼかしの始まる半径 */

    const double twist_radian, const double pivot_radius,

    const int type,          // 0: Accelerator, 1: Uniform Length

    const bool antialias_sw, /*  when true, sampled pixel will be

                                bilinear-interpolated */

    const bool alpha_rendering_sw, const double ellipse_aspect_ratio,

    const double ellipse_angle,

    const double intensity_correlation_with_ellipse) {

  assert(intensity > 0.0);

  radial_ radial(in, out_dim.ly + margin * 2, out_dim.lx + margin * 2, channels,

                 center, antialias_sw, alpha_rendering_sw, intensity, radius,

                 twist_radian, pivot_radius, type, ellipse_aspect_ratio,

                 ellipse_angle * M_PI_180, intensity_correlation_with_ellipse);

  const float* p_in =

      in + margin * (out_dim.lx + margin * 2) * channels + margin * channels;

  float* p_out       = out;

  const float* p_ref = ref;

  for (int yy = margin; yy < out_dim.ly + margin;

       ++yy, p_in += 2 * margin * channels) {

    for (int xx = margin; xx < out_dim.lx + margin;

         ++xx, p_in += channels, p_out += channels) {

      using namespace igs::image::rgba;

      float refVal = (ref) ? *p_ref : 1.f;

      // if the reference value is zero

      if (refVal == 0.f) {

        for (int c = 0; c < channels; ++c) p_out[c] = p_in[c];

        p_ref++;

        continue;

      }

      if (alpha_rendering_sw) {  // blur alpha

        radial.pixel_value(p_in, xx, yy, false, refVal, p_out);

      } else {  // use the src alpha as-is

        p_out[alp] = p_in[alp];

      }

      if (p_out[alp] == 0.f) {

        p_out[red] = p_in[red];

        p_out[gre] = p_in[gre];

        p_out[blu] = p_in[blu];

        if (ref) p_ref++;

        continue;

      }

      // blur RGB channels

      radial.pixel_value(p_in, xx, yy, true, refVal, p_out);

      if (ref) p_ref++;

    }

  }

}

}  // namespace

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

void igs::radial_blur::convert(

    const float* in, float* out, const int margin, /* 参照画像(in)がもつ余白 */

    const TDimension out_dim,             /* 求める画像(out)の大きさ*/

    const int channels, const float* ref, /* outと同じ高さ、幅 */

    const TPointD center, const double twist_radian, const double twist_radius,

    const double intensity, /* 強度。ゼロより大きく2以下 */

    /* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい */

    const double radius, /* 平均値ぼかしの始まる半径 */

    const int type, const bool antialias_sw, const bool alpha_rendering_sw,

    const double ellipse_aspect_ratio, const double ellipse_angle,

    const double intensity_correlation_with_ellipse) {

  /* 強度のないとき */

  if (intensity <= 0.0) {

    igs::image::copy_except_margin(in, margin, out, out_dim.ly, out_dim.lx,

                                   channels);

    return;

  }

  radial_convert(in, out, margin, out_dim, channels, ref, center, intensity,

                 radius, twist_radian, twist_radius, type, antialias_sw,

                 alpha_rendering_sw, ellipse_aspect_ratio, ellipse_angle,

                 intensity_correlation_with_ellipse);

}

//-------------------- comment out end -------------------------

namespace {

double reference_margin_length_(const TPointD center, const double xp,

                                const double yp, const double intensity,

                                const double radius, const double pivot_radius,

                                const int type,

                                const double ellipse_aspect_ratio,

                                const double intensity_correlation_with_ellipse,

                                const QTransform& tr_inv) {

  const QPointF c(center.x, center.y);

  const QPointF p(xp, yp);

  const QVector2D v(p - c);

  const float dist = v.length();

  if (dist <= radius) {

    return 0;

  }

  double length;

  if (ellipse_aspect_ratio == 1.) {

    if (type == Accelerator)

      length = (dist - radius) * intensity;

    else

      length = (pivot_radius - radius) * intensity;

  } else {

    float dist_mod       = QVector2D(tr_inv.map(v.toPointF())).length();

    double distanceRatio = dist / dist_mod;

    if (type == Accelerator)

      length = (dist - radius * distanceRatio) * intensity *

               std::pow(distanceRatio, intensity_correlation_with_ellipse);

    else  // Uniform_Length

      length = (pivot_radius - radius) * intensity *

               std::pow(distanceRatio, intensity_correlation_with_ellipse);

  }

  return length;

}

}  // namespace

int igs::radial_blur::

    reference_margin(/* Twist時は正確ではない... */

                     const int height, const int width, const TPointD center,

                     const double twist_radian, const double twist_radius,

                     const double intensity, /* 強度。ゼロより大きく2以下 */

                     /* radius円境界での平均値ぼかしゼロとするためintensityは2より小さい

                      */

                     const double radius, /* 平均値ぼかしの始まる半径 */

                     const int type, const double ellipse_aspect_ratio,

                     const double ellipse_angle,

                     const double intensity_correlation_with_ellipse) {

  /* 強度のないとき、2以上のとき、なにもしない */

  if (intensity <= 0.0 || 2.0 <= intensity) {

    return 0;

  }

  double margin1 = 0;

  double margin2 = 0;

  QTransform tr_inv;

  if (ellipse_aspect_ratio != 1.0) {

    double axis_x = 2.0 * ellipse_aspect_ratio / (ellipse_aspect_ratio + 1.0);

    double axis_y = axis_x / ellipse_aspect_ratio;

    tr_inv        = QTransform()

                 .rotateRadians(ellipse_angle)

                 .scale(axis_x, axis_y)

                 .inverted();

  }

  /*

  四隅から参照する外部への最大マージンを計算する

  */

  margin1 = reference_margin_length_(

      center, -width / 2.0, -height / 2.0, intensity, radius, twist_radius,

      type, ellipse_aspect_ratio, intensity_correlation_with_ellipse, tr_inv);

  margin2 = reference_margin_length_(

      center, -width / 2.0, height / 2.0, intensity, radius, twist_radius, type,

      ellipse_aspect_ratio, intensity_correlation_with_ellipse, tr_inv);

  if (margin1 < margin2) {

    margin1 = margin2;

  }

  margin2 = reference_margin_length_(

      center, width / 2.0, -height / 2.0, intensity, radius, twist_radius, type,

      ellipse_aspect_ratio, intensity_correlation_with_ellipse, tr_inv);

  if (margin1 < margin2) {

    margin1 = margin2;

  }

  margin2 = reference_margin_length_(

      center, width / 2.0, height / 2.0, intensity, radius, twist_radius, type,

      ellipse_aspect_ratio, intensity_correlation_with_ellipse, tr_inv);

  if (margin1 < margin2) {

    margin1 = margin2;

  }

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

}