#include <cmath>      // sqrt() sin() cos()
#include <vector>     // std::vector
#include <stdexcept>  // std::domain_error()

namespace {
bool inside_polygon_(double radius, int odd_diameter, double xp, double yp,
                     int polygon_num, double degree) {
  double radian = degree * (M_PI / 180), add_radian = 2.0 * M_PI / polygon_num,
         x1 = 0, y1 = 0, x2, y2, xa = -odd_diameter, xb = -odd_diameter;

  for (int ii = 0; ii <= polygon_num;
       ++ii, radian += add_radian, x1 = x2, y1 = y2) {
    /* (x2,y2)は(0,0)を原点とした正多角形の各頂点 */
    /* 注意:参照の図形と 計算結果の図形は左右反転する。
            そのため結果画像は左からCW回転となってしまう。
            数学的開始点(右)及び回転方向(CCW)とするため、
            内部で上下左右反転する。
    */
    x2 = -radius * cos(radian);
    y2 = -radius * sin(radian);

    /* 線分をみるので、loop次へ */
    if (ii <= 0) {
      continue;
    }

    /* ypのy scanline上で交差
    ((y2==y1==yp)の場合も)しない線分なら次へ */
    if (!(((y1 <= yp) && (yp <= y2)) || ((y2 <= yp) && (yp <= y1)))) {
      continue;
    }

    /* 水平線分上 */
    if (y2 == y1) {
      if (((x1 <= xp) && (xp <= x2)) || ((x2 <= xp) && (xp <= x1))) {
        return true;
      } else {
        return false;
      } /* 水平範囲外 */
    }

    /* 初めの交差 */
    if (xa == -odd_diameter) {
      /* (x2 - x1)/(y2 - y1)=(xa - x1)/(yp - y1);
(xa - x1)=(yp - y1)*(x2 - x1)/(y2 - y1); */
      xa = (yp - y1) * (x2 - x1) / (y2 - y1) + x1;
    } else
        /* 2番目の交差 */
        if (xb == -odd_diameter) {
      xb = (yp - y1) * (x2 - x1) / (y2 - y1) + x1;
      if (((xa <= xp) && (xp <= xb)) || ((xb <= xp) && (xp <= xa))) {
        return true;
      } else {
        return false;
      } /* 水平範囲外 */
    }
  }
  return false;
}
/* 円錐型(1==curve)に減衰分布する２次元配列を生成する */
// bool attenuation_distribution_(
void attenuation_distribution_(
    std::vector<std::vector<double>> &lens_matrix,
    std::vector<int> &lens_offsets, std::vector<double *> &lens_starts,
    std::vector<int> &lens_sizes, int &odd_diameter,
    const double radius /* =1.01 円半径(Pixel単位)1より大きい値 */
    ,
    const double curve /* =1. 0:無 0<.<1.0:減 1.0:リニア 1.0<:増 */
    ,
    const int polygon_num /* =2  2:円 3<=:円内接多角形(中心右始) */
    ,
    const double degree /* =0. 多角形開始角度(傾き)(反時計回り) */
    ) {
  /*
  double radius
          影響の半径
          1より大きい浮動小数値
          diameterはこの値から自動設定する
  double curve
          中心から周辺への下り具合(デフォルトは1)
          光り拡散の強弱
          zeroでない浮動小数値
*/
  /*
  lens_offsetsとlens_sizesで影響範囲を表わすmatrixを表わす
  radiusは影響円の半径
  matrix(縦横)サイズ(lens_offsets.size())は
          円(radius)が入る最小整数値でかつ、
          1以上の奇数(1,3,5,...)値

  影響は上下反転、左右反転するので注意、以下説明図
  y

  ^
  |	        |\
  |	        |   \
  |	        |  *   >
  |	|\      |   /
  |	|   \   |/
  |	|  *   >P
  |	|   /   |\
  |	|/      |   \
  |	        |  *   >
  |	        |   /
  |	        |/
--+--------------------------------------> x
  y

  ^
  |	   -----------
  |	    \   *   /
  |	      \   /
  |     ----------P----------
  |      \   *   / \   *   /
  |	 \   /     \   /
  |	   v         v
--+--------------------------------------> x
  |   P       : Target pixel
  |   <,/,\,| : 影響範囲
  |   *       : 描画の結果
*/
  /* --- 散光の力がゼロだと画像に対する変化はない ----- */
  if (0.0 == curve) {
    // return false;
    throw std::domain_error("curve is zero");
  }

  /* --- 直径は、半径(radius)の2倍より大きい整数値 ---- */
  odd_diameter = static_cast<int>(ceil(radius * 2.0));

  /* --- 直径が1Pixel以下は画像に対する変化はない ----- */
  if (odd_diameter <= 1) {
    // return false;
    throw std::domain_error("diameter is equal less than 1");
  }

  /* --- 中心pixelが必要なので、奇数(1,3,5...)にする -- */
  if (0 == (odd_diameter % 2)) {
    ++odd_diameter;
  }

  /* --- Memory確保or再利用 --------------------------- */
  lens_matrix.resize(odd_diameter);
  for (int yy = 0; yy < odd_diameter; ++yy) {
    lens_matrix.at(yy).resize(odd_diameter);
  }
  lens_offsets.resize(odd_diameter);
  lens_starts.resize(odd_diameter);
  lens_sizes.resize(odd_diameter);

  /* ---  1.円のレンズmatrix生成
          2.scanlineスタート位置とスタートポインタのセット
          3.scanlineサイズのセット -------------------- */

  /* --- scanlineスタート位置とscanlineサイズ生成 --- */
  double yp = 0.5 - (odd_diameter / 2.0); /* matrix中心からのy距離 */
  for (int yy = 0; yy < odd_diameter; ++yy, yp += 1.0) {
    lens_offsets.at(yy) = -1;               /* 初期値 */
    lens_starts.at(yy)  = 0;                /* 初期値 */
    lens_sizes.at(yy)   = -1;               /* 初期値 */
    double xp = 0.5 - (odd_diameter / 2.0); /* matrix中心からのx距離 */
    for (int xx = 0; xx < odd_diameter; ++xx, xp += 1.0) {
      const double length = sqrt(xp * xp + yp * yp);
      if ((length <= radius) &&
          ((polygon_num < 3) ||
           inside_polygon_(radius, odd_diameter, xp, yp, polygon_num,
                           degree))) { /* 影響内 */

        /* 1.円のレンズmatrix生成 */

        /* 中心が1で外郭が0その間をリニアに */
        double val = 1.0 - (length / radius);

        /* 光の減衰の様子をGamma曲線で指定する */
        /* 0<.<1.0:光減 1.0:リニア 1.0<:光増 */
        val = pow(val, 1.0 / curve);

        /* 影響内は値を設定 */
        lens_matrix.at(yy).at(xx) = val;

        /* 2.scanlineスタート位置(ポインタ)のセット */

        if (lens_offsets.at(yy) < 0) {
          /* 半径内に入った瞬間その位置を記録 */
          lens_offsets.at(yy) = xx;
          lens_starts.at(yy)  = &lens_matrix.at(yy).at(xx);
        }
      } else { /* 影響外 */
        /* 1.円のレンズmatrix生成 */

        /* 影響外はゼロ */
        lens_matrix.at(yy).at(xx) = 0.0;

        /* 3.scalineサイズのセット */

        if ((0 <= lens_offsets.at(yy)) && (lens_sizes.at(yy) < 0)) {
          /* 半径内を出た瞬間そこまでのサイズを記録 */
          lens_sizes.at(yy) = xx - lens_offsets.at(yy);
        }
      }
    }
    /* scanlineの最後までが半径内のとき */
    if ((0 <= lens_offsets.at(yy)) && (lens_sizes.at(yy) < 0)) {
      lens_sizes.at(yy) = odd_diameter - lens_offsets.at(yy);
    }
  }

  /* matrixの値を正規化 */
  double total = 0.0;
  for (unsigned int yy = 0; yy < lens_matrix.size(); ++yy) {
    for (unsigned int xx = 0; xx < lens_matrix.at(yy).size(); ++xx) {
      total += lens_matrix.at(yy).at(xx);
    }
  }
  for (unsigned int yy = 0; yy < lens_matrix.size(); ++yy) {
    for (unsigned int xx = 0; xx < lens_matrix.at(yy).size(); ++xx) {
      lens_matrix.at(yy).at(xx) /= total;
    }
  }
  // return true;
}
}