content: /*
   * Copyright (C) 2015 Benedikt Lorch
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  package edu.stanford.rsl.conrad.rendering;
  
  import java.util.ArrayList;
  import java.util.HashMap;
  import java.util.LinkedList;
  import java.util.List;
  import java.util.Map;
  import java.util.Queue;
  import java.util.Stack;
  
  import edu.stanford.rsl.conrad.geometry.AbstractCurve;
  import edu.stanford.rsl.conrad.geometry.AbstractShape;
  import edu.stanford.rsl.conrad.geometry.shapes.compound.CompoundShape;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.Edge;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.ProjectPointToLineComparator;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine;
  import edu.stanford.rsl.conrad.geometry.shapes.simple.Triangle;
  import edu.stanford.rsl.conrad.numerics.SimpleOperators;
  import edu.stanford.rsl.conrad.numerics.SimpleVector;
  import edu.stanford.rsl.conrad.physics.PhysicalObject;
  import edu.stanford.rsl.conrad.physics.PhysicalPoint;
  import edu.stanford.rsl.conrad.physics.detector.MaterialPathLengthDetector;
  import edu.stanford.rsl.conrad.physics.materials.Material;
  import edu.stanford.rsl.conrad.utils.CONRAD;
  import edu.stanford.rsl.conrad.utils.Configuration;
  
  public class WatertightRayTracer extends AbstractRayTracer {
  
  	
  	protected int debug = 1;
  	
  	public WatertightRayTracer() {
  		this.comparator = new ProjectPointToLineComparator();
  		if (debug > 0) {
  			System.out.println("This is the WATERTIGHT RAY TRACER");
  		}
  	}
  
  	
  	// Copied from SimpleRayTracer
  	@Override
  	protected ArrayList<PhysicalObject> computeMaterialIntersectionSegments(PhysicalPoint [] rayList){
  		if (Configuration.getGlobalConfiguration().getDetector() instanceof MaterialPathLengthDetector) {
  			return computeMaterialIntersectionSegmentsMatPathLen(rayList);
  		}
  		
  		ArrayList<PhysicalObject> segments = new ArrayList<PhysicalObject>();
  
  		Stack<PhysicalObject> materialStack = new Stack<PhysicalObject>();
  		materialStack.push(rayList[0].getObject());
  
  		// Iterate over hits
  		for (int k=1; k < rayList.length; k++){
  			PhysicalObject obj = new PhysicalObject();
  			Edge edge = new Edge(rayList[k-1], rayList[k]);
  			obj.setShape(edge);
  			if(materialStack.isEmpty()) {
  				obj.setMaterial(scene.getBackgroundMaterial());
  				obj.setNameString("Background");
  			} else {
  				// Take top of stack, but do not remove the top element
  				PhysicalObject current = materialStack.peek();
  				obj.setMaterial(current.getMaterial());
  				obj.setNameString(current.getNameString());
  			}
  			segments.add(obj);
  			PhysicalObject nextObject = rayList[k].getObject();
  
  			if (materialStack.contains(nextObject)){
  				materialStack.remove(nextObject);
  			} else {
  				materialStack.push(nextObject);
  			}
  		}
  
  		return segments;
  	}
  	
  	
  	/**
  	 * Compute material segments for material path length detector
  	 */
  	protected ArrayList<PhysicalObject> computeMaterialIntersectionSegmentsMatPathLen(PhysicalPoint [] rayList){
  		ArrayList<PhysicalObject> segments = new ArrayList<PhysicalObject>();
  		
  		// Assuming that the list of hits is sorted in ray direction
  		// For the Material Path Length Detector, we need to form segments for each material
  		
  		// Sort by material
  		HashMap<Material, List<PhysicalPoint>> hitMap = new HashMap<>();
  		for (int i=0; i<rayList.length; i++) {
  			PhysicalPoint p = rayList[i];
  			if (!hitMap.containsKey(p.getObject().getMaterial())) {
  				hitMap.put(p.getObject().getMaterial(), new LinkedList<PhysicalPoint>());
  			}
  			hitMap.get(p.getObject().getMaterial()).add(p);
  		}
  
  		for (Map.Entry<Material, List<PhysicalPoint>> hitsByMaterial : hitMap.entrySet()) {
  			List<PhysicalPoint> hits = hitsByMaterial.getValue();
  
  			// Iterate over hits
  			for (int k=0; k < hits.size(); k = k + 2){
  				PhysicalObject obj = new PhysicalObject();
  				Edge edge = new Edge(hits.get(k), hits.get(k+1));
  				obj.setShape(edge);
  				obj.setMaterial(hitsByMaterial.getKey());
  				segments.add(obj);
  			}
  		}
  
  		return segments;
  	}
  	
  	
  	/**
  	 * The algorithm assumes that direction's largest absolute value is stored in the z-component
  	 * @param vec direction of ray
  	 * @return array of size 3 containing the permutation of the direction's components such that the third dimension holds the largest absolute value
  	 */
  	private int[] getMaxDirection(SimpleVector vec) {
  		assert(vec.getLen() == 3);
  		int[] idx = new int[3];
  		
  		if (Math.abs(vec.getElement(0)) > Math.abs(vec.getElement(1))) {
  			if (Math.abs(vec.getElement(0)) > Math.abs(vec.getElement(2))) {
  				// x > y and x > z
  				idx[2] = 0;
  				idx[0] = 1;
  				idx[1] = 2;
  			} else {
  				// z > x > y
  				idx[2] = 2;
  				idx[0] = 0;
  				idx[1] = 1;
  			}
  		} else {
  			if (Math.abs(vec.getElement(1)) > Math.abs(vec.getElement(2))) {
  				// y > x and y > z
  				idx[2] = 1;
  				idx[0] = 2;
  				idx[1] = 0;
  			} else {
  				// z > y > x
  				idx[2] = 2;
  				idx[0] = 0;
  				idx[1] = 1;
  			}
  		}
  		
  		// Swap kx and ky dimension if max direction is negative in order to preserve winding direction of triangles
  		if (vec.getElement(idx[2]) < 0.d) {
  			int temp = idx[0];
  			idx[0] = idx[1];
  			idx[1] = temp;
  		}
  		
  		return idx;
  	}
  	
  
  
  	/**
  	 * Implementation of watertight algorithm as proposed by Sven Woop, Carsten Benthin and Ingo Wald
  	 * See http://jcgt.org/published/0002/01/05/paper.pdf and https://github.com/embree/embree/blob/v1.1_watertight/rtcore/triangle/triangle1i_intersector1_watertight.h
  	 */
  	@Override
  	public ArrayList<PhysicalPoint> intersectWithScene(AbstractCurve ray) {
  		
  		StraightLine line = (StraightLine) ray;
  		
  		double[] dir = line.getDirection().copyAsDoubleArray();
  		PointND origin = line.getPoint();
  		// Calculate dimensions where the ray direction is maximal
  		int[] kPerm = getMaxDirection(line.getDirection());
  		
  		// For simplicity, move permutation indices into primitive variables
  		int kx = kPerm[0];
  		int ky = kPerm[1];
  		int kz = kPerm[2];
  		
  		// Calculate shear constants
  		double Sz = 1.0d / dir[kz];
  		double Sx = dir[kx] * Sz;
  		double Sy = dir[ky] * Sz;
  		
  		ArrayList<PhysicalPoint> hits = new ArrayList<>();
  		
  		Queue<AbstractShape> queue;
  		// Iterate over objects of scene
  		for (PhysicalObject obj : scene) {
  			// This algorithm can only process triangles, which are usually embedded in a compound shape
  			queue = new LinkedList<>();
  			queue.add(obj.getShape());
  			while (!queue.isEmpty()) {
  				AbstractShape shape = queue.poll();
  				if (shape instanceof CompoundShape) {
  					CompoundShape cs = (CompoundShape) shape;
  
  					// Test for hits on bounding box. In case the ray does not even hit the bounding box, we will not find an intersection with a triangle
  					if (cs.getHitsOnBoundingBox(ray).size() > 0) {
  						queue.addAll(cs);
  					}
  				}
  				else if (!(shape instanceof Triangle)) {
  					System.err.println("Only triangles can be intersected by the watertight algorithm. As the current shape is not a triangle, it will be skipped.");
  				}
  				else {
  					Triangle triangle = (Triangle) shape;
  					
  					// Calculate vertices relative to ray origin
  					PointND A = new PointND(SimpleOperators.subtract(triangle.getA().getAbstractVector(), origin.getAbstractVector()));
  					PointND B = new PointND(SimpleOperators.subtract(triangle.getB().getAbstractVector(), origin.getAbstractVector()));
  					PointND C = new PointND(SimpleOperators.subtract(triangle.getC().getAbstractVector(), origin.getAbstractVector()));
  					
  					// Perform shear and scale of vertices
  					double[] aCoords = new double[A.getDimension()-1];
  					double[] bCoords = new double[B.getDimension()-1];
  					double[] cCoords = new double[C.getDimension()-1];
  					
  					aCoords[0] = A.get(kx) - Sx * A.get(kz);
  					aCoords[1] = A.get(ky) - Sy * A.get(kz);
  					bCoords[0] = B.get(kx) - Sx * B.get(kz);
  					bCoords[1] = B.get(ky) - Sy * B.get(kz);
  					cCoords[0] = C.get(kx) - Sx * C.get(kz);
  					cCoords[1] = C.get(ky) - Sy * C.get(kz);
  					
  					// Calculate scaled barycentric coordinates
  					double U = cCoords[0] * bCoords[1] - cCoords[1] * bCoords[0];
  					double V = aCoords[0] * cCoords[1] - aCoords[1] * cCoords[0];
  					double W = bCoords[0] * aCoords[1] - bCoords[1] * aCoords[0];
  					
  					// Perform edge tests
  					// All scaled barycentric coordinates need to have the same sign
  					if ((U < 0.d || V < 0.d || W < 0.d) && (U > 0.d || V > 0.d || W > 0.d)) {
  						continue;
  					}
  					
  					// Calculate determinant
  					double det = U + V + W;
  					
  					if (CONRAD.SMALL_VALUE > Math.abs(det)) {
  						// det is almost zero
  						// The ray seems to be co-planar to the triangle surface
  						continue;
  					}
  					
  					// Calculate scaled z-coordinates of vertices and use them to calculate the hit distance
  					double Az = Sz * A.get(kz);
  					double Bz = Sz * B.get(kz);
  					double Cz = Sz * C.get(kz);
  					double T = U * Az + V * Bz + W * Cz;
  					
  					// Handle negative determinant
  					if (det >= 0 && T < 0.d) {
  						continue;
  					} else if (det < 0 && T >= 0.d) {
  						continue;
  					}
  					
  					// Normalize barycentric coordinates and hit distance
  					U = U / det;
  					V = V / det;
  					W = W / det;
  					T = T / det;
  					
  					// Hit has the coordinates [0, 0, T]
  					// Inverting the transformation which was applied to all vertices in the beginning yields:
  					double[] hitCoords = new double[3];
  					hitCoords[0] = T * dir[0] + origin.get(0);
  					hitCoords[1] = T * dir[1] + origin.get(1);
  					hitCoords[2] = T * dir[2] + origin.get(2);
  					
  					PhysicalPoint hit = new PhysicalPoint(hitCoords);
  					// Calculating the dot product of the triangle's normal and the ray direction can help to determine whether the ray entered or left the object at this hit point
  					hit.setHitOrientation(SimpleOperators.multiplyInnerProd(triangle.getNormal(), line.getDirection()));
  					hit.setObject(obj);
  					hits.add(hit);
  				}
  			}
  			// Cleanup for next physical object
  			queue.clear();
  		}
  		
  		if (hits == null || hits.size() == 0) {
  			return null;
  		}
  		
  		return hits;
  
  	}
  
  }
  