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

#define NOMINMAX

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include "opennurbs.h"
#include <hash_map>

// SphereGrid
struct SphereGrid{
	int divides; // divides は 1024 まで
	int type;
	// gidx : Geometrical Index: ix + iy * divides + iz * divides * divides
	// nidx : 点列の順番
	stdext::hash_map<int, int> gidx2nidx;
	ON_SimpleArray<int> nidx2gidx;
	ON_ClassArray<ON_3dPoint> vertices;
	ON_3dPoint center;
	double radius;

	void DecomposeGeomIndex(int gidx, int &ix, int &iy, int &iz){
		div_t dt = std::div(gidx, (divides + 1) * (divides + 1));
		iz = dt.quot;
		dt = std::div(dt.rem, (divides + 1));
		iy = dt.quot, ix = dt.rem;
	}
	int ComposeGeomIndex(int ix, int iy, int iz){
		return ix + iy * (divides + 1) + iz * (divides + 1) * (divides + 1);
	}

		//         4          5
		//          *----------*                ┐Z
		//        ／.        ／|　　　　      ／
		//      ／  .      ／  |            ＊---> X
		//  0 ／    .  1 ／    |            │
		//   *----------*      |            │
		//   |    6 . . | . .7 *            ↓Y
		//   |    .     |    ／
		//   |  .       |  ／
		//  2|         3|／
		//   *----------*
		//

