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 odk.lang.FastMath; 033 034import org.openimaj.math.geometry.point.Point2d; 035import org.openimaj.math.geometry.point.Point2dImpl; 036import org.openimaj.math.geometry.shape.Polygon; 037import org.openimaj.math.geometry.shape.RotatedRectangle; 038 039/** 040 * Rotating calipers algorithm, based on the implementation by <a 041 * 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, lengthAB, lengthBC, 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 } 154 else if (thetaJ <= thetaK && thetaJ <= thetaL) { 155 return thetaJ; 156 } 157 else if (thetaK <= thetaL) { 158 return thetaK; 159 } 160 else { 161 return thetaL; 162 } 163 } 164 165 protected static int getIndex(Polygon convexHull, Corner corner) { 166 167 int index = 0; 168 final Point2d point = convexHull.points.get(index); 169 float px = point.getX(); 170 float py = point.getY(); 171 172 for (int i = 1; i < convexHull.size() - 1; i++) { 173 174 final Point2d temp = convexHull.points.get(i); 175 boolean change = false; 176 177 final float tx = temp.getX(); 178 final float ty = temp.getY(); 179 180 switch (corner) { 181 case UPPER_RIGHT: 182 change = (tx > px || (tx == px && ty > py)); 183 break; 184 case UPPER_LEFT: 185 change = (ty > py || (ty == py && tx < px)); 186 break; 187 case LOWER_LEFT: 188 change = (tx < px || (tx == px && ty < py)); 189 break; 190 case LOWER_RIGHT: 191 change = (ty < py || (ty == py && tx > px)); 192 break; 193 } 194 195 if (change) { 196 index = i; 197 px = tx; 198 py = ty; 199 } 200 } 201 202 return index; 203 } 204 205 protected static class Caliper { 206 207 final static double SIGMA = 0.00000000001; 208 209 final Polygon convexHull; 210 int pointIndex; 211 double currentAngle; 212 213 Caliper(Polygon convexHull, int pointIndex, double currentAngle) { 214 this.convexHull = convexHull; 215 this.pointIndex = pointIndex; 216 this.currentAngle = currentAngle; 217 } 218 219 double getAngleNextPoint() { 220 final Point2d p1 = convexHull.get(pointIndex); 221 final Point2d p2 = convexHull.get((pointIndex + 1) % convexHull.size()); 222 223 final double deltaX = p2.getX() - p1.getX(); 224 final double deltaY = p2.getY() - p1.getY(); 225 226 final double angle = FastMath.atan2(deltaY, deltaX); 227 228 return angle < 0 ? ANGLE_360DEG_IN_RADS + angle : angle; 229 } 230 231 double getConstant() { 232 233 final Point2d p = convexHull.get(pointIndex); 234 235 return p.getY() - (getSlope() * p.getX()); 236 } 237 238 double getDeltaAngleNextPoint() { 239 240 double angle = getAngleNextPoint(); 241 242 angle = angle < 0 ? ANGLE_360DEG_IN_RADS + angle - currentAngle : angle - currentAngle; 243 244 return angle < 0 ? ANGLE_360DEG_IN_RADS : angle; 245 } 246 247 Point2dImpl getIntersection(Caliper that) { 248 249 // the x-intercept of 'this' and 'that': x = ((c2 - c1) / (m1 - m2)) 250 double x; 251 // the y-intercept of 'this' and 'that', given 'x': (m*x) + c 252 double y; 253 254 if (this.isVertical()) { 255 x = convexHull.points.get(pointIndex).getX(); 256 } 257 else if (this.isHorizontal()) { 258 x = that.convexHull.points.get(that.pointIndex).getX(); 259 } 260 else { 261 x = (that.getConstant() - this.getConstant()) / (this.getSlope() - that.getSlope()); 262 } 263 264 if (this.isVertical()) { 265 y = that.getConstant(); 266 } 267 else if (this.isHorizontal()) { 268 y = this.getConstant(); 269 } 270 else { 271 y = (this.getSlope() * x) + this.getConstant(); 272 } 273 274 return new Point2dImpl(x, y); 275 } 276 277 double getSlope() { 278 return Math.tan(currentAngle); 279 } 280 281 boolean isHorizontal() { 282 return (Math.abs(currentAngle) < SIGMA) || (Math.abs(currentAngle - ANGLE_180DEG_IN_RADS) < SIGMA); 283 } 284 285 boolean isVertical() { 286 return (Math.abs(currentAngle - ANGLE_90DEG_IN_RADS) < SIGMA) 287 || (Math.abs(currentAngle - ANGLE_270DEG_IN_RADS) < SIGMA); 288 } 289 290 void rotateBy(double angle) { 291 292 if (this.getDeltaAngleNextPoint() == angle) { 293 pointIndex++; 294 } 295 296 this.currentAngle += angle; 297 } 298 } 299}