#include <cmath>	 // sqrt() sin() cos()</cmath>

#include <vector>	// std::vector</vector>

#include <stdexcept> // std::domain_error()</stdexcept>

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より大きい値 */</int></double></int></std::vector<double>

	,

	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));</int>

	/* --- 直径が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;

}

}