	// divides_exponent: 1 => 2分割, 2 => 4分割, 3 => 8分割, ...
	void Create(const ON_3dPoint &center, double radius, int divides_exponent){
		this->center = center;
		this->radius = radius;
		gidx2nidx.clear();
		nidx2gidx.Empty();
		vertices.Empty();
		divides = divides_exponent <= 0 ? 1 : 2;
		for (int i = 1; i < divides_exponent; ++i) divides *= 2;

		ON_ClassArray<ON_3dVector> directions;
		// まずは頂点から
		for (int i = 0; i < 8; ++i){
			ON_3dVector v((i & 1) ? 1 : -1, (i & 2) ? 1 : -1, (i & 4) ? 1 : -1);
			v.Unitize();
			directions.Append(v);
			int ix = (i & 1) ? divides : 0, iy = (i & 2) ? divides : 0, iz = (i & 4) ? divides : 0;
			int gidx = ComposeGeomIndex(ix, iy, iz);
			gidx2nidx.insert(std::make_pair(gidx, nidx2gidx.Count()));
			nidx2gidx.Append(gidx);
		}
		ON_ClassArray<std::pair<int, int> > nindices;
		for (int k = 0; k < 2; ++k){
			if (k == 0){
				// 次にエッジ
				// Z軸に垂直な4本
				nindices.Append(std::make_pair(3, 7));
				nindices.Append(std::make_pair(2, 6));
				nindices.Append(std::make_pair(1, 5));
				nindices.Append(std::make_pair(0, 4));
				// Y軸に垂直な4本
				nindices.Append(std::make_pair(5, 7));
				nindices.Append(std::make_pair(4, 6));
				nindices.Append(std::make_pair(1, 3));
				nindices.Append(std::make_pair(0, 2));
				// X軸に垂直な4本
				nindices.Append(std::make_pair(6, 7));
				nindices.Append(std::make_pair(4, 5));
				nindices.Append(std::make_pair(2, 3));
				nindices.Append(std::make_pair(0, 1));
			}else if(k == 1){
				// 2本の向かい合うエッジ間を補間
				// nidxが
				//   8                     〜 8 +     (divides - 2) => 頂点 0 - 1 間のエッジ点列
				//   8 +     (divides - 1) 〜 8 + 2 * (divides - 2) => 頂点 2 - 3 間のエッジ点列
				//   8 + 2 * (divides - 1) 〜 8 + 3 * (divides - 2) => 頂点 4 - 5 間のエッジ点列
				//   8 + 3 * (divides - 1) 〜 8 + 4 * (divides - 2) => 頂点 6 - 7 間のエッジ点列
				//   8 + 4 * (divides - 1) 〜 8 + 5 * (divides - 2) => 頂点 0 - 2 間のエッジ点列
				//   8 + 5 * (divides - 1) 〜 8 + 6 * (divides - 2) => 頂点 1 - 3 間のエッジ点列
				//   8 + 6 * (divides - 1) 〜 8 + 7 * (divides - 2) => 頂点 4 - 6 間のエッジ点列
				//   8 + 7 * (divides - 1) 〜 8 + 8 * (divides - 2) => 頂点 5 - 7 間のエッジ点列

				// (0 - 1) => (2 - 3), (2 - 3) => (6 - 7), (6 - 7) => (4 - 5), (4 - 5) => (0 - 1),
				// (1 - 3) => (5 - 7), (4 - 6) => (0 - 2) の順に、2本のエッジ間を補間
				int divi_1 = divides - 1;
				int nidx_edgetop[6][2] = {
					{8,              8 +     divi_1}, // (0 - 1) => (2 - 3)
					{8 +     divi_1, 8 + 3 * divi_1}, // (2 - 3) => (6 - 7)
					{8 + 3 * divi_1, 8 + 2 * divi_1}, // (2 - 3) => (4 - 5)
					{8 + 2 * divi_1, 8             }, // (4 - 5) => (0 - 1)
					{8 + 5 * divi_1, 8 + 7 * divi_1}, // (1 - 3) => (5 - 7)
					{8 + 6 * divi_1, 8 + 4 * divi_1}  // (4 - 6) => (0 - 3)
				};
				for (int j = 5; j >= 0; --j){
					for (int i = 0; i < divi_1; ++i){
						nindices.Append(std::make_pair(nidx_edgetop[j][0] + i, nidx_edgetop[j][1] + i));
					}
				}
			}
			while(nindices.Count()){
				std::pair<int, int> pr = *nindices.Last(); nindices.SetCount(nindices.Count()-1);
				int ix1, iy1, iz1, ix2, iy2, iz2;
				int ixc, iyc, izc;
				DecomposeGeomIndex(nidx2gidx[pr.first], ix1, iy1, iz1);
				DecomposeGeomIndex(nidx2gidx[pr.second], ix2, iy2, iz2);
				if (ix1 != ix2 && std::abs(ix2 - ix1) > 1) ixc = (ix1 + ix2) / 2, iyc = iy1, izc = iz1;
				else if (iy1 != iy2 && std::abs(iy2 - iy1) > 1) ixc = ix1, iyc = (iy1 + iy2) / 2, izc = iz1;
				else if (iz1 != iz2 && std::abs(iz2 - iz1) > 1) ixc = ix1, iyc = iy1, izc = (iz1 + iz2) / 2;
				else continue;
				ON_3dVector vc = directions[pr.first] + directions[pr.second];
				vc.Unitize();
				directions.Append(vc);

				int gidxc = ComposeGeomIndex(ixc, iyc, izc);
				int nidxc = nidx2gidx.Count();
				gidx2nidx[gidxc] = nidxc;
				gidx2nidx.insert(std::make_pair(gidxc, nidxc));
				nidx2gidx.Append(gidxc);
				nindices.Append(std::make_pair(pr.first, nidxc));
				nindices.Append(std::make_pair(nidxc, pr.second));
			}
		}
		vertices.SetCapacity(directions.Count()), vertices.SetCount(vertices.Capacity());
		for (int i = 0; i < directions.Count(); ++i){
			vertices[i] = center + directions[i] * radius;
		}
	}

