/*
 * Polygon_RayTrace
 * Copylight (C) 2013 mocchi
 * mocchi_2003@yahoo.co.jp
 * License: Boost ver.1
 */

#include "PhisicalProperties.h"

#include <map>
#include <algorithm>

#include "opennurbs.h"

#include "rapidjson/rapidjson.h"
#include "rapidjson/document.h"
#include "rj_wrapper.h"

#include "MonteCarlo.h"

static float get_ieee754(uint8_t p[4]){
	return *reinterpret_cast<float *>(p);
}

/// ========== LightSources ==========
// 光源データ
struct LightSourceImpl{
	bool valid;
	virtual ~LightSourceImpl(){
	}
	virtual int64_t Count() const = 0;
	virtual bool Get(int64_t idx, int &blockseed, ON_3dRay &ray, double &flux, double &wavelength) const = 0;
};

struct LightSources::Impl{
	int64_t count;
	int cur_ls;
	ON_SimpleArray<int64_t> accum_count;
	ON_SimpleArray<LightSourceImpl *> impls;
	virtual ~Impl(){
		for (int i = 0; i < impls.Count(); ++i) delete impls[i];
		impls.Destroy();
	}
};

struct ParallelRayImpl : public LightSourceImpl{
	int64_t count;
	double wavelength, total_flux;

	enum Shape{
		Circle
	}shape;
	ON_Plane pln;
	double radius;
	ParallelRayImpl(rj::Value &lsv){
		valid = false;
		rj::Value &origv = lsv["origin"];
		rj::Value &dirv = lsv["direction"];
		if (!origv.IsArray() || origv.Size() != 3 || !dirv.IsArray() || dirv.Size() != 3) return;
		for (rj::SizeType j = 0; j < 3; ++j){
			if (!origv[j].IsNumber() || !dirv[j].IsNumber()) return;
		}

		pln.CreateFromNormal(
			ON_3dPoint(rj_getdouble(origv[0]), rj_getdouble(origv[1]), rj_getdouble(origv[2])),
			ON_3dVector(rj_getdouble(dirv[0]), rj_getdouble(dirv[1]), rj_getdouble(dirv[2]))
		);

		count = rj_getint64(lsv["num_ray"]);
		if (count < 0) return;

		wavelength = rj_getdouble(lsv["wavelength"]);
		total_flux = rj_getdouble(lsv["total_flux"]);
		rj::Value &shapev = lsv["shape"];
		radius = 0;
		if (shapev.IsObject()){
			// とりあえず円決め打ち
			radius = rj_getdouble(shapev["radius"]);
		}
		if (radius <= 0){
			valid = false; return;
		}
		shape = Circle;

		valid = true;
	}
	int64_t Count() const{
		return count;
	}
	bool Get(int64_t idx, int &blockseed, ON_3dRay &ray, double &flux, double &wavelength) const{
		double x, y;
		do{
			x = (CalcHaltonSequence(10 + blockseed, 2) - 0.5) * 2.0;
			y = (CalcHaltonSequence(10 + blockseed, 3) - 0.5) * 2.0;
			++blockseed;
		}while(x * x + y * y > 1.0);
		x *= radius, y *= radius;
		ray.m_P = pln.PointAt(x, y);
		ray.m_V = pln.zaxis;

		flux = total_flux / static_cast<double>(count);
		wavelength = this->wavelength;

		return true;
	}
};

