001/******************************************************************************* 002 * This software is provided as a supplement to the authors' textbooks on digital 003 * image processing published by Springer-Verlag in various languages and editions. 004 * Permission to use and distribute this software is granted under the BSD 2-Clause 005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009package imagingbook.common.hough; 010 011import ij.process.ByteProcessor; 012import ij.process.FloatProcessor; 013import imagingbook.common.geometry.basic.Pnt2d; 014import imagingbook.common.util.ParameterBundle; 015import imagingbook.common.util.progress.ProgressReporter; 016 017import java.util.ArrayList; 018import java.util.Collections; 019import java.util.List; 020 021/** 022 * <p> 023 * This class implements the "classic" Hough Transform for straight lines. This implementation improves the algorithm 024 * described in the textbook in the way the accumulator is updated. Here, the quantity 1 is not added to a single 025 * accumulator cell but gets distributed over two neighboring (radial) cells to reduce aliasing effects. Thus we 026 * accumulate non-integer values and therefore the various accumulators are of type {@code float[][]}. See Sec. 12.2 of 027 * [1] for additional details. 028 * </p> 029 * <p> 030 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 031 * (2022). 032 * </p> 033 * 034 * @author WB 035 * @version 2022/09/10 removed ImageJ progress reporting 036 */ 037public class HoughTransformLines implements ProgressReporter { 038 039 // TODO: revise constructors and parameters, add bias correction 040 041 public static class Parameters implements ParameterBundle<HoughTransformLines> { 042 /** Number of angular steps over [0, pi] */ 043 public int nAng = 256; 044 /** Number of radial steps in each pos/neg direction (accum. size = 2 * nRad + 1) */ 045 public int nRad = 128; 046 } 047 048 private final int nAng; // number of angular steps over [0, pi] 049 private final int nRad; // number of radial steps in each pos/neg direction 050 051 private final int width, height; // size of the reference frame (image) 052 private final double xRef, yRef; // reference point (x/y-coordinate of image center) 053 054 private final double dAng; // increment of angle 055 private final double dRad; // increment of radius 056 private final int cRad; // array index representing the zero radius 057 058 private final int accWidth; // width of the accumulator array (angular direction) 059 private final int accHeight; // height of the accumulator array (radial direction) 060 private final float[][] accumulator; // accumulator array 061 private final float[][] accumulatorExt; // extended accumulator array 062 private final float[][] accumulatorMax; // accumulator, with local maxima only 063 064 private final double[] cosTable; // tabulated cosine values 065 private final double[] sinTable; // tabulated sine values 066 067 // -------------- public constructor(s) ------------------------ 068 069 /** 070 * Constructor, creates a new Hough transform from the binary image. 071 * 072 * @param I input image, relevant (edge) points have pixel values greater 0. 073 * @param params parameter object. 074 */ 075 public HoughTransformLines(ByteProcessor I, Parameters params) { 076 this(I.getWidth(), I.getHeight(), params); 077 this.process(I, accumulator); 078 } 079 080 public HoughTransformLines(ByteProcessor I) { 081 this(I, new Parameters()); 082 } 083 084 /** 085 * Constructor, creates a new Hough transform from a sequence of 2D points. Parameters M, N are only used to specify 086 * the reference point (usually at the center of the image). Use this constructor if the relevant image points are 087 * collected separately. 088 * 089 * @param points an array of 2D points. 090 * @param width width of the corresponding image plane. 091 * @param height height of the corresponding image plane. 092 * @param params parameter object. 093 */ 094 public HoughTransformLines(Pnt2d[] points, int width, int height, Parameters params) { 095 this(width, height, params); 096 this.process(points, accumulator); 097 } 098 099 // Non-public constructor used by public constructors (to initialize all final 100 // members variables). 101 private HoughTransformLines(int width, int height, Parameters params) { 102 if (params == null) 103 params = new Parameters(); 104 this.width = width; 105 this.height = height; 106 this.xRef = width / 2; // integer value 107 this.yRef = height / 2; // integer value 108 this.nAng = params.nAng; 109 this.nRad = params.nRad; 110 this.dAng = Math.PI / nAng; 111 this.dRad = 0.5 * Math.hypot(width, height) / nRad; // nRad radial steps over half the diagonal length 112 this.cRad = nRad; 113 this.accWidth = nAng; 114 this.accHeight = nRad + 1 + nRad; 115 this.accumulator = new float[accWidth][accHeight]; // cells are initialized to zero! 116 this.accumulatorExt = new float[2 * accWidth][accHeight]; 117 this.accumulatorMax = new float[accWidth][accHeight]; 118 this.cosTable = makeCosTable(); 119 this.sinTable = makeSinTable(); 120 } 121 122 // -------------- public methods ------------------------ 123 124 /** 125 * Finds and returns the parameters of the strongest lines with a specified min. pixel count. All objects in the 126 * returned array are valid, but the array may be empty. Note: Could perhaps be implemented more efficiently with 127 * insert-sort. 128 * 129 * @param amin the minimum accumulator value for each line. 130 * @param maxLines maximum number of (strongest) lines to extract. 131 * @return a possibly empty array of {@link HoughLine} objects. 132 */ 133 public HoughLine[] getLines(int amin, int maxLines) { 134 makeAccumulatorExt(); // TODO: should not be here 135 findLocalMaxima(); 136 // collect all lines corresponding to accumulator maxima 137 List<HoughLine> lines = new ArrayList<>(); 138 for (int ri = 0; ri < accHeight; ri++) { 139 for (int ai = 0; ai < accWidth; ai++) { 140 int hcount = (int) accumulatorMax[ai][ri]; 141 if (hcount >= amin) { 142 double angle = angleFromIndex(ai); 143 double radius = radiusFromIndex(ri); 144 lines.add(new HoughLine(angle, radius, xRef, yRef, hcount)); 145 } 146 } 147 } 148 // sort 'lines' by count (highest counts first) 149 Collections.sort(lines); 150 List<HoughLine> slines = lines.subList(0, Math.min(maxLines, lines.size())); 151 // convert the list to an array and return: 152 return slines.toArray(new HoughLine[0]); 153 } 154 155 /** 156 * Returns the reference points x-coordinate. 157 * @return as described. 158 */ 159 public double getXref() { 160 return xRef; 161 } 162 163 /** 164 * Returns the reference points y-coordinate. 165 * @return as described. 166 */ 167 public double getYref() { 168 return yRef; 169 } 170 171 /** 172 * Calculates the actual angle (in radians) for angle index {@code ai} 173 * 174 * @param ai angle index [0,...,nAng-1] 175 * @return Angle [0,...,PI] for angle index ai 176 */ 177 public double angleFromIndex(int ai) { 178 return ai * dAng; 179 } 180 181 /** 182 * Calculates the actual radius for radius index ri. 183 * 184 * @param ri radius index [0,...,nRad-1]. 185 * @return Radius [-maxRad,...,maxRad] with respect to reference point (xc, yc). 186 */ 187 public double radiusFromIndex(int ri) { 188 return (ri - cRad) * dRad; 189 } 190 191 public double radiusToIndex(double rad) { 192 return cRad + rad / dRad; 193 } 194 195 // --------------------------------------------------------------- 196 197 /** 198 * Returns the accumulator as a 2D float-array. 199 * @return the contents of the accumulator 200 */ 201 public float[][] getAccumulator() { 202 return this.accumulator; 203 } 204 205 /** 206 * Returns the local maximum values of the accumulator as a 2D float-array (all non-maximum elements are set to 207 * zero). 208 * 209 * @return the local maximum values of the accumulator 210 */ 211 public float[][] getAccumulatorMax() { 212 return this.accumulatorMax; 213 } 214 215 /** 216 * Returns the extended accumulator as a 2D float-array. 217 * @return the extended accumulator 218 */ 219 public float[][] getAccumulatorExtended() { 220 return this.accumulatorExt; 221 } 222 223 /** 224 * Creates and returns an image of the 2D accumulator array. 225 * @return a {@link FloatProcessor} image of the accumulator 226 */ 227 public FloatProcessor getAccumulatorImage() { 228 return new FloatProcessor(accumulator); 229 } 230 231 /** 232 * Creates and returns an image of the extended 2D accumulator array, produced by adding a vertically mirrored copy 233 * of the accumulator to its right end. 234 * 235 * @return a {@link FloatProcessor} image of the extended accumulator 236 */ 237 public FloatProcessor getAccumulatorImageExtended() { 238 FloatProcessor fp = new FloatProcessor(accumulatorExt); 239 return fp; 240 } 241 242 /** 243 * Creates and returns an image of the local maxima of the 2D accumulator array. 244 * @return a {@link FloatProcessor} image of the accumulator maxima 245 */ 246 public FloatProcessor getAccumulatorMaxImage() { 247 return new FloatProcessor(accumulatorMax); 248 } 249 250 // -------------- nonpublic methods ------------------------ 251 252 private double[] makeCosTable() { 253 double[] cosTab = new double[nAng]; 254 for (int ai = 0; ai < nAng; ai++) { 255 double angle = dAng * ai; 256 cosTab[ai] = Math.cos(angle); 257 } 258 return cosTab; 259 } 260 261 private double[] makeSinTable() { 262 double[] sinTab = new double[nAng]; 263 for (int ai = 0; ai < nAng; ai++) { 264 double angle = dAng * ai; 265 sinTab[ai] = Math.sin(angle); 266 } 267 return sinTab; 268 } 269 270 private void process(ByteProcessor ip, float[][] acc) { 271 for (int v = 0; v < height; v++) { 272 for (int u = 0; u < width; u++) { 273 if (ip.get(u, v) != 0) { // this is a foreground (edge) pixel - use ImageAccessor?? 274 processPoint(u, v, acc); 275 } 276 } 277 } 278 } 279 280 private void process(Pnt2d[] points, float[][] acc) { 281 for (int i = 0; i < points.length; i++) { 282 Pnt2d p = points[i]; 283 if (p != null) { 284 processPoint(p.getX(), p.getY(), acc); 285 } 286 } 287 } 288 289 private void processPoint(double u, double v, float[][] acc) { 290 final double x = u - xRef; 291 final double y = v - yRef; 292 for (int ai = 0; ai < accWidth; ai++) { 293// double theta = dAng * ai; 294// double r = x * Math.cos(theta) + y * Math.sin(theta); 295 double r = x * cosTable[ai] + y * sinTable[ai]; // sin/cos tables improve speed! 296 double ri = radiusToIndex(r); 297 // accumulated quantity (1.0) is distributed to 2 neighboring bins: 298 int r0 = (int) Math.floor(ri); // lower radial bin index 299 int r1 = r0 + 1; // upper radial bin index 300 if (r0 >= 0 && r1 < accHeight) { 301 double alpha = ri - r0; 302 acc[ai][r0] += (1.0 - alpha); 303 acc[ai][r1] += alpha; 304 } 305 } 306 } 307 308 /** 309 * Searches for local maxima in the extended accumulator but enters their positions in 'accumulatorMax'. 310 */ 311 private void findLocalMaxima() { 312 for (int ai = 1; ai <= accWidth; ai++) { // note the range! 313 for (int ri = 1; ri < accHeight - 1; ri++) { 314 float vC = accumulatorExt[ai][ri]; // center value 315 boolean ismax = 316 vC > accumulatorExt[ai + 1][ri] && // 0 317 vC > accumulatorExt[ai + 1][ri - 1] && // 1 318 vC > accumulatorExt[ai] [ri - 1] && // 2 319 vC > accumulatorExt[ai - 1][ri - 1] && // 3 320 vC > accumulatorExt[ai - 1][ri] && // 4 321 vC > accumulatorExt[ai - 1][ri + 1] && // 5 322 vC > accumulatorExt[ai] [ri + 1] && // 6 323 vC > accumulatorExt[ai + 1][ri + 1]; // 7 324 if (ismax) { 325 if (ai < accWidth) 326 accumulatorMax[ai % accWidth][ri] = vC; // take care of ai == accWidth 327 else 328 accumulatorMax[ai % accWidth][accHeight - ri - 1] = vC; 329 } 330 } 331 } 332 } 333 334 /** 335 * Creates the extended 2D accumulator array 336 */ 337 private void makeAccumulatorExt() { 338 for (int ai = 0; ai < accWidth; ai++) { 339 for (int ri = 0; ri < accHeight; ri++) { 340 // insert original accumulator into the left side 341 accumulatorExt[ai][ri] = accumulator[ai][ri]; 342 // insert a vertically flipped copy of accumulator into the right side 343 accumulatorExt[accWidth + ai][ri] = accumulator[ai][accHeight - ri - 1]; 344 } 345 } 346 } 347 348 @Override 349 public double getProgress() { 350 // TODO report progress state 351 return 1; 352 } 353 354} 355 356