package lpv.simulation.propagator;

import javax.vecmath.Point3i;
import javax.vecmath.Vector2f;
import javax.vecmath.Vector3f;
import javax.vecmath.Vector4f;

import lpv.LPV;
import lpv.SH;

public class ExperimentalPropagator implements Propagator {
	// Experimental custom propagator based on the nVidia propagator. 
	
	/*private static final Point3i offsets[] = new Point3i[]{
			new Point3i(1,0,0), new Point3i(-1,0,0),
			new Point3i(0,1,0), new Point3i(0,-1,0),
			new Point3i(0,0,1), new Point3i(0,0,-1)
		};
	private float factor;
	
	public ExperimentalPropagator(float factor) {
		this.factor = factor;
	}
	
	@Override
	public void propagate(LPV in, LPV out, LPV acc, Point3i pos, int iteration) {
		Vector4f o = out.get(pos.x, pos.y, pos.z);
		o.set(0.0f, 0.0f, 0.0f, 0.0f);
		
		for(int neighbor = 0; neighbor < 6; neighbor++) {
			Point3i offset = offsets[neighbor];
			
			Point3i targetPos = new Point3i(pos);
			targetPos.add(offset);
			
			if(targetPos.x >= 0 && targetPos.x < in.getWidth() && 
					targetPos.y >= 0 && targetPos.y < in.getHeight() && 
					targetPos.z >= 0 && targetPos.z < in.getDepth()) {
				
				Vector4f i = in.get(targetPos.x, targetPos.y, targetPos.z);
				
				
				for(int face = 0; face < 6; face++) {
					
					Vector3f dir = new Vector3f(offsets[face].x, offsets[face].y, offsets[face].z);
					dir.scale(0.5f);
					dir.sub(new Vector3f(offset.x, offset.y, offset.z));
					
					float len = dir.length();
					dir.scale(1.0f / len);
					
					float solidAngle;
					if(len <= 0.5f) {
						solidAngle = 0.0f; 
					} else {
						solidAngle = (len >= 1.5f) ? 22.95668f/(4.0f*180.0f) : 24.26083f/(4.0f*180.0f);
					}
					
					Vector4f dirSH = SH.construct(dir.x, dir.y, dir.z);
					
					float flux = factor * solidAngle * Math.max(0.0f, i.dot(dirSH));
					
					Vector4f coeffs = SH.clampedCosineLobe(offsets[face].x, offsets[face].y, offsets[face].z);
					coeffs.scale(flux);
					o.add(coeffs);
				}
			}
		}
		
		if(iteration == 0) {
			o.add(in.get(pos.x, pos.y, pos.z));
		}
		
		acc.get(pos.x, pos.y, pos.z).add(o);
	}*/
	
	
	private float factor;

	public ExperimentalPropagator(float factor) {
		this.factor = factor;
	}

	@Override
	public void propagate(LPV in, LPV out, LPV acc, LPV geometry, Point3i pos, int iteration) {
		Vector4f pixelCoeffs = new Vector4f(0.0f, 0.0f, 0.0f, 0.0f);

		// 6-point axial gathering stencil "cross"
		for(int i = -1; i <= 1; i++) {
			for(int j = -1; j <= 1; j++) {
				for(int k = -1; k <= 1; k++) {
					if(i != 0 || j != 0 || k != 0) {
						IVPropagateDir(pixelCoeffs, in, pos, new Point3i(i, j, k));
					}
				}
			}
		}
		
		pixelCoeffs.scale(factor);
		if(iteration == 0) {
			//pixelCoeffs.add(in.get(pos.x, pos.y, pos.z));
			acc.get(pos.x, pos.y, pos.z).add(in.get(pos.x, pos.y, pos.z));
		}
		out.get(pos.x, pos.y, pos.z).set(pixelCoeffs);
		acc.get(pos.x, pos.y, pos.z).add(pixelCoeffs);
	}

	void IVPropagateDir(Vector4f pixelCoeffs, LPV in, Point3i pos,
			Point3i nOffset) {
		// get adjacent cell's SH coeffs
		Point3i p = new Point3i(pos.x + nOffset.x, pos.y + nOffset.y,
				pos.z + nOffset.z);
		
		Vector4f sampleCoeffs = in.get(p.x, p.y, p.z);
		
		if(sampleCoeffs.x != 0.0f || sampleCoeffs.y != 0.0f || sampleCoeffs.z != 0.0f || sampleCoeffs.w != 0.0f) {
		
		
			// generate function for incoming direction from adjacent cell
			Vector3f normalizedOffset = new Vector3f(nOffset.x, nOffset.y, nOffset.z);
			normalizedOffset.normalize();
			Vector4f shIncomingDirFunction = Cone90Degree(normalizedOffset);
			// integrate incoming radiance with this function
			float incidentLuminance = Math.max(0.0f,
					sampleCoeffs.dot(shIncomingDirFunction));
			
			Vector4f shOutgoing = Cone90Degree(getPrimaryDirection(sampleCoeffs));
			shOutgoing.scale(incidentLuminance);
			
			pixelCoeffs.add(shOutgoing);
		}
	}

	private Vector4f Cone90Degree(Vector3f vcDir) {
		return SHProjectCone(vcDir, (float) Math.toRadians(90.0f));
	}
	
	private Vector4f Cone45Degree(Vector3f vcDir) {
		return SHProjectCone(vcDir, (float) Math.toRadians(45.0f));
	}
	
	private Vector4f Cone180Degree(Vector3f vcDir) {
		return SHProjectCone(vcDir, (float) Math.toRadians(180.0f));
	}
	
	private Vector3f getPrimaryDirection(Vector4f sh) {
		Vector3f dir = new Vector3f(sh.w, sh.y, -sh.z);
		if(dir.x == 0.0f && dir.y == 0.0f && dir.z == 0.0f) {
			return dir;
		} else {
			dir.normalize();
			return dir;
		}
	}

	private Vector4f SHProjectCone(Vector3f vcDir, float angle) {
		Vector2f vZHCoeffs = new Vector2f(
				0.5f * (1.0f - (float) Math.cos(angle)), // 1/2 (1 -
															// Cos[\[Alpha]])
				0.75f * (float) Math.sin(angle) * (float) Math.sin(angle)); // 3/4
																			// Sin[\[Alpha]]^2
		return SHRotate(vcDir, vZHCoeffs);
	}

	private Vector4f SHRotate(Vector3f vcDir, Vector2f vZHCoeffs) {
		// compute sine and cosine of thetta angle
		// beware of singularity when both x and y are 0 (no need to rotate at
		// all)

		Vector2f theta12_cs;
		if(vcDir.x == 0.0f && vcDir.y == 0.0f) {
			theta12_cs = new Vector2f(0.0f, 0.0f);
		} else {
			theta12_cs = new Vector2f(vcDir.x, vcDir.y);
			theta12_cs.normalize();
		}

		// compute sine and cosine of phi angle
		Vector2f phi12_cs = new Vector2f();
		phi12_cs.x = (float) Math.sqrt(1.0f - vcDir.z * vcDir.z);
		phi12_cs.y = vcDir.z;

		Vector4f vResult = new Vector4f();
		// The first band is rotation-independent
		vResult.x = vZHCoeffs.x;
		// rotating the second band of SH
		vResult.y = vZHCoeffs.y * phi12_cs.x * theta12_cs.y;
		vResult.z = -vZHCoeffs.y * phi12_cs.y;
		vResult.w = vZHCoeffs.y * phi12_cs.x * theta12_cs.x;
		return vResult;
	}
}
