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;
031
032import java.util.ArrayList;
033import java.util.List;
034
035import org.openimaj.citation.annotation.Reference;
036import org.openimaj.citation.annotation.ReferenceType;
037import org.openimaj.citation.annotation.References;
038import org.openimaj.math.geometry.point.Point2d;
039import org.openimaj.math.geometry.point.Point2dImpl;
040import org.openimaj.math.geometry.point.PointList;
041import org.openimaj.math.geometry.shape.algorithm.GeneralisedProcrustesAnalysis;
042import org.openimaj.math.geometry.shape.algorithm.ProcrustesAnalysis;
043import org.openimaj.math.matrix.algorithm.pca.PrincipalComponentAnalysis;
044import org.openimaj.math.matrix.algorithm.pca.PrincipalComponentAnalysis.ComponentSelector;
045import org.openimaj.math.matrix.algorithm.pca.SvdPrincipalComponentAnalysis;
046import org.openimaj.util.pair.IndependentPair;
047
048import Jama.Matrix;
049
050/**
051 * A 2d point distribution model learnt from a set of {@link PointList}s with
052 * corresponding points (the ith point in each {@link PointList} is the same
053 * landmark).
054 *
055 * The pdm models the mean shape and the variance from the mean of the top N
056 * principal components. The model is generative and can generate new shapes
057 * from a scaling vector. To ensure that newly generated shapes are plausible,
058 * scaling vectors have {@link Constraint}s applied to them.
059 *
060 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
061 *
062 */
063@References(references = {
064                @Reference(
065                                author = { "Cootes, T. F.", "Taylor, C. J." },
066                                title = "Statistical Models of Appearance for Computer Vision",
067                                type = ReferenceType.Unpublished,
068                                month = "October",
069                                year = "2001",
070                                url = "http://isbe.man.ac.uk/~bim/Models/app_model.ps.gz"
071                ),
072                @Reference(
073                                type = ReferenceType.Inproceedings,
074                                author = { "C. J. Taylor", "D. H. Cooper", "J. Graham" },
075                                title = "Training models of shape from sets of examples",
076                                year = "1992",
077                                booktitle = "Proc. BMVC92, Springer-Verlag",
078                                pages = { "9", "", "18" }
079                )
080})
081public class PointDistributionModel {
082        /**
083         * Interface for modelling constraints applied to the scaling vector of
084         * {@link PointDistributionModel}s so that generated models are plausible.
085         *
086         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
087         */
088        public interface Constraint {
089                /**
090                 * Apply constraints to a scaling vector so that it will generated a
091                 * plausible model and return the new constrained vector.
092                 *
093                 * @param scaling
094                 *            the scaling vector to constrain
095                 * @param lamda
096                 *            the eigenvalues of the {@link PointDistributionModel}
097                 * @return the constrained scaling vector
098                 */
099                public double[] apply(double[] scaling, double[] lamda);
100        }
101
102        /**
103         * A constraint that does nothing.
104         *
105         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
106         *
107         */
108        public static class NullConstraint implements Constraint {
109                @Override
110                public double[] apply(double[] in, double[] lamda) {
111                        return in;
112                }
113        }
114
115        /**
116         * A constraint that ensures that each individual element of the scaling
117         * vector is within +/- x standard deviations of the model.
118         *
119         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
120         *
121         */
122        public static class BoxConstraint implements Constraint {
123                double multiplier;
124
125                /**
126                 * Construct with the given multiplier of the standard deviation.
127                 *
128                 * @param multiplier
129                 */
130                public BoxConstraint(double multiplier) {
131                        this.multiplier = multiplier;
132                }
133
134                @Override
135                public double[] apply(double[] in, double[] lamda) {
136                        final double[] out = new double[in.length];
137
138                        for (int i = 0; i < in.length; i++) {
139                                final double w = multiplier * Math.sqrt(lamda[i]);
140                                out[i] = in[i] > w ? w : in[i] < -w ? -w : in[i];
141                        }
142
143                        return out;
144                }
145        }
146
147        /**
148         * Constrain the scaling vector to a hyper-ellipsoid.
149         *
150         * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
151         *
152         */
153        public static class EllipsoidConstraint implements Constraint {
154                double dmax;
155
156                /**
157                 * Construct with the given maximum normalised ellipsoid radius.
158                 *
159                 * @param dmax
160                 */
161                public EllipsoidConstraint(double dmax) {
162                        this.dmax = dmax;
163                }
164
165                @Override
166                public double[] apply(double[] in, double[] lamda) {
167                        double dmsq = 0;
168                        for (int i = 0; i < in.length; i++) {
169                                dmsq += in[i] * in[i] / lamda[i];
170                        }
171
172                        if (dmsq < dmax * dmax) {
173                                return in;
174                        }
175
176                        final double sc = dmax / Math.sqrt(dmsq);
177                        final double[] out = new double[in.length];
178                        for (int i = 0; i < in.length; i++) {
179                                out[i] = in[i] * sc;
180                        }
181
182                        return out;
183                }
184        }
185
186        protected Constraint constraint;
187        protected PrincipalComponentAnalysis pc;
188        protected PointList mean;
189        protected int numComponents;
190        protected int maxIter = 100;
191
192        /**
193         * Construct a {@link PointDistributionModel} from the given data with a
194         * {@link NullConstraint}.
195         *
196         * @param data
197         */
198        public PointDistributionModel(List<PointList> data) {
199                this(new NullConstraint(), data);
200        }
201
202        /**
203         * Construct a {@link PointDistributionModel} from the given data and
204         * {@link Constraint}.
205         *
206         * @param constraint
207         * @param data
208         */
209        public PointDistributionModel(Constraint constraint, List<PointList> data) {
210                this.constraint = constraint;
211
212                // align
213                mean = GeneralisedProcrustesAnalysis.alignPoints(data, 5, 10);
214
215                // build data matrix
216                final Matrix m = buildDataMatrix(data);
217
218                // perform pca
219                this.pc = new SvdPrincipalComponentAnalysis();
220                pc.learnBasis(m);
221
222                numComponents = this.pc.getEigenValues().length;
223        }
224
225        private Matrix buildDataMatrix(PointList data) {
226                final List<PointList> pls = new ArrayList<PointList>(1);
227                pls.add(data);
228                return buildDataMatrix(pls);
229        }
230
231        private Matrix buildDataMatrix(List<PointList> data) {
232                final int nData = data.size();
233                final int nPoints = data.get(0).size();
234
235                final Matrix m = new Matrix(nData, nPoints * 2);
236                final double[][] mData = m.getArray();
237
238                for (int i = 0; i < nData; i++) {
239                        final PointList pts = data.get(i);
240                        for (int j = 0, k = 0; k < nPoints; j += 2, k++) {
241                                final Point2d pt = pts.points.get(k);
242
243                                mData[i][j] = pt.getX();
244                                mData[i][j + 1] = pt.getY();
245                        }
246                }
247
248                return m;
249        }
250
251        /**
252         * @return the mean shape
253         */
254        public PointList getMean() {
255                return mean;
256        }
257
258        /**
259         * Set the number of components of the PDM
260         *
261         * @param n
262         *            number of components
263         */
264        public void setNumComponents(int n) {
265                pc.selectSubset(n);
266                numComponents = this.pc.getEigenValues().length;
267        }
268
269        /**
270         * Set the number of components of the PDM using a {@link ComponentSelector}
271         * .
272         *
273         * @param selector
274         *            the {@link ComponentSelector} to apply.
275         */
276        public void setNumComponents(ComponentSelector selector) {
277                pc.selectSubset(selector);
278                numComponents = this.pc.getEigenValues().length;
279        }
280
281        /**
282         * Generate a plausible new shape from the scaling vector. The scaling
283         * vector is constrained by the underlying {@link Constraint} before being
284         * used to generate the model.
285         *
286         * @param scaling
287         *            scaling vector.
288         * @return a new shape
289         */
290        public PointList generateNewShape(double[] scaling) {
291                final PointList newShape = new PointList();
292
293                final double[] pts = pc.generate(constraint.apply(scaling, pc.getEigenValues()));
294
295                for (int i = 0; i < pts.length; i += 2) {
296                        final float x = (float) pts[i];
297                        final float y = (float) pts[i + 1];
298
299                        newShape.points.add(new Point2dImpl(x, y));
300                }
301
302                return newShape;
303        }
304
305        /**
306         * Compute the standard deviations of the shape components, multiplied by
307         * the given value.
308         *
309         * @param multiplier
310         *            the multiplier
311         * @return the multiplied standard deviations
312         */
313        public double[] getStandardDeviations(double multiplier) {
314                final double[] rngs = pc.getStandardDeviations();
315
316                for (int i = 0; i < rngs.length; i++) {
317                        rngs[i] = rngs[i] * multiplier;
318                }
319
320                return rngs;
321        }
322
323        /**
324         * Determine the best parameters of the PDM for the given model.
325         *
326         * @param observed
327         *            the observed model.
328         * @return the parameters that best fit the model.
329         */
330        public IndependentPair<Matrix, double[]> fitModel(PointList observed) {
331                double[] model = new double[numComponents];
332                double delta = 1.0;
333                Matrix pose = null;
334
335                final ProcrustesAnalysis pa = new ProcrustesAnalysis(observed);
336                int count = 0;
337                while (delta > 1e-6 && count++ < maxIter) {
338                        final PointList instance = this.generateNewShape(model);
339
340                        pose = pa.align(instance);
341
342                        final PointList projected = observed.transform(pose.inverse());
343
344                        // TODO: tangent space projection???
345
346                        final Matrix y = buildDataMatrix(projected);
347                        final Matrix xbar = new Matrix(new double[][] { pc.getMean() });
348                        double[] newModel = (y.minus(xbar)).times(pc.getBasis()).getArray()[0];
349
350                        newModel = constraint.apply(newModel, pc.getEigenValues());
351
352                        delta = 0;
353                        for (int i = 0; i < newModel.length; i++)
354                                delta += (newModel[i] - model[i]) * (newModel[i] - model[i]);
355                        delta = Math.sqrt(delta);
356
357                        model = newModel;
358                }
359
360                return new IndependentPair<Matrix, double[]>(pose, model);
361        }
362}