	void GetNeighborPoints(int nidx, ON_SimpleArray<int> &nb){
		int gidx = nidx2gidx[nidx];
		int divi1 = divides + 1;
		int divi2 = divi1 * divi1;

		int ix, iy, iz;
		DecomposeGeomIndex(gidx, ix, iy, iz);
		bool yz, zx, xy;
		int ptype = GetPointType(ix, iy, iz, yz, zx, xy);
		if (ptype == 0){
			if (!yz) nb.Append(gidx2nidx[gidx - 1]), nb.Append(gidx2nidx[gidx + 1]);
			if (!zx) nb.Append(gidx2nidx[gidx - divi1]), nb.Append(gidx2nidx[gidx + divi1]);
			if (!xy) nb.Append(gidx2nidx[gidx - divi2]), nb.Append(gidx2nidx[gidx + divi2]);
		}else if (ptype == 1){
			if (!yz){ // X軸に平行なエッジ
				nb.Append(gidx2nidx[gidx - 1]), nb.Append(gidx2nidx[gidx + 1]);
				nb.Append(gidx2nidx[iy == 0 ? gidx + divi1: gidx - divi1]);
				nb.Append(gidx2nidx[iz == 0 ? gidx + divi2 : gidx - divi2]);
			}else if (!zx){ // Y軸に平行なエッジ
				nb.Append(gidx2nidx[ix == 0 ? gidx + 1 : gidx - 1]);
				nb.Append(gidx2nidx[gidx - divi1]), nb.Append(gidx2nidx[gidx + divi1]);
				nb.Append(gidx2nidx[iz == 0 ? gidx + divi2 : gidx - divi2]);
			}else if (!xy){ // Z軸に平行なエッジ
				nb.Append(gidx2nidx[ix == 0 ? gidx + 1 : gidx - 1]);
				nb.Append(gidx2nidx[iy == 0 ? gidx + divi1 : gidx - divi1]);
				nb.Append(gidx2nidx[gidx - divi2]), nb.Append(gidx2nidx[gidx + divi2]);
			}
		}else if (ptype == 2){
			nb.Append(gidx2nidx[ix == 0 ? gidx + 1 : gidx - 1]);
			nb.Append(gidx2nidx[iy == 0 ? gidx + divi1 : gidx - divi1]);
			nb.Append(gidx2nidx[iz == 0 ? gidx + divi2 : gidx - divi2]);
		}
	}
	// type 0:面の中身
	// type 1:縁
	// type 2:頂点
	int GetPointType(int ix, int iy, int iz, bool &yz, bool &zx, bool &xy){
		int cnt = 0;
		yz = zx = xy = false;
		if (ix == 0 || ix == divides) ++cnt, yz = true;
		if (iy == 0 || iy == divides) ++cnt, zx = true;
		if (iz == 0 || iz == divides) ++cnt, xy = true;
		return cnt - 1;
	}
	int Count(){
		return vertices.Count();
	}
	const ON_3dPoint Point(int idx){
		return vertices[idx];
	}
	void UpdatePoint(int idx, const ON_3dPoint &pt_new){
		vertices[idx] = pt_new;
	}
	void Save(){
		ONX_Model m;
		for (int k = 0; k < 6; ++k){
			ON_String xyz; xyz.Format("xyz_%d.txt", k+1);
			FILE *fp = std::fopen(xyz, "w");
			ON_PointCloud pc;
			for (int j = 0; j <= divides; ++j){
				for (int i = 0; i <= divides; ++i){
					int gidx;
					switch(k){
						case 0: gidx = ComposeGeomIndex(i, j, 0); break;
						case 1: gidx = ComposeGeomIndex(i, j, divides); break;
						case 2: gidx = ComposeGeomIndex(i, 0, j); break;
						case 3: gidx = ComposeGeomIndex(i, divides, j); break;
						case 4: gidx = ComposeGeomIndex(0, i, j); break;
						case 5: gidx = ComposeGeomIndex(divides, i, j); break;
					}
					int nidx = gidx2nidx[gidx];
					ON_3dPoint &pt = vertices[nidx];
					pc.AppendPoint(pt);
					int r = ((k+1) & 1) ? 255 : 0, g = ((k+1) & 2) ? 255 : 0, b = ((k+1) & 4) ? 255 : 0;
					std::fprintf(fp, "%f %f %f %d %d %d\n", pt.x, pt.y, pt.z, r, g, b);
				}
			}
			ONX_Model_Object &obj = m.m_object_table.AppendNew();
			obj.m_object = pc.Duplicate(), obj.m_bDeleteObject = true;
			std::fclose(fp);
		}
		m.Write("spheregrid.3dm", 4);
	}
};

int main(){
	SphereGrid sg;
	sg.Create(ON_3dPoint(10, 0, 5), 20, 6);
	sg.Save();
	return 0;
}
