package cz.fidentis.analyst.procrustes; import cz.fidentis.analyst.face.HumanFace; import cz.fidentis.analyst.feature.FeaturePoint; import cz.fidentis.analyst.icp.IcpTransformation; import cz.fidentis.analyst.mesh.core.MeshPoint; import cz.fidentis.analyst.procrustes.exceptions.ProcrustesAnalysisException; import cz.fidentis.analyst.procrustes.utils.ProcrustesAnalysisUtils; import org.ejml.data.Matrix; import org.ejml.simple.SimpleMatrix; import org.ejml.simple.SimpleSVD; import java.util.ArrayList; import java.util.HashMap; import java.util.List; import javax.vecmath.Point3d; /** * @author Jakub Kolman */ public class ProcrustesAnalysis { private ProcrustesAnalysisFaceModel faceModel1; private ProcrustesAnalysisFaceModel faceModel2; private Point3d modelCentroid1; private Point3d modelCentroid2; // based on some example models this is an approximate centroid of feature points // private Point3d centroid = new Point3d(0, 0, 0); /** * Constructor * * @param humanFace1 * @param humanFace2 * @throws ProcrustesAnalysisException */ public ProcrustesAnalysis( HumanFace humanFace1, HumanFace humanFace2) throws ProcrustesAnalysisException { ProcrustesAnalysisFaceModel model1 = new ProcrustesAnalysisFaceModel(humanFace1); ProcrustesAnalysisFaceModel model2 = new ProcrustesAnalysisFaceModel(humanFace2); if (model1.getFeaturePointsMap().values().size() != model2.getFeaturePointsMap().values().size()) { throw new ProcrustesAnalysisException("Lists of feature points do not have the same size"); } // this.orderedFeaturePointList1 = new ArrayList<>(model1.getFeaturePointsMap().values()); // this.orderedFeaturePointList2 = new ArrayList<>(model2.getFeaturePointsMap().values()); if (!featurePointTypesEquivalence( model1.getFeaturePointValues(), model2.getFeaturePointValues())) { throw new ProcrustesAnalysisException("Lists of feature points do not have the same feature point types"); } this.faceModel1 = model1; this.faceModel2 = model2; this.modelCentroid1 = ProcrustesAnalysisUtils.findCentroidOfFeaturePoints(model1.getFeaturePointValues()); this.modelCentroid2 = ProcrustesAnalysisUtils.findCentroidOfFeaturePoints(model2.getFeaturePointValues()); } public ProcrustesAnalysisFaceModel getFaceModel1() { return faceModel1; } public ProcrustesAnalysisFaceModel getFaceModel2() { return faceModel2; } public void analyze() { List<SimpleMatrix> transformation = new ArrayList<>(); this.superImpose(); this.rotate(); } private Point3d computeCentroidsDistance(Point3d centroid1, Point3d centroid2) { double x = (centroid1.x - centroid2.x); double y = (centroid1.y - centroid2.y); double z = (centroid1.z - centroid2.z); Point3d computedCentroid = new Point3d(x, y, z); return computedCentroid; } /** * Imposes two face models (lists of feature points and vertices) over each other */ private void superImpose() { // moves vertices of face models // moveModelToSameOrigin(this.faceModel2.getVertices(), this.faceModel2.getFeaturePointsMap()); normalize(this.faceModel1, this.modelCentroid1); normalize(this.faceModel2, this.modelCentroid2); // moveVertices(this.modelCentroid1, featurePointModel1.getVertices()); // moveModelToPoint(this.modelCentroid1, faceModel1.getVertices(), faceModel1.getFeaturePointsMap()); // moveModelToPoint(this.modelCentroid1, faceModel2.getVertices(), faceModel2.getFeaturePointsMap()); } private void normalize(ProcrustesAnalysisFaceModel faceModel, Point3d centroid) { // float size = ProcrustesAnalysisUtils.countSize(centroid, faceModel.getFeaturePointsMap()); centerToOrigin(faceModel, centroid); } private void centerToOrigin(ProcrustesAnalysisFaceModel faceModel, Point3d centroid) { for (FeaturePoint fp: faceModel.getFeaturePointsMap().values()) { fp.getPosition().x -= centroid.x; fp.getPosition().y -= centroid.y; fp.getPosition().z -= centroid.z; } SimpleMatrix centeredVertices = new SimpleMatrix( faceModel.getVerticesMatrix().numRows(), faceModel.getVerticesMatrix().numCols()); for (int i = 0; i < faceModel.getVerticesMatrix().numRows(); i++) { centeredVertices.set(i, 0, faceModel.getVerticesMatrix().get(i, 0) - centroid.x); centeredVertices.set(i, 1, faceModel.getVerticesMatrix().get(i, 1) - centroid.y); centeredVertices.set(i, 2, faceModel.getVerticesMatrix().get(i, 2) - centroid.z); } faceModel.setVerticesMatrix(centeredVertices); } private void moveModelToSameOrigin(List<MeshPoint> vertices, HashMap<Integer, FeaturePoint> featurePointsMap) { Point3d diff = computeCentroidsDistance(this.modelCentroid1, this.modelCentroid2); moveModelToPoint(diff, vertices, featurePointsMap); } private void moveModelToPoint(Point3d centroid, List<MeshPoint> vertices, HashMap<Integer, FeaturePoint> featurePointsMap) { moveVertices(centroid, vertices); moveFeaturePointsToVertices(centroid, featurePointsMap); } private void moveFeaturePointsToVertices(Point3d centroid, HashMap<Integer, FeaturePoint> featurePointsMap) { for (FeaturePoint fp : featurePointsMap.values()) { fp.getPosition().x = fp.getPosition().x + centroid.x; fp.getPosition().y = fp.getPosition().y + centroid.y; fp.getPosition().z = fp.getPosition().z + centroid.z; } } /** * By rotation of matrices solves orthogonal procrustes problem */ private void rotate() { // There is no reason trying to rotate less than 3 elements if (this.getFaceModel1().getFeaturePointsMap().size() < 3) { throw new ProcrustesAnalysisException("To do procrustes analysis models have to have at least 3 feature points"); } SimpleMatrix primaryMatrix = ProcrustesAnalysisUtils.createFeaturePointMatrix(this.faceModel2.getFeaturePointValues()); SimpleMatrix transposedMatrix = ProcrustesAnalysisUtils.createFeaturePointMatrix( this.faceModel1.getFeaturePointValues()).transpose(); SimpleMatrix svdMatrix = transposedMatrix.mult(primaryMatrix); SimpleSVD<SimpleMatrix> singularValueDecomposition = svdMatrix.svd(); SimpleMatrix transposedU = singularValueDecomposition.getU().transpose(); SimpleMatrix r = singularValueDecomposition.getV().mult(transposedU); primaryMatrix = primaryMatrix.mult(r); this.faceModel2.setFeaturePointsMap( ProcrustesAnalysisUtils.createFeaturePointMapFromMatrix( primaryMatrix, this.faceModel2)); // for (int i = 0; i < this.faceModel2.getFeaturePointTypeCorrespondence().size(); i++) { // FeaturePoint fp = new FeaturePoint( // primaryMatrix.get(i, 0), // primaryMatrix.get(i, 1), // primaryMatrix.get(i, 2), // this.faceModel2.getFeaturePointsMap().get( // this.faceModel2.getFeaturePointTypeCorrespondence().get(i)).getFeaturePointType() // ); // this.faceModel2.getFeaturePointsMap().put(i, fp); // } // rotateVertices(this.faceModel2, primaryMatrix); SimpleMatrix rotatedVertices = this.faceModel2.getVerticesMatrix().mult(r); createListFromMatrix(rotatedVertices, this.faceModel2); } private void createListFromMatrix(SimpleMatrix rotatedVertices, ProcrustesAnalysisFaceModel faceModel) { if (faceModel.getVertices() != null) { for (int i = 0; i < rotatedVertices.numRows(); i++) { faceModel.getVertices().get(i).getPosition().x = rotatedVertices.get(i, 0); faceModel.getVertices().get(i).getPosition().y = rotatedVertices.get(i, 1); faceModel.getVertices().get(i).getPosition().z = rotatedVertices.get(i, 2); } } } /** * Rotates every vertex of model by given matrix * * @param model * @param matrix */ // if rotated vertices are drawn immediately it is better to set them after rotating them all // so it would be drawn just once // private void rotateVertices(ProcrustesAnalysisFaceModel model, SimpleMatrix matrix) { // if (model.getVertices() != null) { // for (MeshPoint v : model.getVertices()) { // v.setPosition(ProcrustesAnalysisUtils.rotateVertex(v, matrix)); // } // } // } /** * Calculates new vertices adjusted to the centroid by method {@see findCenteroidOfFeaturePoints}, * so moved featurepoints would correspond with the same place on face as they did before. * * @param centroid * @param vertices * @return */ private void moveVertices(Point3d centroid, List<MeshPoint> vertices) { if (vertices != null && centroid != null) { // List<MeshPoint> movedVertices = null; for (MeshPoint v : vertices) { v.getPosition().x = (v.getPosition().x + centroid.x); v.getPosition().y = (v.getPosition().y + centroid.y); v.getPosition().z = (v.getPosition().z + centroid.z); } } else { throw new ProcrustesAnalysisException("Could not compute vertices locations after moving the centroid from model to feature points"); } } /** * Creates sorted feature point lists and compares them whether they have * the same feature point types. * * @param fpList1 * @param fpList2 * @return */ private boolean featurePointTypesEquivalence(List<FeaturePoint> fpList1, List<FeaturePoint> fpList2) { return ProcrustesAnalysisUtils.checkFeaturePointsType(fpList1, fpList2); } /** * Calculate scaling ratio of how much the appropriate object corresponding * to the second feature point list has to be scale up or shrunk. * <p> * If returned ratioValue is greater 1 then it means that the second object * should be scaled up ratioValue times. If returned ratioValue is smaller 1 * than the second object should be shrunk. * * @return ratioValue */ // private double calculateScalingValue() { // double[] distancesOfList1 = ProcrustesAnalysisUtils.calculateMeanDistancesFromOrigin( // this.faceModel1.getFeaturePointValues()); // double[] distancesOfList2 = ProcrustesAnalysisUtils.calculateMeanDistancesFromOrigin( // this.faceModel2.getFeaturePointValues()); // // double[] ratioArray = new double[distancesOfList1.length]; // double ratioValue = 0; // // for (int i = 0; i < distancesOfList1.length; i++) { // ratioArray[i] += distancesOfList1[i] / distancesOfList2[i]; // } // // for (double ratio : ratioArray) { // ratioValue += ratio; // } // // return ratioValue / distancesOfList1.length; // } }