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.Collection;
034import java.util.Collections;
035import java.util.Iterator;
036import java.util.List;
037
038import org.openimaj.citation.annotation.Reference;
039import org.openimaj.citation.annotation.ReferenceType;
040import org.openimaj.math.geometry.line.Line2d;
041import org.openimaj.math.geometry.point.Point2d;
042import org.openimaj.math.geometry.point.Point2dImpl;
043import org.openimaj.math.geometry.shape.util.PolygonUtils;
044import org.openimaj.math.geometry.shape.util.RotatingCalipers;
045
046import Jama.Matrix;
047
048/**
049 * A polygon, modelled as a list of vertices. Polygon extends {@link PointList},
050 * so the vertices are the underlying {@link PointList#points}, and they are
051 * considered to be joined in order.
052 * 
053 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
054 */
055public class Polygon extends PointList implements Shape
056{
057        /**
058         * Polygons can contain other polygons. If the polygon is representing a
059         * shape, then the inner polygons can represent holes in the polygon or
060         * other polygons within the polygon.
061         */
062        private List<Polygon> innerPolygons = new ArrayList<Polygon>();
063
064        /** If this polygon is a hole within another polygon, this is set to true */
065        private boolean isHole = false;
066
067        /**
068         * Constructs an empty polygon to which vertices may be added.
069         */
070        public Polygon()
071        {
072                this(false);
073        }
074
075        /**
076         * Constructs an empty polygon to which vertices may be added. The boolean
077         * parameter determines whether this polygon will represent a hole (rather
078         * than a solid).
079         * 
080         * @param representsHole
081         *            Whether the polygon will represent a hole.
082         */
083        public Polygon(boolean representsHole)
084        {
085                this.isHole = representsHole;
086        }
087
088        /**
089         * Construct a Polygon from vertices
090         * 
091         * @param vertices
092         *            the vertices
093         */
094        public Polygon(Point2d... vertices) {
095                super(vertices);
096        }
097
098        /**
099         * Construct a Polygon from vertices
100         * 
101         * @param vertices
102         *            the vertices
103         */
104        public Polygon(Collection<? extends Point2d> vertices) {
105                this(vertices, false);
106        }
107
108        /**
109         * Construct a Polygon from the vertices, possibly copying the vertices
110         * first
111         * 
112         * @param vertices
113         *            the vertices
114         * @param copy
115         *            should the vertices be copied
116         */
117        public Polygon(Collection<? extends Point2d> vertices, boolean copy) {
118                super(vertices, copy);
119        }
120
121        /**
122         * Get the vertices of the polygon
123         * 
124         * @return the vertices
125         */
126        public List<Point2d> getVertices() {
127                return points;
128        }
129
130        /**
131         * Get the number of vertices
132         * 
133         * @return the number of vertices
134         */
135        public int nVertices() {
136                if (isClosed())
137                        return points.size() - 1;
138                return points.size();
139        }
140
141        /**
142         * Is the polygon closed (i.e. is the last vertex equal to the first)?
143         * 
144         * @return true if closed; false if open
145         */
146        public boolean isClosed() {
147                if (points.size() > 0 && points.get(0).getX() == points.get(points.size() - 1).getX()
148                                && points.get(0).getY() == points.get(points.size() - 1).getY())
149                        return true;
150                return false;
151        }
152
153        /**
154         * Close the polygon if it's not already closed
155         */
156        public void close() {
157                if (!isClosed() && points.size() > 0)
158                        points.add(points.get(0));
159        }
160
161        /**
162         * Open the polygon if it's closed
163         */
164        public void open() {
165                if (isClosed() && points.size() > 1)
166                        points.remove(points.size() - 1);
167        }
168
169        /**
170         * Test whether the point p is inside the polygon using the winding rule
171         * algorithm. Also tests whether the point is in any of the inner polygons.
172         * If the inner polygon represents a hole and the point is within that
173         * polygon then the point is not within this polygon.
174         * 
175         * @param point
176         *            the point to test
177         * @return true if the point is inside; false otherwise
178         */
179        @Override
180        public boolean isInside(Point2d point) {
181                final boolean isClosed = isClosed();
182                if (!isClosed)
183                        close();
184
185                boolean isOdd = false;
186
187                for (int pp = 0; pp < getNumInnerPoly(); pp++)
188                {
189                        final List<Point2d> v = getInnerPoly(pp).getVertices();
190                        int j = v.size() - 1;
191                        for (int i = 0; i < v.size(); i++) {
192                                if (v.get(i).getY() < point.getY() && v.get(j).getY() >= point.getY() ||
193                                                v.get(j).getY() < point.getY() && v.get(i).getY() >= point.getY())
194                                {
195                                        if (v.get(i).getX() + (point.getY() - v.get(i).getY()) /
196                                                        (v.get(j).getY() - v.get(i).getY()) * (v.get(j).getX() - v.get(i).getX()) < point.getX())
197                                        {
198                                                isOdd = !isOdd;
199                                        }
200                                }
201                                j = i;
202                        }
203                }
204
205                if (!isClosed)
206                        open();
207                return isOdd;
208        }
209
210        /**
211         * {@inheritDoc}
212         */
213        @Override
214        public Polygon clone() {
215                final Polygon clone = new Polygon();
216                clone.setIsHole(isHole);
217
218                for (final Point2d p : points)
219                        clone.points.add(p.copy());
220
221                for (final Polygon innerPoly : innerPolygons)
222                        clone.addInnerPolygon(innerPoly.clone());
223
224                return clone;
225        }
226
227        /**
228         * Calculates the difference between two polygons and returns a new polygon.
229         * It assumes that the given polygon and this polygon have the same number
230         * of vertices.
231         * 
232         * @param p
233         *            the polygon to subtract.
234         * @return the difference polygon
235         */
236        public Polygon difference(Polygon p)
237        {
238                final List<Point2d> v = new ArrayList<Point2d>();
239
240                for (int i = 0; i < points.size(); i++)
241                        v.add(new Point2dImpl(
242                                        points.get(i).getX() - p.getVertices().get(i).getX(),
243                                        points.get(i).getY() - p.getVertices().get(i).getY()));
244
245                final Polygon p2 = new Polygon(v);
246                for (int i = 0; i < innerPolygons.size(); i++)
247                        p2.addInnerPolygon(innerPolygons.get(i).difference(
248                                        p2.getInnerPoly(i + 1)));
249
250                return p2;
251        }
252
253        @Override
254        public double calculateArea() {
255                return Math.abs(calculateSignedArea());
256        }
257
258        /**
259         * Calculate the area of the polygon. This does not take into account holes
260         * in the polygon.
261         * 
262         * @return the area of the polygon
263         */
264        public double calculateSignedArea() {
265                final boolean closed = isClosed();
266                double area = 0;
267
268                if (!closed)
269                        close();
270
271                // TODO: This does not take into account the winding
272                // rule and therefore holes
273                for (int k = 0; k < points.size() - 1; k++) {
274                        final float ik = points.get(k).getX();
275                        final float jk = points.get(k).getY();
276                        final float ik1 = points.get(k + 1).getX();
277                        final float jk1 = points.get(k + 1).getY();
278
279                        area += ik * jk1 - ik1 * jk;
280                }
281
282                if (!closed)
283                        open();
284
285                return 0.5 * area;
286        }
287
288        /**
289         * Calls {@link Polygon#intersectionArea(Shape, int)} with 1 step per pixel
290         * dimension. Subsequently this function returns the shared whole pixels of
291         * this polygon and that.
292         * 
293         * @param that
294         * @return intersection area
295         */
296        @Override
297        public double intersectionArea(Shape that) {
298                return this.intersectionArea(that, 1);
299        }
300
301        /**
302         * Return an estimate for the area of the intersection of this polygon and
303         * another polygon. For each pixel step 1 is added if the point is inside
304         * both polygons. For each pixel, perPixelPerDimension steps are taken.
305         * Subsequently the intersection is:
306         * 
307         * sumIntersections / (perPixelPerDimension * perPixelPerDimension)
308         * 
309         * @param that
310         * @return normalised intersection area
311         */
312        @Override
313        public double intersectionArea(Shape that, int nStepsPerDimension) {
314                final Rectangle overlapping = this.calculateRegularBoundingBox().overlapping(that.calculateRegularBoundingBox());
315                if (overlapping == null)
316                        return 0;
317                double intersection = 0;
318                final double step = Math.max(overlapping.width, overlapping.height) / (double) nStepsPerDimension;
319                double nReads = 0;
320                for (float x = overlapping.x; x < overlapping.x + overlapping.width; x += step) {
321                        for (float y = overlapping.y; y < overlapping.y + overlapping.height; y += step) {
322                                final boolean insideThis = this.isInside(new Point2dImpl(x, y));
323                                final boolean insideThat = that.isInside(new Point2dImpl(x, y));
324                                nReads++;
325                                if (insideThis && insideThat) {
326                                        intersection++;
327                                }
328                        }
329                }
330
331                return (intersection / nReads) * (overlapping.width * overlapping.height);
332        }
333
334        /**
335         * {@inheritDoc}
336         * 
337         * @see org.openimaj.math.geometry.shape.Shape#asPolygon()
338         */
339        @Override
340        public Polygon asPolygon() {
341                return this;
342        }
343
344        /**
345         * Add a vertex to the polygon
346         * 
347         * @param x
348         *            x-coordinate of the vertex
349         * @param y
350         *            y-coordinate of the vertex
351         */
352        public void addVertex(float x, float y) {
353                points.add(new Point2dImpl(x, y));
354        }
355
356        /**
357         * Add a vertex to the polygon
358         * 
359         * @param pt
360         *            coordinate of the vertex
361         */
362        public void addVertex(Point2d pt) {
363                points.add(pt);
364        }
365
366        /**
367         * Iterates through the vertices and rounds all vertices to round integers.
368         * Side-affects this polygon.
369         * 
370         * @return this polygon
371         */
372        public Polygon roundVertices()
373        {
374                final Iterator<Point2d> i = this.iterator();
375                while (i.hasNext())
376                {
377                        final Point2d p = i.next();
378                        final Point2dImpl p2 = new Point2dImpl((int) p.getX(), (int) p.getY());
379
380                        int xx = -1;
381                        if ((xx = this.points.indexOf(p2)) != -1 &&
382                                        this.points.get(xx) != p)
383                                i.remove();
384                        else
385                        {
386                                p.setX(p2.x);
387                                p.setY(p2.y);
388                        }
389                }
390
391                for (final Polygon pp : innerPolygons)
392                        pp.roundVertices();
393
394                return this;
395        }
396
397        /**
398         * Return whether this polygon has no vertices or not.
399         * 
400         * @return TRUE if this polygon has no vertices
401         */
402        public boolean isEmpty()
403        {
404                return this.points.isEmpty() && innerPolygons.isEmpty();
405        }
406
407        /**
408         * Returns the number of inner polygons in this polygon including this
409         * polygon.
410         * 
411         * @return the number of inner polygons in this polygon.
412         */
413        public int getNumInnerPoly()
414        {
415                return innerPolygons.size() + 1;
416        }
417
418        /**
419         * Get the inner polygon at the given index. Note that index 0 will return
420         * this polygon, while index i+1 will return the inner polygon i.
421         * 
422         * @param index
423         *            the index of the polygon to get
424         * @return The inner polygon at the given index.
425         */
426        public Polygon getInnerPoly(int index)
427        {
428                if (index == 0)
429                        return this;
430                return innerPolygons.get(index - 1);
431        }
432
433        /**
434         * Add an inner polygon to this polygon. If there is no main polygon defined
435         * (the number of vertices is zero) then the given inner polygon will become
436         * the main polygon if the <code>inferOuter</code> boolean is true. If this
437         * variable is false, the inner polygon will be added to the list of inner
438         * polygons for this polygon regardless of whether a main polygon exists.
439         * When the main polygon is inferred from the given polygon, the vertices
440         * are copied into this polygon's vertex list.
441         * 
442         * @param p
443         *            The inner polygon to add
444         * @param inferOuter
445         *            Whether to infer the outer polygon from this inner one
446         */
447        public void addInnerPolygon(Polygon p, boolean inferOuter)
448        {
449                if (!inferOuter)
450                {
451                        this.innerPolygons.add(p);
452                }
453                else
454                {
455                        if (this.points.size() == 0)
456                        {
457                                this.points.addAll(p.points);
458                                this.isHole = p.isHole;
459                        }
460                        else
461                        {
462                                this.addInnerPolygon(p, false);
463                        }
464                }
465        }
466
467        /**
468         * Add an inner polygon to this polygon. If there is no main polygon defined
469         * (the number of vertices is zero) then the given inner polygon will become
470         * the main polygon.
471         * 
472         * @param p
473         *            The inner polygon to add
474         */
475        public void addInnerPolygon(Polygon p)
476        {
477                this.addInnerPolygon(p, true);
478        }
479
480        /**
481         * Returns the list of inner polygons.
482         * 
483         * @return the list of inner polygons
484         */
485        public List<Polygon> getInnerPolys()
486        {
487                return this.innerPolygons;
488        }
489
490        /**
491         * Set whether this polygon represents a hole in another polygon.
492         * 
493         * @param isHole
494         *            Whether this polygon is a whole.
495         */
496        public void setIsHole(boolean isHole)
497        {
498                this.isHole = isHole;
499        }
500
501        /**
502         * Returns whether this polygon is representing a hole in another polygon.
503         * 
504         * @return Whether this polygon is representing a hole in another polygon.
505         */
506        public boolean isHole()
507        {
508                return this.isHole;
509        }
510
511        /**
512         * {@inheritDoc}
513         * 
514         * @see java.lang.Object#equals(java.lang.Object)
515         */
516        @Override
517        public boolean equals(Object obj)
518        {
519                return (obj instanceof Polygon) &&
520                                this.equals((Polygon) obj);
521        }
522
523        /**
524         * Specific equals method for polygons where the polgyons are tested for the
525         * vertices being in the same order. If the vertices are not in the vertex
526         * list in the same manner but are in the same order (when wrapped around)
527         * the method will return true. So the triangles below will return true:
528         * 
529         * {[1,1],[2,2],[1,2]} and {[1,2],[1,1],[2,2]}
530         * 
531         * @param p
532         *            The polygon to test against
533         * @return TRUE if the polygons are the same.
534         */
535        public boolean equals(Polygon p)
536        {
537                if (isHole() != p.isHole())
538                        return false;
539                if (this.points.size() != p.points.size())
540                        return false;
541                if (this.isEmpty() && p.isEmpty())
542                        return true;
543
544                final int i = this.points.indexOf(p.points.get(0));
545                if (i == -1)
546                        return false;
547
548                final int s = this.points.size();
549                for (int n = 0; n < s; n++)
550                {
551                        if (!p.points.get(n).equals(this.points.get((n + i) % s)))
552                                return false;
553                }
554
555                return true;
556        }
557
558        /**
559         * {@inheritDoc}
560         * 
561         * @see java.lang.Object#hashCode()
562         */
563        @Override
564        public int hashCode()
565        {
566                return points.hashCode() * (isHole() ? -1 : 1);
567        }
568
569        /**
570         * Displays the complete list of vertices unless the number of vertices is
571         * greater than 10 - then a sublist is shown of 5 from the start and 5 from
572         * the end separated by ellipses.
573         * 
574         * {@inheritDoc}
575         * 
576         * @see java.lang.Object#toString()
577         */
578        @Override
579        public String toString()
580        {
581                final StringBuffer sb = new StringBuffer();
582                if (isHole())
583                        sb.append("H");
584                final int len = 10;
585                if (points.size() < len)
586                        sb.append(points.toString());
587                else
588                        sb.append(points.subList(0, len / 2).toString() + "..." +
589                                        points.subList(points.size() - len / 2, points.size())
590                                                        .toString());
591
592                if (innerPolygons.size() > 0)
593                {
594                        sb.append("\n    - " + innerPolygons.size() + " inner polygons:");
595                        for (final Polygon ip : innerPolygons)
596                                sb.append("\n       + " + ip.toString());
597                }
598
599                return sb.toString();
600        }
601
602        /**
603         * Returns the intersection of this polygon and the given polygon.
604         * 
605         * @param p2
606         *            The polygon to intersect with.
607         * @return The resulting polygon intersection
608         */
609        public Polygon intersect(Polygon p2)
610        {
611                return new PolygonUtils().intersection(this, p2);
612        }
613
614        /**
615         * Returns the union of this polygon and the given polygon.
616         * 
617         * @param p2
618         *            The polygon to union with.
619         * @return The resulting polygon union
620         */
621        public Polygon union(Polygon p2)
622        {
623                return new PolygonUtils().union(this, p2);
624        }
625
626        /**
627         * Returns the XOR of this polygon and the given polygon.
628         * 
629         * @param p2
630         *            The polygon to XOR with.
631         * @return The resulting polygon
632         */
633        public Polygon xor(Polygon p2)
634        {
635                return new PolygonUtils().xor(this, p2);
636        }
637
638        /**
639         * Reduce the number of vertices in a polygon. This implementation is based
640         * on a modified Ramer-Douglas–Peucker algorithm that is designed to work
641         * with polygons rather than polylines. The implementation searches for two
642         * initialisation points that are approximatatly maximally distant, then
643         * performs the polyline algorithm on the two segments, and finally ensures
644         * that the joins in the segments are consistent with the given epsilon
645         * value.
646         * 
647         * @param eps
648         *            distance value below which a vertex can be ignored
649         * @return new polygon that approximates this polygon, but with fewer
650         *         vertices
651         */
652        public Polygon reduceVertices(double eps)
653        {
654                if (eps == 0 || nVertices() <= 3)
655                        return this.clone();
656
657                int prevStartIndex = 0;
658                int startIndex = 0;
659                for (int init = 0; init < 3; init++) {
660                        double dmax = 0;
661                        prevStartIndex = startIndex;
662
663                        final Point2d first = points.get(startIndex);
664                        for (int i = 0; i < points.size(); i++) {
665                                final double d;
666                                d = Line2d.distance(first, points.get(i));
667
668                                if (d > dmax) {
669                                        startIndex = i;
670                                        dmax = d;
671                                }
672                        }
673                }
674
675                if (prevStartIndex > startIndex) {
676                        final int tmp = prevStartIndex;
677                        prevStartIndex = startIndex;
678                        startIndex = tmp;
679                }
680
681                final List<Point2d> l1 = new ArrayList<Point2d>();
682                l1.addAll(points.subList(prevStartIndex, startIndex + 1));
683
684                final List<Point2d> l2 = new ArrayList<Point2d>();
685                l2.addAll(points.subList(startIndex, points.size()));
686                l2.addAll(points.subList(0, prevStartIndex + 1));
687
688                final Polygon newPoly = new Polygon();
689                final List<Point2d> line1 = ramerDouglasPeucker(l1, eps);
690                final List<Point2d> line2 = ramerDouglasPeucker(l2, eps);
691                newPoly.points.addAll(line1.subList(0, line1.size() - 1));
692                newPoly.points.addAll(line2.subList(0, line2.size() - 1));
693
694                // deal with the joins
695                // if (newPoly.nVertices() > 3) {
696                // Point2d r1 = null, r2 = null;
697                // Point2d p0 = newPoly.points.get(newPoly.points.size() - 1);
698                // Point2d p1 = newPoly.points.get(0);
699                // Point2d p2 = newPoly.points.get(1);
700                //
701                // Line2d l = new Line2d(p0, p2);
702                // Line2d norm = l.getNormal(p1);
703                // Point2d intersect = l.getIntersection(norm).intersectionPoint;
704                // if (intersect != null && Line2d.distance(p1, intersect) <= eps) {
705                // // remove p1
706                // r1 = p1;
707                // }
708                //
709                // p0 = newPoly.points.get(line1.size() - 1);
710                // p1 = newPoly.points.get(line1.size());
711                // p2 = newPoly.points.get((line1.size() + 1) >= newPoly.size() ? 0 :
712                // (line1.size() + 1));
713                //
714                // l = new Line2d(p0, p2);
715                // norm = l.getNormal(p1);
716                // intersect = l.getIntersection(norm).intersectionPoint;
717                // if (intersect != null && Line2d.distance(p1, intersect) <= eps) {
718                // // remove p2
719                // r2 = p2;
720                // }
721                //
722                // if (r1 != null) {
723                // newPoly.points.remove(r1);
724                // }
725                // if (newPoly.nVertices() > 3 && r2 != null) {
726                // newPoly.points.remove(r2);
727                // }
728                // }
729
730                if (!newPoly.isClockwise()) {
731                        Collections.reverse(newPoly.points);
732                }
733
734                for (final Polygon ppp : innerPolygons)
735                        newPoly.addInnerPolygon(ppp.reduceVertices(eps));
736
737                return newPoly;
738        }
739
740        /**
741         * Ramer Douglas Peucker polyline algorithm
742         * 
743         * @param p
744         *            the polyline to simplify
745         * @param eps
746         *            distance value below which a vertex can be ignored
747         * @return the simplified polyline
748         */
749        private static List<Point2d> ramerDouglasPeucker(List<Point2d> p, double eps) {
750                // Find the point with the maximum distance
751                double dmax = 0;
752                int index = 0;
753                final int end = p.size() - 1;
754                final Line2d l = new Line2d(p.get(0), p.get(end));
755                for (int i = 1; i < end - 1; i++) {
756                        final double d;
757
758                        final Point2d p1 = p.get(i);
759                        final Line2d norm = l.getNormal(p1);
760                        final Point2d p2 = l.getIntersection(norm).intersectionPoint;
761                        if (p2 == null)
762                                continue;
763                        d = Line2d.distance(p1, p2);
764
765                        if (d > dmax) {
766                                index = i;
767                                dmax = d;
768                        }
769                }
770
771                final List<Point2d> newPoly = new ArrayList<Point2d>();
772
773                // If max distance is greater than epsilon, recursively simplify
774                if (dmax > eps) {
775                        // Recursively call RDP
776                        final List<Point2d> line1 = ramerDouglasPeucker(p.subList(0, index + 1), eps);
777                        final List<Point2d> line2 = ramerDouglasPeucker(p.subList(index, end + 1), eps);
778                        newPoly.addAll(line1.subList(0, line1.size() - 1));
779                        newPoly.addAll(line2);
780                } else {
781                        newPoly.add(p.get(0));
782                        newPoly.add(p.get(end));
783                }
784
785                // Return the result
786                return newPoly;
787        }
788
789        /**
790         * Apply a 3x3 transform matrix to a copy of the polygon and return it
791         * 
792         * @param transform
793         *            3x3 transform matrix
794         * @return the transformed polygon
795         */
796        @Override
797        public Polygon transform(Matrix transform) {
798                final List<Point2d> newVertices = new ArrayList<Point2d>();
799
800                for (final Point2d p : points) {
801                        final Matrix p1 = new Matrix(3, 1);
802                        p1.set(0, 0, p.getX());
803                        p1.set(1, 0, p.getY());
804                        p1.set(2, 0, 1);
805
806                        final Matrix p2_est = transform.times(p1);
807
808                        final Point2d out = new Point2dImpl((float) (p2_est.get(0, 0) / p2_est.get(2, 0)),
809                                        (float) (p2_est.get(1, 0) / p2_est.get(2, 0)));
810
811                        newVertices.add(out);
812                }
813
814                final Polygon p = new Polygon(newVertices);
815                for (final Polygon pp : innerPolygons)
816                        p.addInnerPolygon(pp.transform(transform));
817                return p;
818        }
819
820        /**
821         * Compute the regular (oriented to the axes) bounding box of the polygon.
822         * 
823         * @return the regular bounding box as [x,y,width,height]
824         */
825        @Override
826        public Rectangle calculateRegularBoundingBox() {
827                float xmin = Float.MAX_VALUE, xmax = -Float.MAX_VALUE, ymin = Float.MAX_VALUE, ymax = -Float.MAX_VALUE;
828
829                for (int pp = 0; pp < getNumInnerPoly(); pp++)
830                {
831                        final Polygon ppp = getInnerPoly(pp);
832                        for (final Point2d p : ppp.getVertices()) {
833                                final float px = p.getX();
834                                final float py = p.getY();
835                                if (px < xmin)
836                                        xmin = px;
837                                if (px > xmax)
838                                        xmax = px;
839                                if (py < ymin)
840                                        ymin = py;
841                                if (py > ymax)
842                                        ymax = py;
843
844                        }
845                }
846
847                return new Rectangle(xmin, ymin, xmax - xmin, ymax - ymin);
848        }
849
850        /**
851         * Translate the polygons position
852         * 
853         * @param x
854         *            x-translation
855         * @param y
856         *            y-translation
857         */
858        @Override
859        public void translate(float x, float y) {
860                for (int pp = 0; pp < getNumInnerPoly(); pp++)
861                {
862                        final Polygon ppp = getInnerPoly(pp);
863                        for (final Point2d p : ppp.getVertices()) {
864                                p.setX(p.getX() + x);
865                                p.setY(p.getY() + y);
866                        }
867                }
868        }
869
870        /**
871         * Scale the polygon by the given amount about (0,0). Scalefactors between 0
872         * and 1 shrink the polygon.
873         * 
874         * @param sc
875         *            the scale factor.
876         */
877        @Override
878        public void scale(float sc) {
879                for (int pp = 0; pp < getNumInnerPoly(); pp++)
880                {
881                        final Polygon ppp = getInnerPoly(pp);
882                        for (final Point2d p : ppp.getVertices()) {
883                                p.setX(p.getX() * sc);
884                                p.setY(p.getY() * sc);
885                        }
886                }
887        }
888
889        /**
890         * Scale the polygon only in the x-direction by the given amount about
891         * (0,0). Scale factors between 0 and 1 will shrink the polygon
892         * 
893         * @param sc
894         *            The scale factor
895         * @return this polygon
896         */
897        @Override
898        public Polygon scaleX(float sc)
899        {
900                for (int pp = 0; pp < getNumInnerPoly(); pp++)
901                {
902                        final Polygon ppp = getInnerPoly(pp);
903                        for (final Point2d p : ppp.getVertices()) {
904                                p.setX(p.getX() * sc);
905                        }
906                }
907                return this;
908        }
909
910        /**
911         * Scale the polygon only in the y-direction by the given amount about
912         * (0,0). Scale factors between 0 and 1 will shrink the polygon
913         * 
914         * @param sc
915         *            The scale factor
916         * @return this polygon
917         */
918        @Override
919        public Polygon scaleY(float sc)
920        {
921                for (int pp = 0; pp < getNumInnerPoly(); pp++)
922                {
923                        final Polygon ppp = getInnerPoly(pp);
924                        for (final Point2d p : ppp.getVertices()) {
925                                p.setY(p.getY() * sc);
926                        }
927                }
928                return this;
929        }
930
931        /**
932         * Scale the polygon by the given amount about (0,0). Scale factors between
933         * 0 and 1 shrink the polygon.
934         * 
935         * @param scx
936         *            the scale factor in the x direction
937         * @param scy
938         *            the scale factor in the y direction.
939         * @return this polygon
940         */
941        @Override
942        public Polygon scaleXY(float scx, float scy)
943        {
944                for (int pp = 0; pp < getNumInnerPoly(); pp++)
945                {
946                        final Polygon ppp = getInnerPoly(pp);
947                        for (final Point2d p : ppp.getVertices()) {
948                                p.setX(p.getX() * scx);
949                                p.setY(p.getY() * scy);
950                        }
951                }
952                return this;
953        }
954
955        /**
956         * Scale the polygon by the given amount about the given point. Scalefactors
957         * between 0 and 1 shrink the polygon.
958         * 
959         * @param centre
960         *            the centre of the scaling operation
961         * @param sc
962         *            the scale factor
963         */
964        @Override
965        public void scale(Point2d centre, float sc) {
966                this.translate(-centre.getX(), -centre.getY());
967                for (int pp = 0; pp < getNumInnerPoly(); pp++)
968                {
969                        final Polygon ppp = getInnerPoly(pp);
970                        for (final Point2d p : ppp.getVertices()) {
971                                p.setX(p.getX() * sc);
972                                p.setY(p.getY() * sc);
973                        }
974                }
975                this.translate(centre.getX(), centre.getY());
976        }
977
978        /**
979         * Calculate the centroid of the polygon.
980         * 
981         * @return calls {@link #calculateFirstMoment()};
982         */
983        @Override
984        public Point2d calculateCentroid() {
985                final double[] pt = calculateFirstMoment();
986                return new Point2dImpl((float) pt[0], (float) pt[1]);
987        }
988
989        /**
990         * Treating the polygon as a continuous piecewise function, calculate
991         * exactly the first moment. This follows working presented by Carsten
992         * Steger in "On the Calculation of Moments of Polygons" ,
993         * 
994         * @return the first moment
995         */
996        @Reference(
997                        author = { "Carsten Steger" },
998                        title = "On the Calculation of Moments of Polygons",
999                        type = ReferenceType.Techreport,
1000                        month = "August",
1001                        year = "1996",
1002                        url = "http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.29.8765&rep=rep1&type=pdf"
1003                        )
1004                        public double[] calculateFirstMoment() {
1005                final boolean closed = isClosed();
1006
1007                if (!closed)
1008                        close();
1009
1010                double area = 0;
1011                double ax = 0;
1012                double ay = 0;
1013                // TODO: This does not take into account the winding
1014                // rule and therefore holes
1015                for (int k = 0; k < points.size() - 1; k++) {
1016                        final float xk = points.get(k).getX();
1017                        final float yk = points.get(k).getY();
1018                        final float xk1 = points.get(k + 1).getX();
1019                        final float yk1 = points.get(k + 1).getY();
1020
1021                        final float shared = xk * yk1 - xk1 * yk;
1022                        area += shared;
1023                        ax += (xk + xk1) * shared;
1024                        ay += (yk + yk1) * shared;
1025                }
1026
1027                if (!closed)
1028                        open();
1029
1030                area *= 0.5;
1031
1032                return new double[] { ax / (6 * area), ay / (6 * area) };
1033        }
1034
1035        /**
1036         * Treating the polygon as a continuous piecewise function, calculate
1037         * exactly the second moment. This follows working presented by Carsten
1038         * Steger in "On the Calculation of Moments of Polygons" ,
1039         * 
1040         * @return the second moment as an array with values: (axx,axy,ayy)
1041         */
1042        @Reference(
1043                        author = { "Carsten Steger" },
1044                        title = "On the Calculation of Moments of Polygons",
1045                        type = ReferenceType.Techreport,
1046                        month = "August",
1047                        year = "1996",
1048                        url = "http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.29.8765&rep=rep1&type=pdf"
1049                        )
1050                        public double[] calculateSecondMoment() {
1051                final boolean closed = isClosed();
1052                final double area = calculateSignedArea();
1053
1054                if (!closed)
1055                        close();
1056
1057                double axx = 0;
1058                double ayy = 0;
1059                double axy = 0;
1060                // TODO: This does not take into account the winding
1061                // rule and therefore holes
1062                for (int k = 0; k < points.size() - 1; k++) {
1063                        final float xk = points.get(k).getX();
1064                        final float yk = points.get(k).getY();
1065                        final float xk1 = points.get(k + 1).getX();
1066                        final float yk1 = points.get(k + 1).getY();
1067
1068                        final float shared = xk * yk1 - xk1 * yk;
1069                        axx += (xk * xk + xk * xk1 + xk1 * xk1) * shared;
1070                        ayy += (yk * yk + yk * yk1 + yk1 * yk1) * shared;
1071                        axy += (2 * xk * yk + xk * yk1 + xk1 * yk + 2 * xk1 * yk1) * shared;
1072                }
1073
1074                if (!closed)
1075                        open();
1076
1077                return new double[] {
1078                                axx / (12 * area),
1079                                axy / (24 * area),
1080                                ayy / (12 * area)
1081                };
1082        }
1083
1084        /**
1085         * Treating the polygon as a continuous piecewise function, calculate
1086         * exactly the centralised second moment. This follows working presented by
1087         * Carsten Steger in "On the Calculation of Moments of Polygons" ,
1088         * 
1089         * @return the second moment as an array with values: (axx,axy,ayy)
1090         */
1091        @Reference(
1092                        author = { "Carsten Steger" },
1093                        title = "On the Calculation of Moments of Polygons",
1094                        type = ReferenceType.Techreport,
1095                        month = "August",
1096                        year = "1996",
1097                        url = "http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.29.8765&rep=rep1&type=pdf"
1098                        )
1099                        public double[] calculateSecondMomentCentralised() {
1100                final double[] firstMoment = this.calculateFirstMoment();
1101                final double[] secondMoment = this.calculateSecondMoment();
1102
1103                return new double[] {
1104                                secondMoment[0] - firstMoment[0] * firstMoment[0],
1105                                secondMoment[1] - firstMoment[0] * firstMoment[1],
1106                                secondMoment[2] - firstMoment[1] * firstMoment[1],
1107                };
1108
1109        }
1110
1111        /**
1112         * Calculates the principle direction of the connected component. This is
1113         * given by
1114         * <code>0.5 * atan( (M<sub>20</sub>-M<sub>02</sub>) / 2 * M<sub>11</sub> )</code>
1115         * so results in an angle between -PI and +PI.
1116         * 
1117         * @return The principle direction (-PI/2 to +PI/2 radians) of the connected
1118         *         component.
1119         */
1120        public double calculateDirection() {
1121                final double[] secondMoment = calculateSecondMomentCentralised();
1122                final double u20 = secondMoment[0];
1123                final double u02 = secondMoment[1];
1124                final double u11 = secondMoment[2];
1125                final double theta = 0.5 * Math.atan2((2 * u11), (u20 - u02));
1126
1127                return theta;
1128        }
1129
1130        /**
1131         * Using
1132         * {@link EllipseUtilities#ellipseFromCovariance(float, float, Matrix, float)}
1133         * and the {@link #calculateSecondMomentCentralised()} return the Ellipse
1134         * best fitting the shape of this polygon.
1135         * 
1136         * @return {@link Ellipse} defining the shape of the polygon
1137         */
1138        public Ellipse toEllipse() {
1139                final double[] secondMoment = calculateSecondMomentCentralised();
1140                final double u20 = secondMoment[0];
1141                final double u11 = secondMoment[1];
1142                final double u02 = secondMoment[2];
1143                final Point2d center = calculateCentroid();
1144                final Matrix sm = new Matrix(new double[][] {
1145                                new double[] { u20, u11 },
1146                                new double[] { u11, u02 },
1147                });
1148                // Used the sqrt(3) as the scale, not sure why. This is not correct.
1149                // Find the correct value!
1150                return EllipseUtilities.ellipseFromCovariance(
1151                                center.getX(),
1152                                center.getY(),
1153                                sm,
1154                                (float) Math.sqrt(3)
1155                                );
1156        }
1157
1158        /**
1159         * Test if the outer polygon is convex.
1160         * 
1161         * @return true if the outer polygon is convex; false if non-convex or less
1162         *         than two vertices
1163         */
1164        @Override
1165        public boolean isConvex() {
1166                final boolean isOriginallyClosed = this.isClosed();
1167                if (isOriginallyClosed)
1168                        this.open();
1169
1170                final int size = size();
1171
1172                if (size < 3)
1173                        return false;
1174
1175                float res = 0;
1176                for (int i = 0; i < size; i++) {
1177                        final Point2d p = points.get(i);
1178                        final Point2d tmp = points.get((i + 1) % size);
1179                        final Point2dImpl v = new Point2dImpl();
1180                        v.x = tmp.getX() - p.getX();
1181                        v.y = tmp.getY() - p.getY();
1182                        final Point2d u = points.get((i + 2) % size);
1183
1184                        if (i == 0) // in first loop direction is unknown, so save it in res
1185                                res = u.getX() * v.y - u.getY() * v.x + v.x * p.getY() - v.y * p.getX();
1186                        else
1187                        {
1188                                final float newres = u.getX() * v.y - u.getY() * v.x + v.x * p.getY() - v.y * p.getX();
1189                                if ((newres > 0 && res < 0) || (newres < 0 && res > 0))
1190                                        return false;
1191                        }
1192                }
1193
1194                if (isOriginallyClosed)
1195                        close();
1196
1197                return true;
1198        }
1199
1200        /**
1201         * Test if this has no inner polygons or sub-parts
1202         * 
1203         * @see Polygon#getInnerPolys()
1204         * 
1205         * @return true if this is polygon has no inner polygons; false otherwise.
1206         */
1207        public boolean hasNoInnerPolygons() {
1208                return innerPolygons == null || innerPolygons.size() == 0;
1209        }
1210
1211        @Override
1212        public double calculatePerimeter() {
1213                final Point2d p1 = points.get(0);
1214                float p1x = p1.getX();
1215                float p1y = p1.getY();
1216
1217                Point2d p2 = points.get(points.size() - 1);
1218                float p2x = p2.getX();
1219                float p2y = p2.getY();
1220
1221                double peri = Line2d.distance(p1x, p1y, p2x, p2y);
1222                for (int i = 1; i < this.points.size(); i++) {
1223                        p2 = points.get(i);
1224                        p2x = p2.getX();
1225                        p2y = p2.getY();
1226                        peri += Line2d.distance(p1x, p1y, p2x, p2y);
1227                        p1x = p2x;
1228                        p1y = p2y;
1229                }
1230
1231                return peri;
1232        }
1233
1234        /**
1235         * Test (outer) polygon path direction
1236         * 
1237         * @return true is the outer path direction is clockwise w.r.t OpenIMAJ
1238         *         coordinate system
1239         */
1240        public boolean isClockwise() {
1241                double signedArea = 0;
1242                for (int i = 0; i < this.points.size() - 1; i++) {
1243                        final float x1 = points.get(i).getX();
1244                        final float y1 = points.get(i).getY();
1245                        final float x2 = points.get(i + 1).getX();
1246                        final float y2 = points.get(i + 1).getY();
1247
1248                        signedArea += (x1 * y2 - x2 * y1);
1249                }
1250
1251                final float x1 = points.get(points.size() - 1).getX();
1252                final float y1 = points.get(points.size() - 1).getY();
1253                final float x2 = points.get(0).getX();
1254                final float y2 = points.get(0).getY();
1255                signedArea += (x1 * y2 - x2 * y1);
1256
1257                return signedArea >= 0;
1258        }
1259
1260        /**
1261         * Calculate convex hull using Melkman's algorithm. This is faster than the
1262         * {@link #calculateConvexHull()} method, but will only work for simple
1263         * polygons (those that don't self-intersect)
1264         * <p>
1265         * Based on http://softsurfer.com/Archive/algorithm_0203/algorithm_0203.htm,
1266         * but modified (fixed) to deal with the problem of the initialisation
1267         * points potentially being co-linear.
1268         * <p>
1269         * Copyright 2001, softSurfer (www.softsurfer.com) This code may be freely
1270         * used and modified for any purpose providing that this copyright notice is
1271         * included with it. SoftSurfer makes no warranty for this code, and cannot
1272         * be held liable for any real or imagined damage resulting from its use.
1273         * Users of this code must verify correctness for their application.
1274         * 
1275         * @return A polygon defining the shape of the convex hull
1276         */
1277        public Polygon calculateConvexHullMelkman() {
1278                if (this.points.size() <= 3)
1279                        return new Polygon(this.points);
1280
1281                final int n = this.points.size();
1282                int i = 1;
1283                while (i + 1 < n && isLeft(points.get(0), points.get(i), points.get(i + 1)) == 0)
1284                        i++;
1285
1286                if (n - i <= 3)
1287                        return new Polygon(this.points);
1288
1289                // initialize a deque D[] from bottom to top so that the
1290                // 1st three vertices of V[] are a counterclockwise triangle
1291                final Point2d[] D = new Point2d[2 * n + 1];
1292                int bot = n - 2, top = bot + 3; // initial bottom and top deque indices
1293                D[bot] = D[top] = points.get(i + 1); // 3rd vertex is at both bot and
1294                                                                                                // top
1295                if (isLeft(points.get(0), points.get(i), points.get(i + 1)) > 0) {
1296                        D[bot + 1] = points.get(0);
1297                        D[bot + 2] = points.get(i); // ccw vertices are: 2,0,1,2
1298                } else {
1299                        D[bot + 1] = points.get(i);
1300                        D[bot + 2] = points.get(0); // ccw vertices are: 2,1,0,2
1301                }
1302
1303                i += 2;
1304
1305                // compute the hull on the deque D[]
1306                for (; i < n; i++) { // process the rest of vertices
1307                        // test if next vertex is inside the deque hull
1308                        if ((isLeft(D[bot], D[bot + 1], points.get(i)) > 0) &&
1309                                        (isLeft(D[top - 1], D[top], points.get(i)) > 0))
1310                                continue; // skip an interior vertex
1311
1312                        // incrementally add an exterior vertex to the deque hull
1313                        // get the rightmost tangent at the deque bot
1314                        while (isLeft(D[bot], D[bot + 1], points.get(i)) <= 0)
1315                                ++bot; // remove bot of deque
1316                        D[--bot] = points.get(i); // insert V[i] at bot of deque
1317
1318                        // get the leftmost tangent at the deque top
1319                        while (isLeft(D[top - 1], D[top], points.get(i)) <= 0)
1320                                --top; // pop top of deque
1321                        D[++top] = points.get(i); // push V[i] onto top of deque
1322                }
1323
1324                // transcribe deque D[] to the output hull array H[]
1325                final Polygon H = new Polygon();
1326                final List<Point2d> vertices = H.getVertices();
1327                for (int h = 0; h <= (top - bot); h++)
1328                        vertices.add(D[bot + h]);
1329
1330                return H;
1331        }
1332
1333        private float isLeft(Point2d P0, Point2d P1, Point2d P2) {
1334                return (P1.getX() - P0.getX()) * (P2.getY() - P0.getY()) - (P2.getX() -
1335                                P0.getX()) * (P1.getY() - P0.getY());
1336        }
1337
1338        /**
1339         * Compute the minimum size rotated bounding rectangle that contains this
1340         * shape using the rotating calipers approach.
1341         * 
1342         * @see org.openimaj.math.geometry.shape.Shape#minimumBoundingRectangle()
1343         * @see RotatingCalipers#getMinimumBoundingRectangle(Polygon, boolean)
1344         */
1345        @Override
1346        public RotatedRectangle minimumBoundingRectangle() {
1347                return RotatingCalipers.getMinimumBoundingRectangle(this, false);
1348        }
1349
1350        /**
1351         * Compute the minimum size rotated bounding rectangle that contains this
1352         * shape using the rotating calipers approach.
1353         * 
1354         * @see RotatingCalipers#getMinimumBoundingRectangle(Polygon, boolean)
1355         * 
1356         * @param assumeSimple
1357         *            can the algorithm assume the polygon is simple and use an
1358         *            optimised (Melkman's) convex hull?
1359         * 
1360         * @return the minimum bounding box
1361         */
1362        public RotatedRectangle minimumBoundingRectangle(boolean assumeSimple) {
1363                return RotatingCalipers.getMinimumBoundingRectangle(this, assumeSimple);
1364        }
1365}