package cz.fidentis.analyst.visitors.mesh;

import cz.fidentis.analyst.mesh.MeshVisitor;
import cz.fidentis.analyst.mesh.core.MeshFacet;
import cz.fidentis.analyst.mesh.core.MeshTriangle;
import cz.fidentis.analyst.mesh.core.TriangleFan;
import java.util.ArrayList;
import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Objects;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;

/**
 * Abstract class for algorithms calculating curvatures of mesh vertices.
 * @see https://computergraphics.stackexchange.com/questions/1718/what-is-the-simplest-way-to-compute-principal-curvature-for-a-mesh-triangle
 * @see http://rodolphe-vaillant.fr/?e=20
 * <p>
 * <b>All curvature algorithms suppose that the triangle vertices are oriented clockwise!</b>
 * </p>
 * 
 * @author Natalia Bebjakova
 * @author Radek Oslejsek
 */
public class Curvature extends MeshVisitor {
    
    private final Map<MeshFacet, List<Double>> gaussian = new HashMap<>();
    private final Map<MeshFacet, List<Double>> mean = new HashMap<>();
    private final Map<MeshFacet, List<Double>> minPrincipal = new HashMap<>();
    private final Map<MeshFacet, List<Double>> maxPrincipal = new HashMap<>();
    
    private final boolean parallel;
    
    /**
     * Constructor.
     * 
     * @param parallel If {@code true}, then the algorithm runs concurrently utilizing all CPU cores
     */
    public Curvature(boolean parallel) {
        this.parallel = parallel;
    }
    
    /**
     * Constructor for parallel computation.
     * 
     */
    public Curvature() {
        this(true);
    }
    
    /**
     * Returns Gaussian curvatures for all inspected mesh facets. The order corresponds to
     * the order of vertices, i.e., i-th value represents the curvature of i-th mesh vertex.
     * 
     * @return Gaussian curvatures of inspected mesh facets.
     */
    public Map<MeshFacet, List<Double>> getGaussianCurvatures() {
        return Collections.unmodifiableMap(gaussian);
    }
    
    /**
     * Returns mean curvatures for all inspected mesh facets. The order corresponds to
     * the order of vertices, i.e., i-th value represents the curvature of i-th mesh vertex.
     * 
     * @return Mean curvatures of inspected mesh facets.
     */
    public Map<MeshFacet, List<Double>> getMeanCurvatures() {
        return Collections.unmodifiableMap(mean);
    }
    
    /**
     * Returns minimum principal curvatures for all inspected mesh facets. The order corresponds to
     * the order of vertices, i.e., i-th value represents the curvature of i-th mesh vertex.
     * 
     * @return Minimum principal curvatures of inspected mesh facets.
     */
    public Map<MeshFacet, List<Double>> getMinPrincipalCurvatures() {
        return Collections.unmodifiableMap(minPrincipal);
    }
    
    /**
     * Returns maximum principal curvatures for all inspected mesh facets. The order corresponds to
     * the order of vertices, i.e., i-th value represents the curvature of i-th mesh vertex.
     * 
     * @return Maximum principal curvatures of inspected mesh facets.
     */
    public Map<MeshFacet, List<Double>> getMaxPrincipalCurvatures() {
        return Collections.unmodifiableMap(maxPrincipal);
    }
    
