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.image.analysis.watershed.feature;
031
032import org.openimaj.image.FImage;
033import org.openimaj.image.pixel.IntValuePixel;
034import org.openimaj.math.geometry.point.Point2d;
035import org.openimaj.math.geometry.shape.Circle;
036import org.openimaj.math.geometry.shape.Ellipse;
037import org.openimaj.math.geometry.shape.EllipseUtilities;
038import org.openimaj.math.geometry.shape.Polygon;
039import org.openimaj.math.util.QuadraticEquation;
040
041import Jama.Matrix;
042
043/**
044 *      Accumulate the values of u11, u20 and u02 required to fit an
045 *      ellipse to the feature.
046 *
047 *      @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
048 *      
049 */
050public class MomentFeature implements ComponentFeature
051{
052        int n = 0;
053        double mx = 0;
054        double my = 0;
055        double Mx2 = 0;
056        double My2 = 0;
057        double sxy = 0;
058        double sx = 0;
059        double sy = 0;
060
061        @Override
062        public void merge(ComponentFeature f) {
063                MomentFeature mf = (MomentFeature) f;
064
065                double dx = mf.mx - mx;
066                double dy = mf.my - my;
067
068                mx = (n*mx + mf.n*mf.mx) / (n + mf.n);
069                my = (n*my + mf.n*mf.my) / (n + mf.n);
070
071                Mx2 += mf.Mx2 + dx*dx*n*mf.n / (n + mf.n);
072                My2 += mf.My2 + dy*dy*n*mf.n / (n + mf.n);
073
074                n += mf.n;
075                sxy += mf.sxy;
076                sx += mf.sx;
077                sy += mf.sy;
078        }
079
080        @Override
081        public void addSample(IntValuePixel p) {
082                n++;
083                double dx = p.x - mx;
084                double dy = p.y - my;
085
086                mx += dx / n;
087                my += dy / n;
088
089                Mx2 += dx * (p.x - mx);
090                My2 += dy * (p.y - my);
091
092                sx += p.x;
093                sy += p.y;
094                sxy += p.x*p.y;
095        }
096
097        /**
098         * Get the number of pixels accumulated
099         * @return the number of pixels
100         */
101        public double n() {
102                return n;
103        }
104
105        /**
106         * Get the value of u11
107         * @return the value of u11
108         */
109        public double u11() {
110                return (sxy - sx*sy / n) / n;
111        }
112
113        /**
114         * Get the value of u20
115         * @return the value of u20
116         */
117        public double u20() {
118                return Mx2 / n;
119        }
120
121        /**
122         * Get the value of u02
123         * @return the value of u02
124         */
125        public double u02() {
126                return My2 / n;
127        }
128
129        /**
130         * Get the value of m10
131         * @return the value of m10
132         */
133        public double m10() {
134                return mx;
135        }
136
137        /**
138         * Get the value of m01
139         * @return the value of m01
140         */
141        public double m01() {
142                return my;
143        }
144
145        @Override
146        public MomentFeature clone() {
147                try {
148                        return (MomentFeature) super.clone();
149                } catch (CloneNotSupportedException e) {
150                        throw new AssertionError(e);
151                }
152        }
153
154        /**
155         * Create an ellipse based on the features parameters
156         * @return an ellipse
157         */
158        public Ellipse getEllipse() {
159                return getEllipse(1);
160        }
161        /**
162         * Create an ellipse based on the features parameters. Scale
163         * the ellipse about its centre with the given scale factor.
164         * @param sf the scale factor
165         * @return an ellipse
166         */
167        public Ellipse getEllipse(float sf) {
168                float u = (float) m10();
169                float v = (float) m01();
170                Matrix sm = new Matrix(new double[][]{
171                                {u20(),u11()},
172                                {u11(),u02()}
173                });
174                return EllipseUtilities.ellipseFromCovariance(u, v, sm, sf);
175        }
176        
177        /**
178         * Create an ellipse based on the features parameters. Scale
179         * the ellipse about its centre with the given scale factor.
180         * @param sf the scale factor
181         * @return an ellipse
182         */
183        public Circle getCircle(float sf) {
184                Ellipse e = getEllipse(sf);
185                Point2d p = e.calculateCentroid();
186                return new Circle(p.getX(),p.getY(),(float)(e.getMajor() + e.getMinor())/2);
187        }
188
189        /**
190         * Create a rotated rectangle that fits around an ellipse
191         * based on the features parameters. Scale the ellipse about 
192         * its centre with the given scale factor.
193         * @param sf the scale factor
194         * @return a polygon representing a rotated rectangle
195         */
196        public Polygon getEllipseBoundingBox(float sf) {
197//              double xx = u20();
198//              double xy = u11();
199//              double yy = u02();
200//
201//              double theta = 0.5 * Math.atan2(2*xy, xx-yy);
202//
203//              double trace = xx + yy;
204//              double det = (xx*yy) - (xy*xy);
205//              double [] eigval = QuadraticEquation.solveGeneralQuadratic(1, -1*trace, det);
206//
207//              double a = Math.sqrt(eigval[1]) * sf * 2;
208//              double b = Math.sqrt(eigval[0]) * sf * 2;
209//
210//              float Ax = (float) (a*Math.cos(theta) + b*Math.cos((Math.PI/2) + theta) + mx);
211//              float Ay = (float) (a*Math.sin(theta) + b*Math.sin((Math.PI/2) + theta) + my);
212//              float Bx = (float) (a*Math.cos(theta) - b*Math.cos((Math.PI/2) + theta) + mx);
213//              float By = (float) (a*Math.sin(theta) - b*Math.sin((Math.PI/2) + theta) + my);
214//              float Cx = (float) (-a*Math.cos(theta) - b*Math.cos((Math.PI/2) + theta) + mx);
215//              float Cy = (float) (-a*Math.sin(theta) - b*Math.sin((Math.PI/2) + theta) + my);
216//              float Dx = (float) (-a*Math.cos(theta) + b*Math.cos((Math.PI/2) + theta) + mx);
217//              float Dy = (float) (-a*Math.sin(theta) + b*Math.sin((Math.PI/2) + theta) + my);
218//
219//              return new Polygon(new Point2dImpl(Ax,Ay), new Point2dImpl(Bx,By), new Point2dImpl(Cx,Cy), new Point2dImpl(Dx,Dy));
220                return this.getEllipse(sf).calculateOrientedBoundingBox();
221        }
222
223        /**
224         * Get the primary orientation of the feature.
225         * @return the orientation
226         */
227        public double getOrientation() {
228                double xx = u20();
229                double xy = u11();
230                double yy = u02();
231
232                return 0.5 * Math.atan2(2*xy, xx-yy);
233        }
234        
235        protected double[] getEllipseBoundingRectsData(double sf) {
236                double xx = u20();
237                double xy = u11();
238                double yy = u02();
239
240                double theta = 0.5 * Math.atan2(2*xy, xx-yy);
241
242                double trace = xx + yy;
243                double det = (xx*yy) - (xy*xy);
244                double [] eigval = QuadraticEquation.solveGeneralQuadratic(1, -1*trace, det);
245
246                double a = 1.0 + Math.sqrt(Math.abs(eigval[1])) * sf * 4;
247                double b = 1.0 + Math.sqrt(Math.abs(eigval[0])) * sf * 4;
248
249                double [] data = {Math.max(a, b), Math.min(a, b), theta};
250                
251                return data;
252        }
253        
254        /**
255         * Extract a rectangular image patch centered on the feature
256         * with the same primary orientation and a given scale factor. 
257         * @param image the image to extract from
258         * @param sf the scale factor
259         * @return a rectangular image patch
260         */
261        public FImage extractEllipsePatch(FImage image, double sf) {
262                double [] data = getEllipseBoundingRectsData(sf);
263                double height = data[1], width = data[0], ori = data[2];                
264                
265                int sx = (int) Math.rint(width);
266                int sy = (int) Math.rint(height);
267
268                FImage patch = new FImage(sx, sy);
269
270                //extract pixels
271                for (int y=0; y<sy; y++) {
272                        for (int x=0; x<sx; x++) {
273                                double xbar = x - sx / 2.0;
274                                double ybar = y - sy / 2.0;
275
276                                double xx = (xbar * Math.cos(ori) - ybar * Math.sin(ori)) + mx;
277                                double yy = (xbar * Math.sin(ori) + ybar * Math.cos(ori)) + my;
278
279                                patch.setPixel(x, y, image.getPixelInterp(xx, yy));
280                        }
281                }
282
283                return patch;
284        }
285}