001/**
002 * Copyright (c) 2011, The University of Southampton and the individual contributors.
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without modification,
006 * are permitted provided that the following conditions are met:
007 *
008 *   *  Redistributions of source code must retain the above copyright notice,
009 *      this list of conditions and the following disclaimer.
010 *
011 *   *  Redistributions in binary form must reproduce the above copyright notice,
012 *      this list of conditions and the following disclaimer in the documentation
013 *      and/or other materials provided with the distribution.
014 *
015 *   *  Neither the name of the University of Southampton nor the names of its
016 *      contributors may be used to endorse or promote products derived from this
017 *      software without specific prior written permission.
018 *
019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
029 */
030package org.openimaj.math.geometry.shape.algorithm;
031
032import java.util.ArrayList;
033import java.util.List;
034
035import org.openimaj.math.geometry.point.Point2d;
036import org.openimaj.math.geometry.point.PointList;
037import org.openimaj.math.geometry.transforms.TransformUtilities;
038
039import Jama.Matrix;
040
041/**
042 * Implementation of shape alignment using Generalised Procrustes Analysis.
043 * 
044 * The alignment process is iterative and stops after a given number
045 * of iterations or when the Procrustes Distance between the current mean and previous
046 * reference shape is less than a threshold (i.e. the rate of change is small).
047 * 
048 * All shapes are aligned inplace. The reference shape is optionally normalised
049 * to a standardised scale and translated to the origin.
050 *  
051 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
052 */
053public class GeneralisedProcrustesAnalysis {
054        
055        private static final int DEFAULT_MAX_ITERS = 10;
056        private float threshold;
057        private int maxIters;
058
059        /**
060         * Construct the {@link GeneralisedProcrustesAnalysis} with the given
061         * parameters. The alignment process is 
062         * iterative and stops after a given number of iterations (last parameter)
063         * or when the Procrustes Distance between the current mean and previous
064         * reference shape is less than a threshold (i.e. the rate of change is small) (first parameter).
065         * 
066         * @param threshold the threshold on the Procrustes Distance at which to stop iterating
067         * @param maxIters the maximum number of iterations.
068         */
069        public GeneralisedProcrustesAnalysis(float threshold, int maxIters) {
070                this.threshold = threshold;
071                this.maxIters = maxIters;
072        }
073        
074        /**
075         * Construct the {@link GeneralisedProcrustesAnalysis} with the given
076         * parameters. The alignment process is 
077         * iterative and stops after a default number of iterations (10)
078         * or when the Procrustes Distance between the current mean and previous
079         * reference shape is less than a threshold (i.e. the rate of change is small)
080         * 
081         * @param threshold the threshold on the Procrustes Distance at which to stop iterating
082         */
083        public GeneralisedProcrustesAnalysis(float threshold) {
084                this(threshold, DEFAULT_MAX_ITERS);
085        }
086        
087        /**
088         * Align the input shapes to the "mean" shape. All shapes are aligned inplace.
089         * The mean shape is returned. 
090         * 
091         * @param shapes The shapes to align 
092         * @return the mean shape.
093         */
094        public PointList align(List<PointList> shapes) {
095                return alignPoints(shapes, threshold, maxIters);
096        }
097        
098        /**
099         * Align the input shapes to the "mean" shape using Generalised Procrustes Analysis.
100         * The alignment process is iterative and stops after a given number
101         * of iterations or when the Procrustes Distance between the current mean and previous
102         * reference shape is less than a threshold (i.e. the rate of change is small).
103         * 
104         * All shapes are aligned inplace. The reference shape is initially normalised
105         * to a standardised scale and translated to the origin. The mean shape is returned.
106         * 
107         * @param inputShapes The shapes to align
108         * @param threshold the threshold on the Procrustes Distance at which to stop iterating
109         * @param maxIters the maximum number of iterations.
110         * @return the mean shape
111         */
112        public static PointList alignPoints(List<PointList> inputShapes, float threshold, int maxIters) {
113                PointList reference = inputShapes.get(0); 
114                
115                List<PointList> workingShapes = new ArrayList<PointList>(inputShapes);  
116                workingShapes.remove(reference);
117                
118                //Use a ProcrustesAnalysis to get the reference scaling and cog (normalised if required)
119                ProcrustesAnalysis referencePa = new ProcrustesAnalysis(reference, true);
120                double referenceScaling = referencePa.scaling;
121                Point2d referenceCog = referencePa.referenceCog;
122        
123                if (workingShapes.size() == 0)
124                        return reference;
125                
126                PointList mean = alignPointsAndAverage(workingShapes, reference, referenceScaling, referenceCog);
127                
128                while (ProcrustesAnalysis.computeProcrustesDistance(reference, mean) > threshold && maxIters-- >= 0) {
129                        reference = mean;
130                        mean = alignPointsAndAverage(inputShapes, reference, referenceScaling, referenceCog);
131                }
132                
133                return mean;
134        }
135        
136        protected static PointList alignPointsAndAverage(List<PointList> shapes, PointList reference, double referenceScaling, Point2d referenceCog) {
137                ProcrustesAnalysis pa = new ProcrustesAnalysis(reference);
138                
139                for (PointList shape : shapes) {
140                        pa.align(shape);
141                }
142                
143                PointList mean = PointList.computeMean(shapes);
144                
145                //normalise translation to reference
146                Point2d cog = mean.calculateCentroid();
147                Matrix trans = TransformUtilities.translateToPointMatrix(cog, referenceCog);
148                mean.translate((float)trans.get(0,2), (float)trans.get(1,2));
149                                
150                //normalise scaling to reference
151                double scale = ProcrustesAnalysis.computeScale(mean, referenceCog.getX(), referenceCog.getY());
152                float sf = (float)(scale / referenceScaling);
153                mean.scale(referenceCog, sf);
154                
155                return mean;
156        }
157        
158}