    @Override
    public void visitMeshFacet(final MeshFacet facet) {
        synchronized (this) {
            if (gaussian.containsKey(facet)) {
                return; // already visited facet
            }
            // prefill the arrays - required for concurrent computation.
            int numVert = facet.getNumberOfVertices();
            gaussian.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
            mean.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
            minPrincipal.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
            maxPrincipal.put(facet, new ArrayList<Double>(Collections.nCopies(numVert, Double.NaN)));
        }
        
        final Cache cache = new Cache(facet);
        
        if (parallel) {
            ExecutorService executor = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors());
            for (int vertIndexA = 0; vertIndexA < facet.getNumberOfVertices(); vertIndexA++) {
                final int index = vertIndexA;
                Runnable worker = new Runnable() {
                    @Override
                    public void run() {
                        computeCurvature(facet, cache, index);
                    }
                };
                executor.execute(worker);
            }
            executor.shutdown();
            while (!executor.isTerminated()){}
        } else {
            for (int vertIndexA = 0; vertIndexA < facet.getNumberOfVertices(); vertIndexA++) {
                computeCurvature(facet, cache, vertIndexA);
            }
        }
    }
    
    protected void computeCurvature(MeshFacet facet, Cache cache, int vertIndexA) {
        TriangleFan oneRing = facet.getOneRingNeighborhood(vertIndexA);

        if (oneRing.isBoundary()) {
            synchronized (this) {
                this.gaussian.get(facet).set(vertIndexA, 0.0);
                this.mean.get(facet).set(vertIndexA, 0.0);
                this.minPrincipal.get(facet).set(vertIndexA, 0.0);
                this.maxPrincipal.get(facet).set(vertIndexA, 0.0);
                return;
            }
        }

        double sumAngles = 0.0;
        double sumArea = 0.0;
        Vector3d pointSum = new Vector3d();

        Point3d vertA = facet.getVertex(vertIndexA).getPosition();
        for (int i = 0; i < oneRing.getVertices().size() - 1; i++) {
            Point3d vertB = facet.getVertex(oneRing.getVertices().get(i)).getPosition();
            Point3d vertC = facet.getVertex(oneRing.getVertices().get(i + 1)).getPosition();

            // Angles:
            sumAngles += cache.getAngle(vertA, vertB, vertC);

            // Area:
            double alpha = cache.getAngle(vertA, vertB, vertC);
            double beta = cache.getAngle(vertB, vertA, vertC);
            double gamma = cache.getAngle(vertC, vertA, vertB);

            Vector3d ab = new Vector3d(vertB);
            ab.sub(vertA);

            Vector3d ac = new Vector3d(vertC);
            ac.sub(vertA);

            double lenSqrtAB = ab.lengthSquared();
            double lenSqrtAC = ac.lengthSquared();

            double piHalf = Math.PI / 2.0; // 90 degrees
            if (alpha >= piHalf || beta >= piHalf || gamma >= piHalf) { // check for obtuse angle
                double area = 0.5 * lenSqrtAB * lenSqrtAC * Math.sin(alpha);
                sumArea += (alpha > piHalf) ? area / 2.0 : area / 4.0;
            } else {
                double cotGamma = cache.getCotan(vertC, vertA, vertB);
                double cotBeta = cache.getCotan(vertB, vertA, vertC);
                sumArea += (lenSqrtAB * cotGamma + lenSqrtAC * cotBeta) / 8.0;
            }

            // Laplace:
            Point3d vertD;
            if (i == (oneRing.getVertices().size() - 2)) { // last triangle
                vertD = facet.getVertex(oneRing.getVertices().get(1)).getPosition();
            } else {
                vertD = facet.getVertex(oneRing.getVertices().get(i + 2)).getPosition();
            }
            //ac.scale(cache2.get(tri).getCotan(vertB) + cache2.get(triNext).getCotan(vertD));
            ac.scale(cache.getCotan(vertB, vertA, vertC) + cache.getCotan(vertD, vertA, vertC));
            pointSum.add(ac);
        }

        double gaussVal = (2.0 * Math.PI - sumAngles) / sumArea;
        double meanVal = 0.25 * sumArea * pointSum.length();
        double delta = Math.max(0, Math.pow(meanVal, 2) - gaussVal);

        synchronized (this) {
            this.gaussian.get(facet).set(vertIndexA, gaussVal);
            this.mean.get(facet).set(vertIndexA, meanVal);
            this.minPrincipal.get(facet).set(vertIndexA, meanVal - Math.sqrt(delta));
            this.maxPrincipal.get(facet).set(vertIndexA, meanVal + Math.sqrt(delta));
        }
    }

    /**
     * Cache key.
     * @author Radek Oslejsek
     */
    private class CacheKey {
        
        private final Point3d vCentral;
        private final Point3d v2;
        private final Point3d v3;
        
        /**
         * Constructor.
         * @param vCentral The vertex for which the values are stored
         * @param v2 Other vertex in the triangle
         * @param v3 Other vertex in the triangle
         */
        CacheKey(Point3d vCentral, Point3d v2, Point3d v3) {
            this.vCentral = vCentral;
            this.v2 = v2;
            this.v3 = v3;
        }

        @Override
        public int hashCode() {
            int hash = 7;
            hash = 11 * hash + Objects.hashCode(this.vCentral);
            hash = 11 * hash + Objects.hashCode(this.v2);
            hash = 11 * hash + Objects.hashCode(this.v3);
            return hash;
        }

        @Override
        public boolean equals(Object obj) {
            if (this == obj) {
                return true;
            }
            if (obj == null) {
                return false;
            }
            if (getClass() != obj.getClass()) {
                return false;
            }
            final CacheKey other = (CacheKey) obj;
            if (!Objects.equals(this.vCentral, other.vCentral)) {
                return false;
            }
            if (!Objects.equals(this.v2, other.v2)) {
                return false;
            }
            if (!Objects.equals(this.v3, other.v3)) {
                return false;
            }
            return true;
        }

        @Override
        public String toString() {
            return vCentral + " | " + v2 + " | " + v3;
        }
        
    }
    
    /**
     * Helper class that caches triangle characteristics used multiples times during the curvature computation.
     * @author Radek Oslejsek
     */
    private class Cache {
        
        private final Map<CacheKey, List<Double>> cache = new HashMap<>();
        
        /**
         * Constructor.
         * @param facet Mesh facet 
         */
        Cache(MeshFacet facet) {
            for (MeshTriangle tri: facet.getTriangles()) {
                computeValues(facet, tri);
            }
        }
        
        /**
         * Returns cached angle in the vertex {@code vCentral} at the triangle
         * containing vertices {@code v2} and {@code v3}. 
         * 
         * @param vCentral Central vertex
         * @param v2 Another vertex in the triangle
         * @param v3 Another vertex in the triangle
         * @return Angle at vertex {@code vCentral}
         */
        public double getAngle(Point3d vCentral, Point3d v2, Point3d v3) {
            CacheKey key = new CacheKey(vCentral, v2, v3);
            if (cache.containsKey(key)) {
                return cache.get(key).get(0);
            } 
            
            key = new CacheKey(vCentral, v3, v2);
            if (cache.containsKey(key)) {
                return cache.get(key).get(0);
            } 
            
            throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]");
        }
        
        /**
         * Returns cached cotangent of the angle in the vertex {@code vCentral} at the triangle
         * containing vertices {@code v2} and {@code v3}. 
         * 
         * @param vCentral Central vertex
         * @param v2 Another vertex in the triangle
         * @param v3 Another vertex in the triangle
         * @return Cotangent of the angle at vertex {@code vCentral}
         */
        public double getCotan(Point3d vCentral, Point3d v2, Point3d v3) {
            CacheKey key = new CacheKey(vCentral, v2, v3);
            if (cache.containsKey(key)) {
                return cache.get(key).get(1);
            }
            
            key = new CacheKey(vCentral, v3, v2);
            if (cache.containsKey(key)) {
                return cache.get(key).get(1);
            } 
            
            throw new IllegalArgumentException("[" + vCentral + ", " + v2 + ", " + v3 + "]");
        }
        
        private void computeValues(MeshFacet facet, MeshTriangle tri) {
            Vector3d ab = new Vector3d(facet.getVertex(tri.index2).getPosition());
            ab.sub(facet.getVertex(tri.index1).getPosition());
            
            Vector3d ac = new Vector3d(facet.getVertex(tri.index3).getPosition());
            ac.sub(facet.getVertex(tri.index1).getPosition());
            
            Vector3d bc = new Vector3d(facet.getVertex(tri.index3).getPosition());
            bc.sub(facet.getVertex(tri.index2).getPosition());
            
            ab.normalize();
            ac.normalize();
            bc.normalize();
            
            double cosA = ab.dot(ac);
            
            ab.scale(-1.0);
            double cosB = ab.dot(bc);
            
            ac.scale(-1.0);
            bc.scale(-1.0);
            double cosC = ac.dot(bc);
            
            List<Double> list1 = new ArrayList<>();
            List<Double> list2 = new ArrayList<>();
            List<Double> list3 = new ArrayList<>();
            
            list1.add(Math.acos(cosA));
            list2.add(Math.acos(cosB));
            list3.add(Math.acos(cosC));
            
            list1.add(1.0 / Math.tan(list1.get(0)));
            list2.add(1.0 / Math.tan(list2.get(0)));
            list3.add(1.0 / Math.tan(list3.get(0)));
            
            cache.put(new CacheKey(tri.getVertex1(), tri.getVertex2(), tri.getVertex3()), list1);
            cache.put(new CacheKey(tri.getVertex2(), tri.getVertex1(), tri.getVertex3()), list2);
            cache.put(new CacheKey(tri.getVertex3(), tri.getVertex2(), tri.getVertex1()), list3);
        }
    }    
}
