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.util; 031 032import org.openimaj.math.geometry.point.Point2d; 033import org.openimaj.math.geometry.point.Point2dImpl; 034import org.openimaj.math.geometry.shape.Polygon; 035import org.openimaj.math.geometry.shape.RotatedRectangle; 036 037import net.jafama.FastMath; 038 039/** 040 * Rotating calipers algorithm, based on the implementation by 041 * <a href="https://github.com/bkiers/RotatingCalipers">Bart Kiers</a>. 042 * <p> 043 * Modified to only use radians for angles and fit better within the OpenIMAJ 044 * geometry classes. 045 * 046 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 047 * @author Bart Kiers 048 * 049 */ 050public final class RotatingCalipers { 051 private static final double ANGLE_0DEG_IN_RADS = 0; 052 private static final double ANGLE_90DEG_IN_RADS = Math.PI / 2; 053 private static final double ANGLE_180DEG_IN_RADS = Math.PI; 054 private static final double ANGLE_270DEG_IN_RADS = 3 * Math.PI / 2; 055 private static final double ANGLE_360DEG_IN_RADS = 2 * Math.PI; 056 057 protected enum Corner { 058 UPPER_RIGHT, UPPER_LEFT, LOWER_LEFT, LOWER_RIGHT 059 } 060 061 private static double getArea(Point2dImpl[] rectangle) { 062 063 final double deltaXAB = rectangle[0].x - rectangle[1].x; 064 final double deltaYAB = rectangle[0].y - rectangle[1].y; 065 066 final double deltaXBC = rectangle[1].x - rectangle[2].x; 067 final double deltaYBC = rectangle[1].y - rectangle[2].y; 068 069 final double lengthAB = FastMath.sqrt((deltaXAB * deltaXAB) + (deltaYAB * deltaYAB)); 070 final double lengthBC = FastMath.sqrt((deltaXBC * deltaXBC) + (deltaYBC * deltaYBC)); 071 072 return lengthAB * lengthBC; 073 } 074 075 /** 076 * Use the rotating calipers algorithm to optimally find the minimum sized 077 * rotated rectangle that encompasses the outer shell of the given polygon. 078 * 079 * @param poly 080 * the polygon 081 * @param assumeSimple 082 * can the algorithm assume the polygon is simple and use an 083 * optimised (Melkman's) convex hull? 084 * @return the minimum enclosing rectangle 085 */ 086 public static RotatedRectangle getMinimumBoundingRectangle(Polygon poly, boolean assumeSimple) { 087 final Polygon convexHull = assumeSimple ? poly.calculateConvexHullMelkman() : poly.calculateConvexHull(); 088 089 if (convexHull.size() < 3) { 090 // FIXME 091 return new RotatedRectangle(convexHull.calculateRegularBoundingBox(), 0); 092 } 093 094 Point2dImpl[] minimum = null; 095 double minimumAngle = 0; 096 double area = Double.MAX_VALUE; 097 098 final Caliper I = new Caliper(convexHull, getIndex(convexHull, Corner.UPPER_RIGHT), ANGLE_90DEG_IN_RADS); 099 final Caliper J = new Caliper(convexHull, getIndex(convexHull, Corner.UPPER_LEFT), ANGLE_180DEG_IN_RADS); 100 final Caliper K = new Caliper(convexHull, getIndex(convexHull, Corner.LOWER_LEFT), ANGLE_270DEG_IN_RADS); 101 final Caliper L = new Caliper(convexHull, getIndex(convexHull, Corner.LOWER_RIGHT), ANGLE_0DEG_IN_RADS); 102 103 while (L.currentAngle < ANGLE_90DEG_IN_RADS) { 104 final Point2dImpl[] rectangle = new Point2dImpl[] { 105 L.getIntersection(I), 106 I.getIntersection(J), 107 J.getIntersection(K), 108 K.getIntersection(L) 109 }; 110 111 final double tempArea = getArea(rectangle); 112 if (minimum == null || tempArea < area) { 113 minimum = rectangle; 114 minimumAngle = L.currentAngle; 115 area = tempArea; 116 } 117 118 final double smallestTheta = getSmallestTheta(I, J, K, L); 119 I.rotateBy(smallestTheta); 120 J.rotateBy(smallestTheta); 121 K.rotateBy(smallestTheta); 122 L.rotateBy(smallestTheta); 123 } 124 125 return makeRotated(minimum, minimumAngle); 126 } 127 128 private static RotatedRectangle makeRotated(Point2dImpl[] rectangle, double angle) { 129 final double deltaXAB = rectangle[0].x - rectangle[1].x; 130 final double deltaYAB = rectangle[0].y - rectangle[1].y; 131 132 final double deltaXBC = rectangle[1].x - rectangle[2].x; 133 final double deltaYBC = rectangle[1].y - rectangle[2].y; 134 135 final double lengthAB = FastMath.sqrt((deltaXAB * deltaXAB) + (deltaYAB * deltaYAB)); 136 final double lengthBC = FastMath.sqrt((deltaXBC * deltaXBC) + (deltaYBC * deltaYBC)); 137 138 final double cx = (rectangle[0].x + rectangle[1].x + rectangle[2].x + rectangle[3].x) / 4; 139 final double cy = (rectangle[0].y + rectangle[1].y + rectangle[2].y + rectangle[3].y) / 4; 140 141 return new RotatedRectangle(cx, cy, lengthBC, lengthAB, angle); 142 } 143 144 private static double getSmallestTheta(Caliper I, Caliper J, Caliper K, Caliper L) { 145 146 final double thetaI = I.getDeltaAngleNextPoint(); 147 final double thetaJ = J.getDeltaAngleNextPoint(); 148 final double thetaK = K.getDeltaAngleNextPoint(); 149 final double thetaL = L.getDeltaAngleNextPoint(); 150 151 if (thetaI <= thetaJ && thetaI <= thetaK && thetaI <= thetaL) { 152 return thetaI; 153 } else if (thetaJ <= thetaK && thetaJ <= thetaL) { 154 return thetaJ; 155 } else if (thetaK <= thetaL) { 156 return thetaK; 157 } else { 158 return thetaL; 159 } 160 } 161 162 protected static int getIndex(Polygon convexHull, Corner corner) { 163 164 int index = 0; 165 final Point2d point = convexHull.points.get(index); 166 float px = point.getX(); 167 float py = point.getY(); 168 169 for (int i = 1; i < convexHull.size() - 1; i++) { 170 171 final Point2d temp = convexHull.points.get(i); 172 boolean change = false; 173 174 final float tx = temp.getX(); 175 final float ty = temp.getY(); 176 177 switch (corner) { 178 case UPPER_RIGHT: 179 change = (tx > px || (tx == px && ty > py)); 180 break; 181 case UPPER_LEFT: 182 change = (ty > py || (ty == py && tx < px)); 183 break; 184 case LOWER_LEFT: 185 change = (tx < px || (tx == px && ty < py)); 186 break; 187 case LOWER_RIGHT: 188 change = (ty < py || (ty == py && tx > px)); 189 break; 190 } 191 192 if (change) { 193 index = i; 194 px = tx; 195 py = ty; 196 } 197 } 198 199 return index; 200 } 201 202 protected static class Caliper { 203 204 final static double SIGMA = 0.00000000001; 205 206 final Polygon convexHull; 207 int pointIndex; 208 double currentAngle; 209 210 Caliper(Polygon convexHull, int pointIndex, double currentAngle) { 211 this.convexHull = convexHull; 212 this.pointIndex = pointIndex; 213 this.currentAngle = currentAngle; 214 } 215 216 double getAngleNextPoint() { 217 final Point2d p1 = convexHull.get(pointIndex); 218 final Point2d p2 = convexHull.get((pointIndex + 1) % convexHull.size()); 219 220 final double deltaX = p2.getX() - p1.getX(); 221 final double deltaY = p2.getY() - p1.getY(); 222 223 final double angle = FastMath.atan2(deltaY, deltaX); 224 225 return angle < 0 ? ANGLE_360DEG_IN_RADS + angle : angle; 226 } 227 228 double getConstant() { 229 230 final Point2d p = convexHull.get(pointIndex); 231 232 return p.getY() - (getSlope() * p.getX()); 233 } 234 235 double getDeltaAngleNextPoint() { 236 237 double angle = getAngleNextPoint(); 238 239 angle = angle < 0 ? ANGLE_360DEG_IN_RADS + angle - currentAngle : angle - currentAngle; 240 241 return angle < 0 ? ANGLE_360DEG_IN_RADS : angle; 242 } 243 244 Point2dImpl getIntersection(Caliper that) { 245 246 // the x-intercept of 'this' and 'that': x = ((c2 - c1) / (m1 - m2)) 247 double x; 248 // the y-intercept of 'this' and 'that', given 'x': (m*x) + c 249 double y; 250 251 if (this.isVertical()) { 252 x = convexHull.points.get(pointIndex).getX(); 253 } else if (this.isHorizontal()) { 254 x = that.convexHull.points.get(that.pointIndex).getX(); 255 } else { 256 x = (that.getConstant() - this.getConstant()) / (this.getSlope() - that.getSlope()); 257 } 258 259 if (this.isVertical()) { 260 y = that.getConstant(); 261 } else if (this.isHorizontal()) { 262 y = this.getConstant(); 263 } else { 264 y = (this.getSlope() * x) + this.getConstant(); 265 } 266 267 return new Point2dImpl(x, y); 268 } 269 270 double getSlope() { 271 return Math.tan(currentAngle); 272 } 273 274 boolean isHorizontal() { 275 return (Math.abs(currentAngle) < SIGMA) || (Math.abs(currentAngle - ANGLE_180DEG_IN_RADS) < SIGMA); 276 } 277 278 boolean isVertical() { 279 return (Math.abs(currentAngle - ANGLE_90DEG_IN_RADS) < SIGMA) 280 || (Math.abs(currentAngle - ANGLE_270DEG_IN_RADS) < SIGMA); 281 } 282 283 void rotateBy(double angle) { 284 285 if (this.getDeltaAngleNextPoint() == angle) { 286 pointIndex++; 287 } 288 289 this.currentAngle += angle; 290 } 291 } 292}