struct OsramRayImpl : public LightSourceImpl{
	int64_t count, count_infile;
	double wavelength, total_flux;
	struct Item{
		ON_3dRay ray;
		double flux, wavelength;
	};
	ON_ClassArray<Item> raies;
	ON_Plane axis;
	OsramRayImpl(rj::Value &lsv){
		valid = false;

		// 原点座標の取得
		rj::Value &origv = rj_getmember(lsv, "origin");
		if (!origv.IsArray() || origv.Size() != 3) return;
		ON_3dPoint origin(rj_getdouble(origv[0]), rj_getdouble(origv[1]), rj_getdouble(origv[2]));

		// 方向の取得
		// "x_dir", "y_dir", "z_dir" のうち、2つの方向を指定していること
		rj::Value &xdirv = rj_getmember(lsv, "x_dir");
		rj::Value &ydirv = rj_getmember(lsv, "y_dir");
		rj::Value &zdirv = rj_getmember(lsv, "z_dir");

		int dir_cnt = 0;
		int dir_exist = 0;
		if (xdirv.IsArray() && xdirv.Size() == 3) ++dir_cnt, dir_exist |= 1;
		if (ydirv.IsArray() && ydirv.Size() == 3) ++dir_cnt, dir_exist |= 2;
		if (zdirv.IsArray() && zdirv.Size() == 3) ++dir_cnt, dir_exist |= 4;
		if (dir_cnt != 2) return;

		// optical の方向のベクトルを指定していること
		const char *opt = rj_getstring(rj_getmember(lsv, "optical"));
		int dir_opt = 0;
		if ((*opt == 'x' && (dir_opt = 1)) && (dir_exist & 1) == 0) return;
		if ((*opt == 'y' && (dir_opt = 2)) && (dir_exist & 2) == 0) return;
		if ((*opt == 'z' && (dir_opt = 4)) && (dir_exist & 4) == 0) return;
		if (dir_opt == 0) return;

		ON_3dVector xdir, ydir, zdir;
		if (dir_exist & 1) xdir.Set(rj_getdouble(xdirv[0]), rj_getdouble(xdirv[1]), rj_getdouble(xdirv[2]));
		if (dir_exist & 2) ydir.Set(rj_getdouble(ydirv[0]), rj_getdouble(ydirv[1]), rj_getdouble(ydirv[2]));
		if (dir_exist & 4) zdir.Set(rj_getdouble(zdirv[0]), rj_getdouble(zdirv[1]), rj_getdouble(zdirv[2]));

		// 残りの方向の計算、および opticalではない方向を垂直にする計算
		if (dir_opt == 1){
			if (dir_exist & 2){
				zdir = ON_CrossProduct(xdir, ydir);
				ydir = ON_CrossProduct(zdir, xdir);
			}else if (dir_exist & 4){
				ydir = ON_CrossProduct(zdir, xdir);
				zdir = ON_CrossProduct(xdir, ydir);
			}
		}else if (dir_opt == 2){
			if (dir_exist & 4){
				xdir = ON_CrossProduct(ydir, zdir);
				zdir = ON_CrossProduct(xdir, ydir);
			}else if (dir_exist & 1){
				zdir = ON_CrossProduct(xdir, ydir);
				xdir = ON_CrossProduct(ydir, zdir);
			}
		}else if (dir_opt == 4){
			if (dir_exist & 1){
				ydir = ON_CrossProduct(zdir, xdir);
				xdir = ON_CrossProduct(ydir, zdir);
			}else if (dir_exist & 2){
				xdir = ON_CrossProduct(ydir, zdir);
				ydir = ON_CrossProduct(zdir, xdir);
			}
		}
		xdir.Unitize(), ydir.Unitize(), zdir.Unitize();
		axis.CreateFromFrame(origin, xdir, ydir);

		count = rj_getint64(lsv["num_ray"]);
		if (count < 0) return;

		total_flux = rj_getdouble(lsv["total_flux"]);

		const char *rayfile = rj_getstring(lsv["path"]);
		FILE *fp = std::fopen(rayfile, "rb");
		if (!fp) return;
		std::fseek(fp, 0, SEEK_END);
		size_t sz = std::ftell(fp);
		std::rewind(fp);
		count_infile = 0;
		double flux_raies = 0;
		if (sz > 320){
			count_infile = (sz - 320) / 32;
			raies.SetCapacity(count_infile);
			if (count > count_infile) count = count_infile;

			for (int64_t i = 0; i < count; ++i){
				int64_t idx = static_cast<int64_t>(
					CalcHaltonSequence(i, 9) * static_cast<double>(count_infile)
				);
				if (idx >= count_infile) idx = count_infile-1;
				std::fseek(fp, idx * 32 + 320, SEEK_SET);
				uint8_t buf[32];
				std::fread(buf, 1, 32, fp);
				Item &item = raies.AppendNew();
				item.ray.m_P = axis.PointAt(get_ieee754(buf   ), get_ieee754(buf+ 4), get_ieee754(buf+ 8));
				item.ray.m_V = axis.xaxis * get_ieee754(buf+12) + axis.yaxis * get_ieee754(buf+16) + axis.zaxis * get_ieee754(buf+20);
				item.flux = get_ieee754(buf+24);
				item.wavelength = get_ieee754(buf+28);
				flux_raies += item.flux;
			}
			// raies の flux の合計値が total_flux になるように比率を変更
			double r = total_flux / flux_raies;
			for (int i = 0; i < raies.Count(); ++i) raies[i].flux *= r;
		}
		std::fclose(fp);
		if (count_infile <= 0) return;

		valid = true;
	}
	int64_t Count() const{
		return count;
	}
	bool Get(int64_t idx, int &blockseed, ON_3dRay &ray, double &flux, double &wavelength) const{
		if (idx < 0 || idx >= raies.Count()) return false;
		ray = raies[idx].ray;
		flux = raies[idx].flux;
		wavelength = raies[idx].wavelength;
		return true;
	}
};

