package lpv.simulation.propagator;

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

import lpv.LPV;

public class BenjaminThautPropagator implements Propagator {
	// Original Author: Benjamin Thaut
	// /home/dt07jo1/Documents/lpv-release/code/shader/lpv_propagation.xml
	// http://3d.benjamin-thaut.de/?p=16
	
	private static final Vector4f SH_c = new Vector4f(0.282094792f,-0.488602512f,0.488602512f,0.488602512f);
	private static final Vector4f SH_cosLobe_c = new Vector4f(0.25f,-0.5f,0.5f,0.5f);
	
	private static final float directFaceSubtendedSolidAngle = 0.12753712f;
	private static final float sideFaceSubtendedSolidAngle = 0.13478556f;
	
	private static final float side1 = 0.447213595f;//1 / sqrt(5)
	private static final float side2 = 0.894427190f;//2 / sqrt(5)
	
	private float factor;
	
	public BenjaminThautPropagator(float factor) {
		this.factor = factor;
	}
	
	@Override
	public void propagate(LPV in, LPV out, LPV acc, LPV geometry, Point3i pos, int iteration) {
		Vector4f o = out.get(pos.x, pos.y, pos.z);
		o.set(0.0f, 0.0f, 0.0f, 0.0f);
		
		Matrix3f[] dirs = new Matrix3f[] {
				new Matrix3f( // Z+
					1.0f, 0.0f, 0.0f,
					0.0f, 1.0f, 0.0f,
					0.0f, 0.0f, 1.0f
				),
				new Matrix3f( // Z-
					-1.0f, 0.0f, 0.0f,
					 0.0f, 1.0f, 0.0f,
					 0.0f, 0.0f,-1.0f
				),
				new Matrix3f( // X+
					 0.0f, 0.0f, 1.0f,
					 0.0f, 1.0f, 0.0f,
					-1.0f, 0.0f, 0.0f
				),
				new Matrix3f( // X-
					0.0f, 0.0f,-1.0f,
					0.0f, 1.0f, 0.0f,
					1.0f, 0.0f, 0.0f
				),
				new Matrix3f( // Y+
					1.0f, 0.0f, 0.0f,
					0.0f, 0.0f, 1.0f,
					0.0f,-1.0f, 0.0f
				),
				new Matrix3f( // Y-
					1.0f, 0.0f, 0.0f,
					0.0f, 0.0f,-1.0f,
					0.0f, 1.0f, 0.0f
				),
		};
		
		for(Matrix3f dir : dirs) {
			propagate(in, pos, dir, o);
		}
		
		acc.get(pos.x, pos.y, pos.z).add(o);
	}
	
	private Vector4f SH_evaluate(Vector3f dir){
		return new Vector4f(SH_c.x, SH_c.y * dir.y, SH_c.z * dir.z, SH_c.w * dir.x);
	}
	
	private Vector4f SH_evaluateCosineLobe_direct(Vector3f dir) {
		return new Vector4f(SH_cosLobe_c.x, SH_cosLobe_c.y * dir.y, SH_cosLobe_c.z * dir.z, SH_cosLobe_c.w * dir.x);
	}

	private void propagate(LPV input, Point3i pos, Matrix3f orientation, Vector4f output) {
		
		Vector3f mainDirection = new Vector3f(0.0f, 0.0f, 1.0f);
		orientation.transform(mainDirection);
		
		Vector3f dirRangeMul = new Vector3f(mainDirection);
		dirRangeMul.scale(factor);
		Vector4f mainDirectionSH = SH_evaluate(dirRangeMul);
		
		Vector4f mainDirectionHemi = SH_evaluateCosineLobe_direct(mainDirection);
		
		Point3i loadPos = new Point3i(
				pos.x - (int)mainDirection.x, 
				pos.y - (int)mainDirection.y, 
				pos.z - (int)mainDirection.z);
		
		if(loadPos.x >= 0 && loadPos.x < input.getWidth() && 
				loadPos.y >= 0 && loadPos.y < input.getHeight() &&
				loadPos.z >= 0 && loadPos.z < input.getDepth()) {
			
			Vector4f neighbour = input.get(loadPos.x, loadPos.y, loadPos.z);
			
			// Evaluate main face
			float flux = Math.max(0.0f, mainDirectionSH.dot(neighbour));
			Vector4f result = new Vector4f(mainDirectionHemi);
			result.scale(flux * directFaceSubtendedSolidAngle);
			
			
			
			
			// Evaluate side faces
			Vector3f sideDirection = new Vector3f(side1, 0.0f, side2);
			orientation.transform(sideDirection);
			
			Vector3f reproDirection = new Vector3f(1.0f, 0.0f, 0.0f);
			orientation.transform(reproDirection);
			
			Vector4f sideDirectionSH = SH_evaluate(sideDirection);
			sideDirectionSH.scale(factor);
			Vector4f sideDirectionHemi = SH_evaluateCosineLobe_direct(reproDirection);
			
			flux = Math.max(0.0f, sideDirectionSH.dot(neighbour));
			
			Vector4f contribution = new Vector4f(sideDirectionHemi);
			contribution.scale(flux * sideFaceSubtendedSolidAngle);
			result.add(contribution);
			
			
			
			
			sideDirection = new Vector3f(-side1, 0.0f, side2);
			orientation.transform(sideDirection);
			
			reproDirection = new Vector3f(-1.0f, 0.0f, 0.0f);
			orientation.transform(reproDirection);
			
			sideDirectionSH = SH_evaluate(sideDirection);
			sideDirectionSH.scale(factor);
			sideDirectionHemi = SH_evaluateCosineLobe_direct(reproDirection);
			
			flux = Math.max(0.0f, sideDirectionSH.dot(neighbour));
			
			contribution = new Vector4f(sideDirectionHemi);
			contribution.scale(flux * sideFaceSubtendedSolidAngle);
			result.add(contribution);
			
			
			
			
			sideDirection = new Vector3f(0.0f, side1, side2);
			orientation.transform(sideDirection);
			
			reproDirection = new Vector3f(0.0f, 1.0f, 0.0f);
			orientation.transform(reproDirection);
			
			sideDirectionSH = SH_evaluate(sideDirection);
			sideDirectionSH.scale(factor);
			sideDirectionHemi = SH_evaluateCosineLobe_direct(reproDirection);
			
			flux = Math.max(0.0f, sideDirectionSH.dot(neighbour));
			
			contribution = new Vector4f(sideDirectionHemi);
			contribution.scale(flux * sideFaceSubtendedSolidAngle);
			result.add(contribution);
			
			
			
			
			sideDirection = new Vector3f(0.0f, -side1, side2);
			orientation.transform(sideDirection);
			
			reproDirection = new Vector3f(0.0f, -1.0f, 0.0f);
			orientation.transform(reproDirection);
			
			sideDirectionSH = SH_evaluate(sideDirection);
			sideDirectionSH.scale(factor);
			sideDirectionHemi = SH_evaluateCosineLobe_direct(reproDirection);
			
			flux = Math.max(0.0f, sideDirectionSH.dot(neighbour));
			
			contribution = new Vector4f(sideDirectionHemi);
			contribution.scale(flux * sideFaceSubtendedSolidAngle);
			result.add(contribution);
			
			
			
			
			output.add(result);
		}
	}
}
