mpi/fruitfly/registration/CrossCorrelation2D.java

   1 package mpi.fruitfly.registration;
   2
   3 /**
   4  * <p>Title: CrossCorrelation2D</p>
   5  *
   6  * <p>Description: </p>
   7  *
   8  * <p>Copyright: Copyright (c) 2007</p>
   9  *
  10  * <p>Company: </p>
  11  *
  12  * <p>License: GPL
  13  *
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License 2
  16  * as published by the Free Software Foundation.
  17  *
  18  * This program is distributed in the hope that it will be useful,
  19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21  * GNU General Public License for more details.
  22  *
  23  * You should have received a copy of the GNU General Public License
  24  * along with this program; if not, write to the Free Software
  25  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  26  *
  27  * @author Stephan Preibisch
  28  * @version 1.0
  29  */
  30 import ij.io.Opener;
  31 import ij.ImagePlus;
  32 import ij.process.*;
  33
  34 import mpi.fruitfly.math.datastructures.FloatArray2D;
  35 import mpi.fruitfly.general.MultiThreading;
  36
  37 import static mpi.fruitfly.general.ImageArrayConverter.*;
  38 import static mpi.fruitfly.math.General.*;
  39
  40 import java.awt.Point;
  41 import java.util.concurrent.atomic.AtomicInteger;
  42 import java.lang.annotation.Annotation;
  43
  44 public class CrossCorrelation2D
  45 {
  46     public FloatArray2D img1, img2;
  47     public boolean showImages = false;
  48
  49     double maxR = -2;
  50     int displaceX = 0, displaceY = 0;
  51     final Double regCoef = new Double(-2);
  52
  53     public CrossCorrelation2D(String image1, String image2, boolean showImages)
  54     {
  55         // load images
  56         ImagePlus img1 = new Opener().openImage(image1);
  57         ImagePlus img2 = new Opener().openImage(image2);
  58
  59         if (showImages)
  60         {
  61             img1.show();
  62             img2.show();
  63         }
  64
  65         ImageProcessor ip1 = img1.getProcessor();
  66         ImageProcessor ip2 = img2.getProcessor();
  67
  68         this.img1 = ImageToFloatArray2D(ip1);
  69         this.img2 = ImageToFloatArray2D(ip2);
  70         this.showImages = showImages;
  71
  72         //computeCrossCorrelation(ImageToFloatArray2D(ip1), ImageToFloatArray2D(ip2), showImages);
  73     }
  74
  75     public CrossCorrelation2D(ImageProcessor ip1, ImageProcessor ip2, boolean showImages)
  76     {
  77         if (showImages)
  78         {
  79             ImagePlus img1 = new ImagePlus("Image 1", ip1);
  80             ImagePlus img2 = new ImagePlus("Image 2", ip2);
  81
  82             img1.show();
  83             img2.show();
  84         }
  85
  86         this.img1 = ImageToFloatArray2D(ip1);
  87         this.img2 = ImageToFloatArray2D(ip2);
  88         this.showImages = showImages;
  89         //computeCrossCorrelation(ImageToFloatArray2D(ip1), ImageToFloatArray2D(ip2), showImages);
  90     }
  91
  92     /**
  93      * Computes a translational registration with the help of the cross correlation measure. <br>
  94      * Limits the overlap to 30% and restricts the shift furthermore by a factor you can tell him
  95      * (this is useful if you f. ex. know that the vertical shift is much less than the horizontal).
  96      *
  97      * NOTE: Works multithreaded
  98      *
  99      * @param relMinOverlapX double - if you want to scan for less possible translations seen from a direct overlay,
 100      * give the relative factor here (e.g. 0.3 means DONOT scan the outer 30%)
 101      * NOTE: Below 0.05 does not really make sense as you then compare only very few pixels (even one) on the edges
 102      * which gives then an R of 1 (perfect match)
 103      * @param relMinOverlapY double - if you want to scan for less possible translations seen from a direct overlay,
 104      * give the relative factor here (e.g. 0.3 means DONOT scan the outer 30%)
 105      * NOTE: Below 0.05 does not really make sense as you then compare only very few pixels (even one) on the edges
 106      * which gives then an R of 1 (perfect match)
 107      * @param showImages boolean - Show the result of the cross correlation translation
 108      * @return double[] return a double array containing {displaceX, displaceY, R}
 109      */
 110     public double[] computeCrossCorrelationMT(final double relMinOverlapX, final double relMinOverlapY, final boolean showImages)
 111     {
 112         //final double relMinOverlap = 0.30;
 113
 114         final int w1 = img1.width;
 115         final int w2 = img2.width;
 116
 117         final int h1 = img1.height;
 118         final int h2 = img2.height;
 119
 120         final int min_border_w = (int) (w1 < w2 ? w1 * relMinOverlapX + 0.5 : w2 * relMinOverlapX + 0.5);
 121         final int min_border_h = (int) (h1 < h2 ? h1 * relMinOverlapY + 0.5 : h2 * relMinOverlapY + 0.5);
 122
 123         /*System.out.println(w1 + " " + h1 + " " + w2 + " " + h2);
 124         System.out.println(min_border_w + " " + min_border_h);
 125         System.out.println("Y [" + (-h2 + min_border_h) + " , " + (h1 - min_border_h) + "]");*/
 126
 127         final AtomicInteger ai = new AtomicInteger(-w1 + min_border_w);
 128
 129         Thread[] threads = MultiThreading.newThreads();
 130
 131         for (int ithread = 0; ithread < threads.length; ++ithread)
 132             threads[ithread] = new Thread(new Runnable()
 133             {
 134                 public void run()
 135                 {
 136                     for (int moveX = ai.getAndIncrement(); moveX < w2 - min_border_w; moveX = ai.getAndIncrement())
 137                     {
 138                         for (int moveY = -h1 + min_border_h; moveY < h2 - min_border_h; moveY++)
 139                         {
 140                             // compute average
 141                             double avg1 = 0, avg2 = 0;
 142                             int count = 0;
 143                             double value1, value2;
 144
 145                             // iterate over the area which overlaps
 146
 147                             // if moveX < 0 we have to start just at -moveX because otherwise x2 is negative....
 148                             // same for moveX > 0 and the right side
 149
 150                             //for (int x1 = 0; x1 < w1; x1++)
 151                             for (int x1 = -min(0, moveX); x1 < min(w1, w2 - moveX); x1++)
 152                             {
 153                                 int x2 = x1 + moveX;
 154
 155                                 /*if (x2 < 0 || x2 > w2 - 1)
 156                                     continue;*/
 157
 158                                 //for (int y1 = 0; y1 < h1; y1++)
 159                                 for (int y1 = -min(0, moveY); y1 < min(h1, h2 - moveY); y1++)
 160                                 {
 161                                     int y2 = y1 + moveY;
 162
 163                                     /*if (y2 < 0 || y2 > h2 - 1)
 164                                         continue;*/
 165
 166                                     value1 = img1.get(x1, y1);
 167                                     if (value1 == -1)
 168                                         continue;
 169
 170                                     value2 = img2.get(x2, y2);
 171                                     if (value2 == -1)
 172                                         continue;
 173
 174                                     avg1 += value1;
 175                                     avg2 += value2;
 176                                     count++;
 177
 178                                 }
 179                             }
 180
 181                             if (0 == count)continue;
 182
 183                             avg1 /= (double) count;
 184                             avg2 /= (double) count;
 185
 186                             double var1 = 0, var2 = 0;
 187                             double coVar = 0;
 188                             double dist1, dist2;
 189
 190                             // compute variances and co-variance
 191                             //for (int x1 = 0; x1 < w1; x1++)
 192                             for (int x1 = -min(0, moveX); x1 < min(w1, w2 - moveX); x1++)
 193                             {
 194                                 int x2 = x1 + moveX;
 195
 196                                 /*if (x2 < 0 || x2 > w2 - 1)
 197                                     continue;*/
 198
 199                                 //for (int y1 = 0; y1 < h1; y1++)
 200                                 for (int y1 = -min(0, moveY); y1 < min(h1, h2 - moveY); y1++)
 201                                 {
 202                                     int y2 = y1 + moveY;
 203
 204                                     /*if (y2 < 0 || y2 > h2 - 1)
 205                                         continue;*/
 206
 207                                     value1 = img1.get(x1, y1);
 208                                     if (value1 == -1)
 209                                         continue;
 210
 211                                     value2 = img2.get(x2, y2);
 212                                     if (value2 == -1)
 213                                         continue;
 214
 215                                     dist1 = value1 - avg1;
 216                                     dist2 = value2 - avg2;
 217
 218                                     coVar += dist1 * dist2;
 219                                     var1 += dist1 * dist1;
 220                                     var2 += dist2 * dist2;
 221                                 }
 222                             }
 223
 224                             var1 /= (double) count;
 225                             var2 /= (double) count;
 226                             coVar /= (double) count;
 227
 228                             double stDev1 = Math.sqrt(var1);
 229                             double stDev2 = Math.sqrt(var2);
 230
 231                             // compute correlation coeffienct
 232                             double R = coVar / (stDev1 * stDev2);
 233
 234                             compareAndSetR(R, moveX, moveY);
 235
 236                             if (R < -1 || R > 1)
 237                             {
 238                                 System.out.println("BIG ERROR! R =" + R);
 239                             }
 240                         }
 241                         //System.out.println(moveX + " [" + (-w2 + min_border_w) + ", " + (w1 - min_border_w) + "] + best R: " + maxR + " " + displaceX + " " + displaceY);
 242                     }
 243
 244                 }
 245             });
 246         MultiThreading.startAndJoin(threads);
 247
 248         if (showImages)
 249         {
 250             System.out.println( -displaceX + " " + -displaceY + " " + maxR);
 251             FloatArray2D result = drawTranslatedImages(img1, img2, new Point( -displaceX, -displaceY), DRAWTYPE_OVERLAP);
 252             FloatArrayToImagePlus(result, "result", 0, 0).show();
 253         }
 254
 255         return new double[]
 256                 { -displaceX, -displaceY, maxR};
 257     }
 258
 259     /**
 260      * Synchronized method to compare local correlation coefficient to the globally best one and updating
 261      * the global one if necessary
 262      *
 263      * @param R double Correlation Coefficient
 264      * @param moveX int Shift in X direction
 265      * @param moveY int Shift in Y direction
 266      */
 267     private synchronized void compareAndSetR(final double R, final int moveX, final int moveY)
 268     {
 269         if (R > maxR)
 270         {
 271             maxR = R;
 272             displaceX = moveX;
 273             displaceY = moveY;
 274         }
 275
 276     }
 277
 278     /**
 279      * Computes a translational registration with the help of the cross correlation measure. <br>
 280      * Limits the overlap to 30% and restricts the vertical shift furthermore by a factor of 16.
 281      *
 282      * (NOTE: this method is only single threaded, use computeCrossCorrelationMT instead)
 283      *
 284      * @deprecated This method is only single threaded, use computeCrossCorrelationMT instead
 285      *
 286      * @param relMinOverlapX double - if you want to scan for less possible translations seen from a direct overlay,
 287      * give the relative factor here (e.g. 0.3 means DONOT scan the outer 30%)
 288      * NOTE: Below 0.05 does not really make sense as you then compare only very few pixels (even one) on the edges
 289      * which gives then an R of 1 (perfect match)
 290      * @param relMinOverlapY double - if you want to scan for less possible translations seen from a direct overlay,
 291      * give the relative factor here (e.g. 0.3 means DONOT scan the outer 30%)
 292      * NOTE: Below 0.05 does not really make sense as you then compare only very few pixels (even one) on the edges
 293      * which gives then an R of 1 (perfect match)
 294      * @param showImages boolean - Show the result of the cross correlation translation
 295      * @return double[] return a double array containing {displaceX, displaceY, R}
 296      */
 297     @Deprecated
 298     public double[] computeCrossCorrelation(final double relMinOverlapX, final double relMinOverlapY, boolean showImages)
 299     {
 300         double maxR = -2;
 301         int displaceX = 0, displaceY = 0;
 302
 303         int w1 = img1.width;
 304         int w2 = img2.width;
 305
 306         int h1 = img1.height;
 307         int h2 = img2.height;
 308
 309         //int factorY = 16;
 310
 311         final int min_border_w = (int) (w1 < w2 ? w1 * relMinOverlapX + 0.5 : w2 * relMinOverlapX + 0.5);
 312         final int min_border_h = (int) (h1 < h2 ? h1 * relMinOverlapY + 0.5 : h2 * relMinOverlapY + 0.5);
 313
 314         //System.out.println("Y [" + (-h2 + min_border_h)/factor + " , " + (h1 - min_border_h)/factor + "]");
 315
 316         for (int moveX = (-w1 + min_border_w); moveX < (w2 - min_border_w); moveX++)
 317         {
 318             for (int moveY = (-h1 + min_border_h); moveY < (h2 - min_border_h); moveY++)
 319             {
 320                 // compute average
 321                 double avg1 = 0, avg2 = 0;
 322                 int count = 0;
 323                 double value1, value2;
 324
 325                 // iterate over the area which overlaps
 326
 327                 // if moveX < 0 we have to start just at -moveX because otherwise x2 is negative....
 328                 // same for moveX > 0 and the right side
 329
 330                 //for (int x1 = 0; x1 < w1; x1++)
 331                 for (int x1 = -min(0, moveX); x1 < min(w1, w2 - moveX); x1++)
 332                 {
 333                     int x2 = x1 + moveX;
 334
 335                     /*if (x2 < 0 || x2 > w2 - 1)
 336                         continue;*/
 337
 338                     //for (int y1 = 0; y1 < h1; y1++)
 339                     for (int y1 = -min(0, moveY); y1 < min(h1, h2 - moveY); y1++)
 340                     {
 341                         int y2 = y1 + moveY;
 342
 343                         /*if (y2 < 0 || y2 > h2 - 1)
 344                             continue;*/
 345
 346                         value1 = img1.get(x1, y1);
 347                         if (value1 == -1)
 348                             continue;
 349
 350                         value2 = img2.get(x2, y2);
 351                         if (value2 == -1)
 352                             continue;
 353
 354                         avg1 += value1;
 355                         avg2 += value2;
 356                         count++;
 357
 358                     }
 359                 }
 360
 361                 //System.exit(0);
 362
 363                 if (0 == count) continue;
 364
 365                 avg1 /= (double) count;
 366                 avg2 /= (double) count;
 367
 368                 double var1 = 0, var2 = 0;
 369                 double coVar = 0;
 370                 double dist1, dist2;
 371
 372                 // compute variances and co-variance
 373                 //for (int x1 = 0; x1 < w1; x1++)
 374                 for (int x1 = -min(0, moveX); x1 < min(w1, w2 - moveX); x1++)
 375                 {
 376                     int x2 = x1 + moveX;
 377
 378                     /*if (x2 < 0 || x2 > w2 - 1)
 379                         continue;*/
 380
 381                     //for (int y1 = 0; y1 < h1; y1++)
 382                     for (int y1 = -min(0, moveY); y1 < min(h1, h2 - moveY); y1++)
 383                     {
 384                         int y2 = y1 + moveY;
 385
 386                         /*if (y2 < 0 || y2 > h2 - 1)
 387                             continue;*/
 388
 389                         value1 = img1.get(x1, y1);
 390                         if (value1 == -1)
 391                             continue;
 392
 393                         value2 = img2.get(x2, y2);
 394                         if (value2 == -1)
 395                             continue;
 396
 397                         dist1 = value1 - avg1;
 398                         dist2 = value2 - avg2;
 399
 400                         coVar += dist1 * dist2;
 401                         var1 += dist1 * dist1;
 402                         var2 += dist2 * dist2;
 403                     }
 404                 }
 405
 406                 var1 /= (double) count;
 407                 var2 /= (double) count;
 408                 coVar /= (double) count;
 409
 410                 double stDev1 = Math.sqrt(var1);
 411                 double stDev2 = Math.sqrt(var2);
 412
 413                 // compute correlation coeffienct
 414                 double R = coVar / (stDev1 * stDev2);
 415
 416                 if (R > maxR)
 417                 {
 418                     maxR = R;
 419                     displaceX = moveX;
 420                     displaceY = moveY;
 421                 }
 422
 423                 if (R < -1 || R > 1)
 424                 {
 425                     System.out.println("BIG ERROR! R =" + R);
 426                 }
 427             }
 428             //System.out.println(moveX + " [" + (-w2 + min_border_w) + ", " + (w1 - min_border_w) + "] + best R: " + maxR);
 429         }
 430
 431         if (showImages)
 432         {
 433             System.out.println( -displaceX + " " + -displaceY);
 434             FloatArray2D result = drawTranslatedImages(img1, img2, new Point( -displaceX, -displaceY), DRAWTYPE_OVERLAP);
 435             FloatArrayToImagePlus(result, "result", 0, 0).show();
 436         }
 437
 438         return new double[]{-displaceX, -displaceY, maxR};
 439
 440     }
 441 }