LightSources::LightSources(rj::Value &lsv){
	pimpl = new Impl();
	pimpl->count = 0;
	pimpl->cur_ls = 0;
	pimpl->accum_count.Destroy();
	pimpl->accum_count.Append(0);
	if (!lsv.IsArray()) return;
	for (rj::SizeType k = 0; k < lsv.Size(); ++k){
		LightSourceImpl *&cimpl = pimpl->impls.AppendNew();
		cimpl = 0;
		rj::Value &lsiv = lsv[k];
		if (!lsiv.IsObject()) continue;
		const char *type = rj_getstring(lsiv["type"], "");
		if (std::strcmp(type, "parallel") == 0){
			cimpl = new ParallelRayImpl(lsiv);
		}else if (std::strcmp(type, "osram") == 0){
			cimpl = new OsramRayImpl(lsiv);
		}
		if (cimpl && !cimpl->valid){
			delete cimpl; cimpl = 0;
		}
		if (cimpl) pimpl->count += cimpl->Count();
		pimpl->accum_count.Append(*pimpl->accum_count.Last() + pimpl->count);
	}
}
LightSources::~LightSources(){
	delete pimpl; pimpl = 0;
}
int LightSources::LSCount() const {
	return pimpl->impls.Count();
}
int64_t LightSources::RayCount() const{
	return pimpl->count;
}
int64_t LightSources::RayCount(int lsidx) const{
	if (lsidx < 0 || lsidx >= pimpl->impls.Count()) return 0;
	if (!pimpl->impls[lsidx]) return 0;
	return pimpl->impls[lsidx]->Count();
}

bool LightSources::Get(int64_t idx, int &blockseed, int &lsidx, int64_t &rayidx, ON_3dRay &ray, double &flux, double &wavelength) const{
	int64_t ls_start = pimpl->accum_count[pimpl->cur_ls];
	LightSourceImpl *cur_impl = 0;
	if (ls_start <= idx && idx < pimpl->accum_count[pimpl->cur_ls+1]){
		cur_impl = pimpl->impls[pimpl->cur_ls];
		lsidx = pimpl->cur_ls;
	}else{
		lsidx = static_cast<int>(
			std::upper_bound(pimpl->accum_count.First(), pimpl->accum_count.Last()+1, idx)
			- pimpl->accum_count.First()) - 1;
		if (lsidx < 0 || lsidx >= pimpl->impls.Count()){
			lsidx = -1;
			return false;
		}
		pimpl->cur_ls = lsidx;
		ls_start = pimpl->accum_count[lsidx];
	}
	rayidx = idx - ls_start;
	return pimpl->impls[lsidx] ? pimpl->impls[lsidx]->Get(idx - ls_start, blockseed, ray, flux, wavelength) : false;
}

/// ========== Materials ==========
// 材質データ
struct Materials::Impl{
	struct Material{
		ON_String name;
		double refl_rate, transmittance, refr_index, absorpt_coef;
		Material(){
			name = "";
			refl_rate = 0, transmittance = 0, refr_index = 0;
			absorpt_coef = 0;
		}
	};
	ON_ClassArray<Material> mats;
};

Materials::Materials(rj::Value &mtv, rj::Value &shv, ON_SimpleArray<int> &shape2matidx){
	std::map<ON_String, rj::SizeType> matname2matidx;
	pimpl = new Impl();
	if (!mtv.IsArray()) return;
	for (rj::SizeType k = 0; k < mtv.Size(); ++k){
		Impl::Material &mat = pimpl->mats.AppendNew();
		rj::Value &matv = mtv[k];
		if (!matv.IsObject()) continue;
		mat.name          = rj_getstring(matv["name"]);
		mat.refl_rate     = rj_getdouble(matv["refl_rate"]);
		mat.transmittance = rj_getdouble(matv["transmittance"]);
		mat.refr_index    = rj_getdouble(rj_getmember(matv, "refr_index"));
		mat.absorpt_coef  = rj_getdouble(rj_getmember(matv, "absorpt_coef"));
		matname2matidx.insert(std::make_pair(mat.name, k));
	}

	if (!shv.IsArray()) return;

	shape2matidx.Destroy();
	for (rj::SizeType k = 0; k < shv.Size(); ++k){
		ON_String matname = rj_getstring(shv[k]["material"]);
		std::map<ON_String, rj::SizeType>::iterator iter = matname2matidx.find(matname);
		if (iter == matname2matidx.end()) shape2matidx.Append(-1);
		else shape2matidx.Append(iter->second);
	}
}

Materials::~Materials(){
	delete pimpl;
}

int Materials::Count() const{
	return pimpl->mats.Count();
}

bool Materials::GetParams(int midx, double wavelength, double &refl_rate, double &refr_index, double &transmittance, double &absorpt_coef) const{
	if (midx < 0 || midx >= Count()) return false;
	// Todo: 波長依存性の対応
	Impl::Material &mat = pimpl->mats[midx];
	refl_rate = mat.refl_rate;
	refr_index = mat.refr_index;
	transmittance = mat.transmittance;
	absorpt_coef = mat.absorpt_coef;
	return true;
}
bool Materials::CalcBSDF(const ON_3dVector &emitdir, bool transmitted, int &seed, ON_3dVector &scattered_dir) const{
	// Todo: 実装
	ON_3dVector dir = emitdir;
	scattered_dir = dir;
	return true;
}
