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

namespace
{
#ifndef M_PI /* for vc2005 */
#define M_PI 3.14159265358979323846
#endif
bool inside_polygon_(
	double radius, int odd_diameter, double xp, double yp, int polygon_num, double degree)
{
	double radian = degree * M_PI / 180.0,
		   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;
}
}