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 org.openimaj.math.statistics.distribution.CachingMultivariateGaussian; 033import org.openimaj.math.statistics.distribution.MultivariateGaussian; 034import org.openimaj.math.util.QuadraticEquation; 035 036import Jama.EigenvalueDecomposition; 037import Jama.Matrix; 038 039/** 040 * An elliptical shape 041 * 042 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 043 * 044 */ 045public class EllipseUtilities { 046 /*** 047 * Construct an ellipse using a parametric ellipse equation, namely: 048 * 049 * X(t) = centerX + major * cos(t) * cos(rotation) - minor * sin(t) * 050 * sin(rotation) Y(t) = centerY + major * cos(t) * cos(rotation) + minor * 051 * sin(t) * sin(rotation) 052 * 053 * @param centerX 054 * @param centerY 055 * @param major 056 * @param minor 057 * @param rotation 058 * @return an ellipse 059 */ 060 public static Ellipse 061 ellipseFromEquation(double centerX, double centerY, double major, double minor, double rotation) 062 { 063 return new Ellipse(centerX, centerY, major, minor, rotation); 064 } 065 066 /** 067 * Construct ellipse from second moment matrix and centroid. 068 * 069 * @param x 070 * x-ordinate of centroid 071 * @param y 072 * y-ordinate of centroid 073 * @param secondMoments 074 * second moments matrix 075 * @return an ellipse 076 */ 077 public static Ellipse ellipseFromSecondMoments(float x, float y, Matrix secondMoments) { 078 return EllipseUtilities.ellipseFromSecondMoments(x, y, secondMoments, 1); 079 } 080 081 /** 082 * Construct ellipse from second moment matrix, scale-factor and centroid. 083 * 084 * @param x 085 * x-ordinate of centroid 086 * @param y 087 * y-ordinate of centroid 088 * @param secondMoments 089 * second moments matrix 090 * @param scaleFactor 091 * the scale factor 092 * @return an ellipse 093 */ 094 public static Ellipse ellipseFromSecondMoments(float x, float y, Matrix secondMoments, double scaleFactor) { 095 final double divFactor = 1 / Math.sqrt(secondMoments.det()); 096 final EigenvalueDecomposition rdr = secondMoments.times(divFactor).eig(); 097 double d1, d2; 098 if (rdr.getD().get(0, 0) == 0) 099 d1 = 0; 100 else 101 d1 = 1.0 / Math.sqrt(rdr.getD().get(0, 0)); 102 if (rdr.getD().get(1, 1) == 0) 103 d2 = 0; 104 else 105 d2 = 1.0 / Math.sqrt(rdr.getD().get(1, 1)); 106 107 final double scaleCorrectedD1 = d1 * scaleFactor; 108 final double scaleCorrectedD2 = d2 * scaleFactor; 109 110 final Matrix eigenMatrix = rdr.getV(); 111 112 final double rotation = Math.atan2(eigenMatrix.get(1, 0), eigenMatrix.get(0, 0)); 113 return ellipseFromEquation(x, y, scaleCorrectedD1, scaleCorrectedD2, rotation); 114 } 115 116 /** 117 * Construct ellipse from covariance matrix, scale-factor and centroid. 118 * 119 * @param x 120 * x-ordinate of centroid 121 * @param y 122 * y-ordinate of centroid 123 * @param sm 124 * covariance matrix 125 * @param sf 126 * scale-factor 127 * @return an ellipse 128 */ 129 public static Ellipse ellipseFromCovariance(float x, float y, Matrix sm, float sf) { 130 final double xy = sm.get(1, 0); 131 final double xx = sm.get(0, 0); 132 final double yy = sm.get(1, 1); 133 final double theta = 0.5 * Math.atan2(2 * xy, xx - yy); 134 135 final double trace = xx + yy; 136 final double det = (xx * yy) - (xy * xy); 137 final double[] eigval = QuadraticEquation.solveGeneralQuadratic(1, -trace, det); 138 139 final double a = Math.sqrt(eigval[1]) * sf; 140 final double b = Math.sqrt(eigval[0]) * sf; 141 return ellipseFromEquation(x, y, a, b, theta); 142 } 143 144 /** 145 * Create the covariance matrix of an ellipse. 146 * 147 * @param e 148 * the ellipse 149 * @return the corresponding covariance matrix 150 */ 151 public static Matrix ellipseToCovariance(Ellipse e) { 152 final Matrix transform = e.transformMatrix(); 153 final Matrix Q = transform.getMatrix(0, 1, 0, 1); 154 return Q.times(Q.transpose()); 155 // double sinrot = Math.sin(e.getRotation()); 156 // double cosrot = Math.cos(e.getRotation()); 157 // double cosrot2 = cosrot * cosrot; 158 // double sinrot2 = sinrot * sinrot; 159 // double a2 = e.getMajor() * e.getMajor(); 160 // double b2 = e.getMinor() * e.getMinor(); 161 // Matrix Q = new Matrix(new double[][]{ 162 // {cosrot2 / a2 + sinrot2 / b2 , sinrot*cosrot*(1/a2 - 1/b2)}, 163 // {sinrot*cosrot*(1/a2 - 1/b2) , sinrot2 / a2 + cosrot2 / b2} 164 // }); 165 // return Q.inverse(); 166 } 167 168 /** 169 * Create an ellipse. 170 * 171 * @param U 172 * @param x 173 * @param y 174 * @param scale 175 * @return an ellipse 176 */ 177 public static Ellipse fromTransformMatrix2x2(Matrix U, float x, float y, float scale) { 178 Matrix uVal, uVec; 179 final EigenvalueDecomposition ueig = U.eig(); 180 uVal = ueig.getD(); 181 uVec = ueig.getV(); 182 183 // Normalize min eigenvalue to 1 to expand patch in the direction of min 184 // eigenvalue of U.inv() 185 final double uval1 = uVal.get(0, 0); 186 final double uval2 = uVal.get(1, 1); 187 188 if (Math.abs(uval1) < Math.abs(uval2)) 189 { 190 uVal.set(0, 0, 1); 191 uVal.set(1, 1, uval2 / uval1); 192 193 } else 194 { 195 uVal.set(1, 1, 1); 196 uVal.set(0, 0, uval1 / uval2); 197 } 198 199 final float ax1 = (float) (1 / Math.abs(uVal.get(1, 1)) * scale); 200 final float ax2 = (float) (1 / Math.abs(uVal.get(0, 0)) * scale); 201 final double phi = Math.atan(uVec.get(1, 1) / uVec.get(0, 1)); 202 203 return new Ellipse(x, y, ax1, ax2, phi); 204 } 205 206 /** 207 * Construct an ellipse that encompasses the shape of a 208 * {@link CachingMultivariateGaussian}. 209 * 210 * @param gaussian 211 * the {@link CachingMultivariateGaussian} 212 * @param scale 213 * the relative size of the ellipse 214 * @return the ellipse 215 */ 216 public static Ellipse ellipseFromGaussian(MultivariateGaussian gaussian, float scale) { 217 final Matrix mean = gaussian.getMean(); 218 final float x = (float) mean.get(0, 0); 219 final float y = (float) mean.get(0, 1); 220 final Matrix covar = gaussian.getCovariance(); 221 return ellipseFromCovariance(x, y, covar, scale); 222 } 223}