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.algorithm; 031 032import java.util.Comparator; 033 034import org.openimaj.image.FImage; 035import org.openimaj.image.analyser.ImageAnalyser; 036import org.openimaj.image.pixel.FValuePixel; 037import org.openimaj.image.processing.algorithm.FourierCorrelation; 038import org.openimaj.math.geometry.shape.Rectangle; 039import org.openimaj.math.util.FloatArrayStatsUtils; 040 041/** 042 * Basic template matching for {@link FImage}s. Template matching is 043 * performed in the frequency domain using an FFT. 044 * <p> 045 * The implementation is heavily inspired by the OpenCV code. 046 * 047 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 048 */ 049public class FourierTemplateMatcher implements ImageAnalyser<FImage> { 050 /** 051 * Different algorithms for comparing templates to images. 052 * 053 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 054 */ 055 public enum Mode { 056 /** 057 * Compute the score at a point as the sum-squared difference between the image 058 * and the template. 059 * 060 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 061 */ 062 SUM_SQUARED_DIFFERENCE { 063 @Override 064 public boolean scoresAscending() { 065 return false; //smaller scores are better 066 } 067 068 @Override 069 public void processCorrelationMap(FImage img, FImage template, FImage corr) { 070 SummedSqAreaTable sum = new SummedSqAreaTable(); 071 img.analyseWith(sum); 072 073 float templateMean = FloatArrayStatsUtils.mean(template.pixels); 074 float templateStdDev = FloatArrayStatsUtils.std(template.pixels); 075 076 float templateNorm = templateStdDev * templateStdDev; 077 float templateSum2 = templateNorm + templateMean * templateMean; 078 079 templateNorm = templateSum2; 080 081 double invArea = 1.0 / ((double)template.width * template.height); 082 templateSum2 /= invArea; 083 templateNorm = (float) Math.sqrt(templateNorm); 084 templateNorm /= Math.sqrt(invArea); 085 086 final float[][] pix = corr.pixels; 087 088 for( int y = 0; y < corr.height; y++ ) { 089 for( int x = 0; x < corr.width; x++ ) { 090 double num = pix[y][x]; 091 double wndSum2 = 0; 092 093 double t = sum.calculateSqSumArea(x, y, x+template.width, y+template.height); 094 wndSum2 += t; 095 096 num = wndSum2 - 2*num + templateSum2; 097 098 pix[y][x] = (float)num; 099 } 100 } 101 } 102 }, 103 /** 104 * Compute the normalised score at a point as the sum-squared difference between the image 105 * and the template. 106 * 107 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 108 */ 109 NORM_SUM_SQUARED_DIFFERENCE { 110 @Override 111 public boolean scoresAscending() { 112 return false; //smaller scores are better 113 } 114 115 @Override 116 public void processCorrelationMap(FImage img, FImage template, FImage corr) { 117 SummedSqAreaTable sum = new SummedSqAreaTable(); 118 img.analyseWith(sum); 119 120 float templateMean = FloatArrayStatsUtils.mean(template.pixels); 121 float templateStdDev = FloatArrayStatsUtils.std(template.pixels); 122 123 float templateNorm = templateStdDev * templateStdDev; 124 float templateSum2 = templateNorm + templateMean * templateMean; 125 126 templateNorm = templateSum2; 127 128 double invArea = 1.0 / ((double)template.width * template.height); 129 templateSum2 /= invArea; 130 templateNorm = (float) Math.sqrt(templateNorm); 131 templateNorm /= Math.sqrt(invArea); 132 133 final float[][] pix = corr.pixels; 134 135 for( int y = 0; y < corr.height; y++ ) { 136 for( int x = 0; x < corr.width; x++ ) { 137 double num = pix[y][x]; 138 double wndMean2 = 0, wndSum2 = 0; 139 140 double t = sum.calculateSqSumArea(x, y, x+template.width, y+template.height); 141 wndSum2 += t; 142 143 num = wndSum2 - 2*num + templateSum2; 144 145 t = Math.sqrt( Math.max(wndSum2 - wndMean2, 0) ) * templateNorm; 146 num /= t; 147 148 pix[y][x] = (float)num; 149 } 150 } 151 } 152 }, 153 /** 154 * Compute the score at a point as the summed product between the image 155 * and the template. 156 * 157 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 158 */ 159 CORRELATION { 160 @Override 161 public boolean scoresAscending() { 162 return true; //bigger scores are better 163 } 164 165 @Override 166 public void processCorrelationMap(FImage img, FImage template, FImage corr) { 167 // Do nothing - image is already 168 } 169 }, 170 /** 171 * Compute the normalised score at a point as the summed product between the image 172 * and the template. 173 * 174 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 175 */ 176 NORM_CORRELATION { 177 @Override 178 public boolean scoresAscending() { 179 return true; //bigger scores are better 180 } 181 182 @Override 183 public void processCorrelationMap(FImage img, FImage template, FImage corr) { 184 SummedSqAreaTable sum = new SummedSqAreaTable(); 185 img.analyseWith(sum); 186 187 float templateMean = FloatArrayStatsUtils.mean(template.pixels); 188 float templateStdDev = FloatArrayStatsUtils.std(template.pixels); 189 190 float templateNorm = templateStdDev * templateStdDev; 191 templateNorm += templateMean * templateMean; 192 193 double invArea = 1.0 / ((double)template.width * template.height); 194 templateNorm = (float) Math.sqrt(templateNorm); 195 templateNorm /= Math.sqrt(invArea); 196 197 final float[][] pix = corr.pixels; 198 199 for( int y = 0; y < corr.height; y++ ) { 200 for( int x = 0; x < corr.width; x++ ) { 201 double num = pix[y][x]; 202 double wndMean2 = 0, wndSum2 = 0; 203 204 double t = sum.calculateSqSumArea(x, y, x+template.width, y+template.height); 205 wndSum2 += t; 206 207 t = Math.sqrt( Math.max(wndSum2 - wndMean2, 0) ) * templateNorm; 208 num /= t; 209 210 pix[y][x] = (float)num; 211 } 212 } 213 } 214 }, 215 /** 216 * Compute the score at a point as the summed product between the mean-centered image patch 217 * and the mean-centered template. 218 * 219 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 220 */ 221 CORRELATION_COEFFICIENT { 222 @Override 223 public boolean scoresAscending() { 224 return true; //bigger scores are better 225 } 226 227 @Override 228 public void processCorrelationMap(FImage img, FImage template, FImage corr) { 229 SummedAreaTable sum = new SummedAreaTable(); 230 img.analyseWith(sum); 231 232 final float templateMean = FloatArrayStatsUtils.mean(template.pixels); //TODO: cache this 233 final float[][] pix = corr.pixels; 234 235 for( int y = 0; y < corr.height; y++ ) { 236 for( int x = 0; x < corr.width; x++ ) { 237 double num = pix[y][x]; 238 double t = sum.calculateArea(x, y, x+template.width, y+template.height); 239 240 num -= t * templateMean; 241 242 pix[y][x] = (float)num; 243 } 244 } 245 } 246 }, 247 /** 248 * Compute the normalised score at a point as the summed product between the mean-centered image patch 249 * and the mean-centered template. 250 * 251 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk) 252 */ 253 NORM_CORRELATION_COEFFICIENT { 254 @Override 255 public boolean scoresAscending() { 256 return true; //bigger scores are better 257 } 258 259 @Override 260 public void processCorrelationMap(FImage img, FImage template, FImage corr) { 261 SummedSqAreaTable sum = new SummedSqAreaTable(); 262 img.analyseWith(sum); 263 264 float templateMean = FloatArrayStatsUtils.mean(template.pixels); 265 float templateStdDev = FloatArrayStatsUtils.std(template.pixels); 266 267 float templateNorm = templateStdDev; 268 269 if( templateNorm == 0 ) 270 { 271 corr.fill(1); 272 return; 273 } 274 275 double invArea = 1.0 / ((double)template.width * template.height); 276 templateNorm /= Math.sqrt(invArea); 277 278 final float[][] pix = corr.pixels; 279 280 for( int y = 0; y < corr.height; y++ ) { 281 for( int x = 0; x < corr.width; x++ ) { 282 double num = pix[y][x]; 283 284 double t = sum.calculateSumArea(x, y, x+template.width, y+template.height); 285 double wndMean2 = t * t * invArea; 286 num -= t * templateMean; 287 288 double wndSum2 = sum.calculateSqSumArea(x, y, x+template.width, y+template.height); 289 290 t = Math.sqrt( Math.max(wndSum2 - wndMean2, 0) ) * templateNorm; 291 num /= t; 292 293 pix[y][x] = (float)num; 294 } 295 } 296 } 297 } 298 ; 299 300 /** 301 * Are the scores ascending (i.e. bigger is better) or descending (smaller is better)? 302 * @return true is bigger scores are better; false if smaller scores are better. 303 */ 304 public abstract boolean scoresAscending(); 305 306 /** 307 * Process the cross-correlation image to the contain the relevant output values for the 308 * chosen mode. 309 * @param img the image 310 * @param template the template 311 * @param corr the cross correlation map produced by {@link FourierCorrelation}. 312 */ 313 public abstract void processCorrelationMap(FImage img, FImage template, FImage corr); 314 } 315 316 private FourierCorrelation correlation; 317 private Mode mode; 318 private Rectangle searchBounds; 319 private FImage responseMap; 320 private int templateWidth; 321 private int templateHeight; 322 323 /** 324 * Default constructor with the template to match. When matching is 325 * performed by {@link #analyseImage(FImage)}, the whole image 326 * will be searched. 327 * 328 * @param template The template. 329 * @param mode The matching mode. 330 */ 331 public FourierTemplateMatcher(FImage template, Mode mode) { 332 this.correlation = new FourierCorrelation(template); 333 this.mode = mode; 334 this.templateWidth = template.width; 335 this.templateHeight = template.height; 336 } 337 338 /** 339 * Construct with the template to match and the bounds rectangle in which 340 * to search. The search bounds rectangle is defined with respect 341 * to the centre of the template. 342 * 343 * @param template The template 344 * @param bounds The bounding box for search. 345 * @param mode The matching mode. 346 */ 347 public FourierTemplateMatcher(FImage template, Rectangle bounds, Mode mode) { 348 this(template, mode); 349 this.searchBounds = bounds; 350 } 351 352 /** 353 * @return the search bound rectangle 354 */ 355 public Rectangle getSearchBounds() { 356 return searchBounds; 357 } 358 359 /** 360 * Set the search bounds rectangle. The search bounds rectangle 361 * is defined with respect to the centre of the template. 362 * Setting to <code>null</code> results in the entire image 363 * being searched. 364 * 365 * @param searchBounds the search bounds to set 366 */ 367 public void setSearchBounds(Rectangle searchBounds) { 368 this.searchBounds = searchBounds; 369 } 370 371 /** 372 * Perform template matching. If a bounds rectangle 373 * is has not been set or is null, then the whole 374 * image will be searched. Otherwise the area of the image 375 * which lies in the previously set search bounds will be 376 * searched. 377 * 378 * @see org.openimaj.image.analyser.ImageAnalyser#analyseImage(org.openimaj.image.Image) 379 */ 380 @Override 381 public void analyseImage(FImage image) { 382 FImage subImage; 383 384 if (this.searchBounds != null) { 385 final int halfWidth = templateWidth / 2; 386 final int halfHeight = templateHeight / 2; 387 388 int x = (int) Math.max(searchBounds.x - halfWidth, 0); 389 int width = (int) searchBounds.width + templateWidth; 390 if (searchBounds.x - halfWidth < 0) { 391 width += (searchBounds.x - halfWidth); 392 } 393 if (x + width > image.width) 394 width = image.width; 395 396 int y = (int) Math.max(searchBounds.y - halfHeight, 0); 397 int height = (int) searchBounds.height + templateHeight; 398 if (searchBounds.y - halfHeight < 0) { 399 height += (searchBounds.y - halfHeight); 400 } 401 if (y + height > image.height) 402 height = image.height; 403 404 //FIXME: this is doing an additional copy; should be rolled into FFT data prep step in FourierTransform 405 subImage = image.extractROI( 406 x, 407 y, 408 width, 409 height 410 ); 411 } else { 412 subImage = image.clone(); 413 } 414 415 responseMap = subImage.process(correlation); 416 responseMap.height = responseMap.height - correlation.template.height + 1; 417 responseMap.width = responseMap.width - correlation.template.width + 1; 418 419 mode.processCorrelationMap(subImage, correlation.template, responseMap); 420 } 421 422 /** 423 * Get the top-N "best" responses found by the template matcher. 424 * 425 * @param numResponses The number of responses 426 * @return the best responses found 427 */ 428 public FValuePixel[] getBestResponses(int numResponses) { 429 Comparator<FValuePixel> comparator = mode.scoresAscending() ? FValuePixel.ReverseValueComparator.INSTANCE : FValuePixel.ValueComparator.INSTANCE; 430 431 return TemplateMatcher.getBestResponses(numResponses, responseMap, getXOffset(), getYOffset(), comparator); 432 } 433 434 /** 435 * @return The x-offset of the top-left of the response map 436 * returned by {@link #getResponseMap()} to the original image 437 * analysed by {@link #analyseImage(FImage)}. 438 */ 439 public int getXOffset() { 440 final int halfWidth = templateWidth / 2; 441 442 if (this.searchBounds == null) 443 return halfWidth; 444 else 445 return (int) Math.max(searchBounds.x - halfWidth, halfWidth); 446 } 447 448 /** 449 * @return The y-offset of the top-left of the response map 450 * returned by {@link #getResponseMap()} to the original image 451 * analysed by {@link #analyseImage(FImage)}. 452 */ 453 public int getYOffset() { 454 final int halfHeight = templateHeight / 2; 455 456 if (this.searchBounds == null) 457 return halfHeight; 458 else 459 return (int) Math.max(searchBounds.y - halfHeight, halfHeight); 460 } 461 462 /** 463 * @return The responseMap generated from the last call to {@link #analyseImage(FImage)} 464 */ 465 public FImage getResponseMap() { 466 return responseMap; 467 } 468}