mpi/fruitfly/registration/ImageFilter.java

   1 package mpi.fruitfly.registration;
   2
   3 /**
   4  * <p>Title: ImageFilter</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
  31 import mpi.fruitfly.math.datastructures.*;
  32 import static mpi.fruitfly.math.General.*;
  33 import java.util.Random;
  34
  35 public class ImageFilter
  36 {
  37     /**
  38      * Does Kaiser-Bessel-Windowing to prevent the fourier spectra from getting infinite numbers,
  39      * the border will be faded to black.
  40      *
  41      * @param FloatArray2D  image Input image stored in row-major order
  42      * @param inverse If true, fading to black in the middle
  43      * @return FloatArray2D The resulting image
  44      */
  45     public static void filterKaiserBessel(FloatArray2D img, boolean inverse)
  46     {
  47         int x, y;
  48         int w = img.width;
  49         int h = img.height;
  50
  51         float twoPiDivWidth = (float)(2.0 * Math.PI) / (float)w;
  52         float twoPiDivHeight = (float)(2.0 * Math.PI) / (float)h;
  53         float kb, val;
  54
  55         // pre-compute filtering values
  56         float[] kbx = new float[w];
  57         float[] kby = new float[h];
  58
  59         for (x = 0; x < w; x++)
  60             kbx[x] = (float)(0.40243 - (0.49804 * Math.cos(twoPiDivWidth * (float) x)) +
  61                      (0.09831 * Math.cos(twoPiDivWidth * 2.0 * (float) x)) -
  62                      (0.00122 * Math.cos(twoPiDivWidth * 3.0 * (float) x)));
  63
  64         for (y = 0; y < h; y++)
  65             kby[y] = (float)(0.40243 - (0.49804 * Math.cos(twoPiDivHeight * (float) y )) +
  66                      (0.09831 * Math.cos(twoPiDivHeight * 2.0 * (float) y )) -
  67                      (0.00122 * Math.cos(twoPiDivHeight * 3.0 * (float) y )));
  68
  69         for (y = 0; y < h; y++)
  70             for (x = 0; x < w; x++)
  71             {
  72                 val = (float) img.get(x,y);
  73                 kb = kbx[x];
  74
  75                 if (inverse)
  76                     kb = 1.0f - kb;
  77
  78                 img.set(kb * val, x, y);
  79             }
  80
  81         for (x = 0; x < w; x++)
  82             for (y = 0; y < h; y++)
  83             {
  84                 val = (float) img.get(x,y);
  85                 kb = kby[y];
  86
  87                 if (inverse)
  88                     kb = 1.0f - kb;
  89
  90                 img.set(kb * val, x, y);
  91             }
  92     }
  93
  94     public static void exponentialWindow(FloatArray2D img)
  95     {
  96         double a = 1000;
  97
  98         // create lookup table
  99         double weightsX[] = new double[img.width];
 100         double weightsY[] = new double[img.height];
 101
 102         for (int x = 0; x < img.width; x++)
 103         {
 104             double relPos = (double)x / (double)(img.width-1);
 105
 106             if (relPos <= 0.5)
 107                 weightsX[x] = 1.0-(1.0/(Math.pow(a,(relPos*2))));
 108             else
 109                 weightsX[x] = 1.0-(1.0/(Math.pow(a,((1-relPos)*2))));
 110         }
 111
 112         for (int y = 0; y < img.height; y++)
 113         {
 114             double relPos = (double)y / (double)(img.height-1);
 115
 116             if (relPos <= 0.5)
 117                 weightsY[y] = 1.0-(1.0/(Math.pow(a,(relPos*2))));
 118             else
 119                 weightsY[y] = 1.0-(1.0/(Math.pow(a,((1-relPos)*2))));
 120         }
 121
 122
 123         for (int y = 0; y < img.height; y++)
 124             for (int x = 0; x < img.width; x++)
 125                 img.set((float) (img.get(x, y) * weightsX[x] * weightsY[y]), x, y);
 126     }
 127
 128     public static void exponentialWindow(FloatArray3D img)
 129     {
 130         double a = 1000;
 131
 132         // create lookup table
 133         double weightsX[] = new double[img.width];
 134         double weightsY[] = new double[img.height];
 135         double weightsZ[] = new double[img.depth];
 136
 137         for (int x = 0; x < img.width; x++)
 138         {
 139             double relPos = (double)x / (double)(img.width-1);
 140
 141             if (relPos <= 0.5)
 142                 weightsX[x] = 1.0-(1.0/(Math.pow(a,(relPos*2))));
 143             else
 144                 weightsX[x] = 1.0-(1.0/(Math.pow(a,((1-relPos)*2))));
 145         }
 146
 147         for (int y = 0; y < img.height; y++)
 148         {
 149             double relPos = (double)y / (double)(img.height-1);
 150
 151             if (relPos <= 0.5)
 152                 weightsY[y] = 1.0-(1.0/(Math.pow(a,(relPos*2))));
 153             else
 154                 weightsY[y] = 1.0-(1.0/(Math.pow(a,((1-relPos)*2))));
 155         }
 156
 157         for (int z = 0; z < img.depth; z++)
 158         {
 159             double relPos = (double)z / (double)(img.depth-1);
 160
 161             if (relPos <= 0.5)
 162                 weightsZ[z] = 1.0-(1.0/(Math.pow(a,(relPos*2))));
 163             else
 164                 weightsZ[z] = 1.0-(1.0/(Math.pow(a,((1-relPos)*2))));
 165         }
 166
 167
 168         for (int z = 0; z < img.depth; z++)
 169             for (int y = 0; y < img.height; y++)
 170                 for (int x = 0; x < img.width; x++)
 171                     img.set((float) (img.get(x, y, z) * weightsX[x] * weightsY[y]* weightsZ[z]), x, y, z);
 172     }
 173
 174     /**
 175      * This class creates a gaussian kernel
 176      *
 177      * @param sigma Standard Derivation of the gaussian function
 178      * @param normalize Normalize integral of gaussian function to 1 or not...
 179      * @return float[] The gaussian kernel
 180      *
 181      * @author   Stephan Saalfeld
 182      */
 183     public static float[] createGaussianKernel1D(float sigma, boolean normalize)
 184     {
 185         int size = 3;
 186         float[] gaussianKernel;
 187
 188         if (sigma <= 0)
 189         {
 190          gaussianKernel = new float[3];
 191          gaussianKernel[1] = 1;
 192         }
 193         else
 194         {
 195          size = max(3, (int)(2*(int)(3*sigma + 0.5)+1));
 196
 197          float two_sq_sigma = 2*sigma*sigma;
 198          gaussianKernel = new float[size];
 199
 200          for (int x = size/2; x >= 0; --x)
 201          {
 202              float val = (float)Math.exp(-(float)(x*x)/two_sq_sigma);
 203
 204              gaussianKernel[size/2-x] = val;
 205              gaussianKernel[size/2+x] = val;
 206          }
 207      }
 208
 209      if (normalize)
 210      {
 211          float sum = 0;
 212          for (float value : gaussianKernel)
 213              sum += value;
 214
 215          for (int i = 0; i < gaussianKernel.length; i++)
 216              gaussianKernel[i] /= sum;
 217      }
 218
 219
 220         return gaussianKernel;
 221     }
 222
 223     public static void normalize(FloatArray img)
 224     {
 225         // compute average
 226         /*float avg = 0;
 227         for (float value : img.data)
 228             avg += value;
 229
 230         avg /= (float)img.data.length;*/
 231
 232         float avg = computeTurkeyBiMedian(img.data, 5);
 233
 234         // compute stdev
 235         float stDev = 0;
 236         for (float value : img.data)
 237             stDev += (value - avg) * (value - avg);
 238
 239         stDev /= (float)(img.data.length - 1);
 240         stDev = (float)Math.sqrt(stDev);
 241
 242         float min = Float.MAX_VALUE;
 243         for (int i = 0; i < img.data.length; i++)
 244         {
 245             img.data[i] = (img.data[i] - avg) / stDev;
 246             if (img.data[i] < min)
 247                 min = img.data[i];
 248         }
 249
 250         for (int i = 0; i < img.data.length; i++)
 251             img.data[i] -= min;
 252
 253     }
 254
 255     public static FloatArray2D createGaussianKernel2D(float sigma, boolean normalize)
 256     {
 257         int size = 3;
 258         FloatArray2D gaussianKernel;
 259
 260         if (sigma <= 0)
 261         {
 262          gaussianKernel = new FloatArray2D(3, 3);
 263          gaussianKernel.data[4] = 1;
 264         }
 265         else
 266         {
 267          size = max(3, (int)(2*(int)(3*sigma + 0.5)+1));
 268
 269          float two_sq_sigma = 2*sigma*sigma;
 270          gaussianKernel = new FloatArray2D(size, size);
 271
 272          for (int y = size/2; y >= 0; --y)
 273          {
 274              for (int x = size/2; x >= 0; --x)
 275              {
 276               float val = (float)Math.exp(-(float)(y*y+x*x)/two_sq_sigma);
 277
 278               gaussianKernel.set(val, size/2-x, size/2-y);
 279               gaussianKernel.set(val, size/2-x, size/2+y);
 280               gaussianKernel.set(val, size/2+x, size/2-y);
 281               gaussianKernel.set(val, size/2+x, size/2+y);
 282              }
 283          }
 284         }
 285
 286         if (normalize)
 287         {
 288             float sum = 0;
 289             for (float value : gaussianKernel.data)
 290                 sum += value;
 291
 292             for (int i = 0; i < gaussianKernel.data.length; i++)
 293                 gaussianKernel.data[i] /= sum;
 294         }
 295
 296
 297         return gaussianKernel;
 298     }
 299
 300     /*
 301     ** create a normalized gaussian impulse with appropriate size and offset center
 302     */
 303     static public FloatArray2D create_gaussian_kernel_2D_offset(
 304                     float sigma,
 305                     float offset_x,
 306                     float offset_y,
 307                     boolean normalize)
 308     {
 309             int size = 3;
 310             FloatArray2D gaussian_kernel;
 311             if (sigma == 0)
 312             {
 313                     gaussian_kernel = new FloatArray2D(3 ,3);
 314                     gaussian_kernel.data[4] = 1;
 315             }
 316             else
 317             {
 318                     size = Math.max(3, (int)( 2 * Math.round( 3 * sigma ) + 1 ) );
 319                     float two_sq_sigma = 2*sigma*sigma;
 320                     // float normalization_factor = 1.0/(float)M_PI/two_sq_sigma;
 321                     gaussian_kernel = new FloatArray2D( size, size );
 322                     for ( int x = size - 1; x >= 0; --x )
 323                     {
 324                             float fx = (float)( x - size / 2 );
 325                             for ( int y = size-1; y >= 0; --y )
 326                             {
 327                                     float fy = (float)(y-size/2);
 328                                     float val = (float)( Math.exp( -( Math.pow( fx - offset_x, 2)+Math.pow(fy-offset_y, 2))/two_sq_sigma));
 329                                     gaussian_kernel.set(val, x, y);
 330                             }
 331                     }
 332             }
 333             if (normalize)
 334             {
 335                     float sum = 0;
 336                     for (float value : gaussian_kernel.data)
 337                     sum += value;
 338
 339                     for (int i = 0; i < gaussian_kernel.data.length; i++)
 340                     gaussian_kernel.data[i] /= sum;
 341             }
 342             return gaussian_kernel;
 343     }
 344
 345     public static FloatArray3D createGaussianKernel3D(float sigma, boolean normalize)
 346     {
 347         int size = 3;
 348         FloatArray3D gaussianKernel;
 349
 350         if (sigma <= 0)
 351         {
 352          gaussianKernel = new FloatArray3D(3, 3, 3);
 353          gaussianKernel.set(1f, 1, 1, 1);
 354         }
 355         else
 356         {
 357          size = max(3, (int)(2*(int)(3*sigma + 0.5)+1));
 358
 359          float two_sq_sigma = 2*sigma*sigma;
 360          gaussianKernel = new FloatArray3D(size, size, size);
 361
 362          for (int z = size/2; z >= 0; --z)
 363          {
 364              for (int y = size / 2; y >= 0; --y)
 365              {
 366                  for (int x = size / 2; x >= 0; --x)
 367                  {
 368                      float val = (float) Math.exp( -(float) (z * z + y * y + x * x) / two_sq_sigma);
 369
 370                      gaussianKernel.set(val, size / 2 - x, size / 2 - y, size / 2 - z);
 371
 372                      gaussianKernel.set(val, size / 2 - x, size / 2 - y, size / 2 + z);
 373                      gaussianKernel.set(val, size / 2 - x, size / 2 + y, size / 2 - z);
 374                      gaussianKernel.set(val, size / 2 + x, size / 2 - y, size / 2 - z);
 375
 376                      gaussianKernel.set(val, size / 2 + x, size / 2 + y, size / 2 - z);
 377                      gaussianKernel.set(val, size / 2 + x, size / 2 - y, size / 2 + z);
 378                      gaussianKernel.set(val, size / 2 - x, size / 2 + y, size / 2 + z);
 379
 380                      gaussianKernel.set(val, size / 2 + x, size / 2 + y, size / 2 + z);
 381                  }
 382              }
 383          }
 384         }
 385
 386         if (normalize)
 387         {
 388             float sum = 0;
 389             for (float value : gaussianKernel.data)
 390                 sum += value;
 391
 392             for (int i = 0; i < gaussianKernel.data.length; i++)
 393                 gaussianKernel.data[i] /= sum;
 394         }
 395
 396
 397         return gaussianKernel;
 398     }
 399
 400     public static FloatArray2D computeIncreasingGaussianX(FloatArray2D input, float stDevStart, float stDevEnd)
 401     {
 402         FloatArray2D output = new FloatArray2D(input.width, input.height);
 403
 404         int width = input.width;
 405         float changeFilterSize = (float)(stDevEnd - stDevStart)/(float)width;
 406         float sigma;
 407         int filterSize;
 408
 409         float avg;
 410
 411         for (int x = 0; x < input.width; x++)
 412         {
 413             sigma = stDevStart + changeFilterSize * (float) x;
 414             FloatArray2D kernel = createGaussianKernel2D(sigma, true);
 415             filterSize = kernel.width;
 416
 417             for (int y = 0; y < input.height; y++)
 418             {
 419                 avg = 0;
 420
 421                 for (int fx = -filterSize / 2; fx <= filterSize / 2; fx++)
 422                     for (int fy = -filterSize / 2; fy <= filterSize / 2; fy++)
 423                     {
 424                         try
 425                         {
 426                             avg += input.get(x + fx, y + fy) * kernel.get(fx + filterSize/2, fy + filterSize/2);
 427                         }
 428                         catch (Exception e)
 429                         {}
 430                         ;
 431
 432                     }
 433
 434                 output.set(avg, x, y);
 435             }
 436         }
 437         return output;
 438     }
 439
 440     public static FloatArray2D computeGaussian(FloatArray2D input, float sigma)
 441     {
 442         FloatArray2D output = new FloatArray2D(input.width, input.height);
 443
 444         float avg, kernelsum;
 445
 446         FloatArray2D kernel = createGaussianKernel2D(sigma, true);
 447         int filterSize = kernel.width;
 448
 449         for (int x = 0; x < input.width; x++)
 450         {
 451             for (int y = 0; y < input.height; y++)
 452             {
 453                 avg = 0;
 454                 kernelsum = 0;
 455
 456                 for (int fx = -filterSize / 2; fx <= filterSize / 2; fx++)
 457                     for (int fy = -filterSize / 2; fy <= filterSize / 2; fy++)
 458                     {
 459                         try
 460                         {
 461                             avg += input.get(x + fx, y + fy) * kernel.get(fx + filterSize/2, fy + filterSize/2);
 462                             kernelsum += kernel.get(fx + filterSize/2, fy + filterSize/2);
 463                         }
 464                         catch (Exception e)
 465                         {};
 466
 467                     }
 468
 469                 output.set(avg/kernelsum, x, y);
 470             }
 471         }
 472         return output;
 473     }
 474
 475     public static FloatArray3D computeGaussian(FloatArray3D input, float sigma)
 476     {
 477         FloatArray3D output = new FloatArray3D(input.width, input.height, input.depth);
 478
 479         float avg, kernelsum;
 480
 481         FloatArray3D kernel = createGaussianKernel3D(sigma, true);
 482         int filterSize = kernel.width;
 483
 484         for (int x = 0; x < input.width; x++)
 485         {
 486             for (int y = 0; y < input.height; y++)
 487             {
 488                 for (int z = 0; z < input.depth; z++)
 489                 {
 490                     avg = 0;
 491                     kernelsum = 0;
 492
 493                     for (int fx = -filterSize / 2; fx <= filterSize / 2; fx++)
 494                         for (int fy = -filterSize / 2; fy <= filterSize / 2; fy++)
 495                             for (int fz = -filterSize / 2; fz <= filterSize / 2; fz++)
 496                             {
 497                                 try
 498                                 {
 499                                     avg += input.get(x + fx, y + fy, z + fz) * kernel.get(fx + filterSize / 2, fy + filterSize / 2, fz + filterSize / 2);
 500                                     kernelsum += kernel.get(fx + filterSize / 2, fy + filterSize / 2, fz + filterSize / 2);
 501                                 } catch (Exception e)
 502                                 {}
 503                                 ;
 504
 505                             }
 506
 507                     output.set(avg / kernelsum, x, y, z);
 508                 }
 509             }
 510         }
 511         return output;
 512     }
 513
 514     public static FloatArray3D computeGaussianFast(FloatArray3D input, float sigma)
 515     {
 516         FloatArray3D output = new FloatArray3D(input.width, input.height, input.depth);
 517
 518         float avg, kernelsum;
 519
 520         float[] kernel = createGaussianKernel1D(sigma, true);
 521
 522         int filterSize = kernel.length;
 523
 524         // fold in x
 525         for (int x = 0; x < input.width; x++)
 526             for (int y = 0; y < input.height; y++)
 527                 for (int z = 0; z < input.depth; z++)
 528                 {
 529                     avg = 0;
 530                     kernelsum = 0;
 531
 532                     for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 533                     {
 534                         if (x+f >= 0 && x+f < input.width)
 535                         {
 536                             avg += input.get(x + f, y, z) * kernel[f + filterSize / 2];
 537                             kernelsum += kernel[f + filterSize / 2];
 538                         }
 539                     }
 540
 541                     output.set(avg / kernelsum, x, y, z);
 542
 543                 }
 544
 545         // fold in y
 546         for (int x = 0; x < input.width; x++)
 547             for (int z = 0; z < input.depth; z++)
 548             {
 549                 float[] temp = new float[input.height];
 550
 551                 for (int y = 0; y < input.height; y++)
 552                 {
 553                     avg = 0;
 554                     kernelsum = 0;
 555
 556                     for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 557                     {
 558                         if (y+f >= 0 && y+f < input.height)
 559                         {
 560                             avg += output.get(x, y + f, z) * kernel[f + filterSize / 2];
 561                             kernelsum += kernel[f + filterSize / 2];
 562                         }
 563                     }
 564                     temp[y] = avg / kernelsum;
 565                 }
 566
 567                 for (int y = 0; y < input.height; y++)
 568                     output.set(temp[y], x, y, z);
 569             }
 570
 571         // fold in z
 572         for (int x = 0; x < input.width; x++)
 573             for (int y = 0; y < input.height; y++)
 574             {
 575                 float[] temp = new float[input.depth];
 576
 577                 for (int z = 0; z < input.depth; z++)
 578                 {
 579                     avg = 0;
 580                     kernelsum = 0;
 581
 582                     for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 583                     {
 584                         if (z+f >= 0 && z+f < input.depth)
 585                         {
 586                             avg += output.get(x, y, z + f) * kernel[f + filterSize / 2];
 587                             kernelsum += kernel[f + filterSize / 2];
 588                         }
 589                     }
 590                     temp[z] = avg / kernelsum;
 591                 }
 592
 593                 for (int z = 0; z < input.depth; z++)
 594                     output.set(temp[z], x, y, z);
 595
 596             }
 597         return output;
 598     }
 599
 600     public static FloatArray2D computeGaussianFastMirror(final FloatArray2D input, final float sigma)
 601     {
 602         final FloatArray2D output = new FloatArray2D(input.width, input.height);
 603
 604         float avg, kernelsum = 0;
 605         final float[] kernel = createGaussianKernel1D(sigma, true);
 606         final int filterSize = kernel.length;
 607
 608         // get kernel sum
 609         for (double value : kernel)
 610             kernelsum += value;
 611
 612         // fold in x
 613         for (int x = 0; x < input.width; x++)
 614             for (int y = 0; y < input.height; y++)
 615                 {
 616                     avg = 0;
 617
 618                     if (x -filterSize / 2 >= 0 && x + filterSize / 2 < input.width)
 619                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 620                             avg += input.get(x + f, y) * kernel[f + filterSize / 2];
 621                     else
 622                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 623                             avg += input.getMirror(x + f, y) * kernel[f + filterSize / 2];
 624
 625                     output.set(avg / kernelsum, x, y);
 626
 627                 }
 628
 629         // fold in y
 630         final float[] temp = new float[input.height];
 631         for (int x = 0; x < input.width; x++)
 632             {
 633
 634                 for (int y = 0; y < input.height; y++)
 635                 {
 636                     avg = 0;
 637
 638                     if (y -filterSize / 2 >= 0 && y + filterSize / 2 < input.height)
 639                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 640                             avg += output.get(x, y + f) * kernel[f + filterSize / 2];
 641                      else
 642                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 643                             avg += output.getMirror(x, y + f) * kernel[f + filterSize / 2];
 644
 645                     temp[y] = avg / kernelsum;
 646                 }
 647
 648                 for (int y = 0; y < input.height; y++)
 649                     output.set(temp[y], x, y);
 650             }
 651
 652         return output;
 653     }
 654
 655     /**
 656      * This class does the gaussian filtering of an image. On the edges of
 657      * the image it does mirror the pixels. It also uses the seperability of
 658      * the gaussian convolution.
 659      *
 660      * @param input FloatProcessor which should be folded (will not be touched)
 661      * @param sigma Standard Derivation of the gaussian function
 662      * @return FloatProcessor The folded image
 663      *
 664      * @author   Stephan Preibisch
 665      */
 666
 667     public static FloatArray3D computeGaussianFastMirror(FloatArray3D input, float sigma)
 668     {
 669         FloatArray3D output = new FloatArray3D(input.width, input.height, input.depth);
 670
 671         float avg, kernelsum = 0;
 672         float[] kernel = createGaussianKernel1D(sigma, true);
 673         int filterSize = kernel.length;
 674
 675         // get kernel sum
 676         for (double value : kernel)
 677             kernelsum += value;
 678
 679         // fold in x
 680         for (int x = 0; x < input.width; x++)
 681             for (int y = 0; y < input.height; y++)
 682                 for (int z = 0; z < input.depth; z++)
 683                 {
 684                     avg = 0;
 685
 686                     if (x -filterSize / 2 >= 0 && x + filterSize / 2 < input.width)
 687                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 688                             avg += input.get(x + f, y, z) * kernel[f + filterSize / 2];
 689                     else
 690                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 691                             avg += input.getMirror(x + f, y, z) * kernel[f + filterSize / 2];
 692
 693                     output.set(avg / kernelsum, x, y, z);
 694
 695                 }
 696
 697         // fold in y
 698         for (int x = 0; x < input.width; x++)
 699             for (int z = 0; z < input.depth; z++)
 700             {
 701                 float[] temp = new float[input.height];
 702
 703                 for (int y = 0; y < input.height; y++)
 704                 {
 705                     avg = 0;
 706
 707                     if (y -filterSize / 2 >= 0 && y + filterSize / 2 < input.height)
 708                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 709                             avg += output.get(x, y + f, z) * kernel[f + filterSize / 2];
 710                     else
 711                        for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 712                            avg += output.getMirror(x, y + f, z) * kernel[f + filterSize / 2];
 713
 714                     temp[y] = avg / kernelsum;
 715                 }
 716
 717                 for (int y = 0; y < input.height; y++)
 718                     output.set(temp[y], x, y, z);
 719             }
 720
 721         // fold in z
 722         for (int x = 0; x < input.width; x++)
 723             for (int y = 0; y < input.height; y++)
 724             {
 725                 float[] temp = new float[input.depth];
 726
 727                 for (int z = 0; z < input.depth; z++)
 728                 {
 729                     avg = 0;
 730
 731                     if (z -filterSize / 2 >= 0 && z + filterSize / 2 < input.depth)
 732                         for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 733                             avg += output.get(x, y, z + f) * kernel[f + filterSize / 2];
 734                     else
 735                        for (int f = -filterSize / 2; f <= filterSize / 2; f++)
 736                            avg += output.getMirror(x, y, z + f) * kernel[f + filterSize / 2];
 737
 738                     temp[z] = avg / kernelsum;
 739                 }
 740
 741                 for (int z = 0; z < input.depth; z++)
 742                     output.set(temp[z], x, y, z);
 743
 744             }
 745         return output;
 746     }
 747
 748     public static FloatArray2D distortSamplingX(FloatArray2D input)
 749     {
 750         FloatArray2D output = new FloatArray2D(input.width, input.height);
 751
 752         int filterSize = 3;
 753         float avg;
 754
 755         Random rnd = new Random(353245632);
 756
 757         for (int x = 0; x < input.width; x++)
 758         {
 759             FloatArray2D kernel = new FloatArray2D(3,1);
 760
 761             float random = (rnd.nextFloat()-0.5f)*2;
 762             float val1, val2, val3;
 763
 764             if (random < 0)
 765             {
 766                 val1 = -random;
 767                 val2 = 1+random;
 768                 val3 = 0;
 769             }
 770             else
 771             {
 772                 val3 = random;
 773                 val2 = 1-random;
 774                 val1 = 0;
 775             }
 776
 777             kernel.set(val1, 0, 0);
 778             kernel.set(val2, 1, 0);
 779             kernel.set(val3, 2, 0);
 780
 781             for (int y = 0; y < input.height; y++)
 782             {
 783                 avg = 0;
 784
 785                 for (int fx = -filterSize / 2; fx <= filterSize / 2; fx++)
 786                 {
 787                     try
 788                     {
 789                         avg += input.get(x + fx, y) * kernel.get(fx + filterSize / 2, 0);
 790                     } catch (Exception e)
 791                     {}
 792                     ;
 793                 }
 794
 795                 output.set(avg, x, y);
 796             }
 797         }
 798         return output;
 799     }
 800
 801     public static FloatArray2D distortSamplingY(FloatArray2D input)
 802     {
 803         FloatArray2D output = new FloatArray2D(input.width, input.height);
 804
 805         int filterSize = 3;
 806         float avg;
 807
 808         Random rnd = new Random(7893469);
 809
 810         for (int y = 0; y < input.height; y++)
 811         {
 812             FloatArray2D kernel = new FloatArray2D(1,3);
 813
 814             float random = (rnd.nextFloat()-0.5f)*2;
 815             float val1, val2, val3;
 816
 817             if (random < 0)
 818             {
 819                 val1 = -random;
 820                 val2 = 1+random;
 821                 val3 = 0;
 822             }
 823             else
 824             {
 825                 val3 = random;
 826                 val2 = 1-random;
 827                 val1 = 0;
 828             }
 829
 830             kernel.set(val1, 0, 0);
 831             kernel.set(val2, 0, 1);
 832             kernel.set(val3, 0, 2);
 833
 834             for (int x = 0; x < input.width; x++)
 835             {
 836                 avg = 0;
 837
 838                 for (int fy = -filterSize / 2; fy <= filterSize / 2; fy++)
 839                 {
 840                     try
 841                     {
 842                         avg += input.get(x, y + fy) * kernel.get(0, fy + filterSize / 2);
 843                     } catch (Exception e)
 844                     {}
 845                     ;
 846                 }
 847
 848                 output.set(avg, x, y);
 849             }
 850         }
 851         return output;
 852     }
 853
 854     public static FloatArray3D computeSobelFilter(FloatArray3D input)
 855     {
 856         FloatArray3D output = new FloatArray3D(input.width, input.height, input.depth);
 857
 858         float sobelX, sobelY, sobelZ;
 859         float v1, v2, v3, v4, v5, v6, v7, v8, v9, v10;
 860
 861         for (int z = 1; z < input.depth -1 ; z++)
 862             for (int y = 1; y < input.height -1 ; y++)
 863                 for (int x = 1; x < input.width -1; x++)
 864                 {
 865                     // Sobel in X
 866                     v1 = input.get(x-1,y,z-1);
 867                     v2 = input.get(x+1,y,z-1);
 868
 869                     v3 = input.get(x-1,y-1,z);
 870                     v4 = input.get(x-1,y,  z);
 871                     v5 = input.get(x-1,y+1,z);
 872
 873                     v6 = input.get(x+1,y-1,z);
 874                     v7 = input.get(x+1,y,  z);
 875                     v8 = input.get(x+1,y+1,z);
 876
 877                     v9 = input.get(x-1,y,z+1);
 878                     v10 = input.get(x+1,y,z+1);
 879
 880                     sobelX = v1    + (-v2)   +
 881                              v3    + (2*v4)  + v5 +
 882                              (-v6) + (-2*v7) + (-v8) +
 883                              v9    + (-v10);
 884
 885                     // Sobel in Y
 886                     v1 = input.get(x,y-1,z-1);
 887                     v2 = input.get(x,y+1,z-1);
 888
 889                     v3 = input.get(x-1,y-1,z);
 890                     v4 = input.get(x,  y-1,z);
 891                     v5 = input.get(x+1,y-1,z);
 892
 893                     v6 = input.get(x-1,y+1,z);
 894                     v7 = input.get(x  ,y+1,z);
 895                     v8 = input.get(x+1,y+1,z);
 896
 897                     v9 =  input.get(x,y-1,z+1);
 898                     v10 = input.get(x,y+1,z+1);
 899
 900                     sobelY = v1    + (-v2)   +
 901                              v3    + (2*v4)  + v5 +
 902                              (-v6) + (-2*v7) + (-v8) +
 903                              v9    + (-v10);
 904
 905
 906                     // Sobel in Z
 907                     v1 = input.get(x  ,y-1,z-1);
 908                     v2 = input.get(x-1,y,  z-1);
 909                     v3 = input.get(x  ,y,  z-1);
 910                     v4 = input.get(x+1,y,  z-1);
 911                     v5 = input.get(x  ,y+1,z-1);
 912
 913                     v6 =  input.get(x  ,y-1,z+1);
 914                     v7 =  input.get(x-1,y,  z+1);
 915                     v8 =  input.get(x  ,y,  z+1);
 916                     v9 =  input.get(x+1,y,  z+1);
 917                     v10 = input.get(x  ,y+1,z+1);
 918
 919                     sobelZ =         v1      +
 920                              v2    + (2*v3)  + v4    +
 921                                      v5      +
 922                                      (-v6)   +
 923                              (-v7) + (-2*v8) + (-v9) +
 924                                      (-v10);
 925
 926
 927                     output.set((float)Math.sqrt(Math.pow(sobelX,2) + Math.pow(sobelY,2)  + Math.pow(sobelZ,2)), x, y, z);
 928                 }
 929
 930         return output;
 931     }
 932
 933     public static FloatArray3D computeBinaryFilter3(FloatArray3D input, float threshold)
 934     {
 935         FloatArray3D output = new FloatArray3D(input.width, input.height, input.depth);
 936
 937         for (int z = 1; z < input.depth - 1; z++)
 938             for (int y = 1; y < input.height - 1; y++)
 939                 for (int x = 1; x < input.width - 1; x++)
 940                 {
 941                     float avg = 0;
 942
 943                     for (int xs = x -1; xs <= x + 1; xs++)
 944                         for (int ys = y -1; ys <= y + 1; ys++)
 945                             for (int zs = z -1; zs <= z + 1; zs++)
 946                                 avg += input.get(xs, ys, zs);
 947
 948                     avg /= 27f;
 949
 950                     if (avg > threshold)
 951                         output.set(1, x, y, z);
 952                     else
 953                         output.set(0, x, y, z);
 954                 }
 955
 956         return output;
 957     }
 958
 959     public static FloatArray3D computeBinaryPlusSobelFilter3(FloatArray3D input, float threshold)
 960     {
 961         FloatArray3D output = computeSobelFilter(input);
 962
 963         // find maximum
 964         float max = 0;
 965         float min = Float.MAX_VALUE;
 966
 967         for (float value : output.data)
 968         {
 969             if (value > max)
 970                 max = value;
 971
 972             if (value < min)
 973                 min = value;
 974         }
 975
 976         // norm to 0...1
 977         for (int i = 0; i < output.data.length; i++)
 978             output.data[i] = (output.data[i]-min)/(max-min)*2;
 979
 980         // add
 981         for (int z = 1; z < input.depth - 1; z++)
 982             for (int y = 1; y < input.height - 1; y++)
 983                 for (int x = 1; x < input.width - 1; x++)
 984                 {
 985                     float avg = 0;
 986
 987                     for (int xs = x -1; xs <= x + 1; xs++)
 988                         for (int ys = y -1; ys <= y + 1; ys++)
 989                             for (int zs = z -1; zs <= z + 1; zs++)
 990                                 avg += input.get(xs, ys, zs);
 991
 992                     avg /= 27f;
 993
 994                     if (avg > threshold)
 995                         output.set(output.get(x,y,z) + 1, x, y, z);
 996                 }
 997
 998         return output;
 999     }
1000
1001     public static FloatArray2D computeLaPlaceFilter3(FloatArray2D input)
1002     {
1003         FloatArray2D output = new FloatArray2D(input.width, input.height);
1004
1005         float derivX, derivY;
1006         float x1, x2, x3;
1007         float y1, y2, y3;
1008
1009         for (int y = 1; y < input.height -1 ; y++)
1010             for (int x = 1; x < input.width -1; x++)
1011             {
1012                 x1 = input.get(x-1,y);
1013                 x2 = input.get(x,y);
1014                 x3 = input.get(x+1,y);
1015
1016                 derivX = x1 - 2*x2 + x3;
1017
1018                 y1 = input.get(x,y-1);
1019                 y2 = input.get(x,y);
1020                 y3 = input.get(x,y+1);
1021
1022                 derivY = y1 - 2*y2 + y3;
1023
1024                 output.set((float)Math.sqrt(Math.pow(derivX,2) + Math.pow(derivY,2)), x, y);
1025             }
1026
1027         return output;
1028     }
1029
1030     /*public static void computeLaPlaceFilterInPlace3(FloatArray2D input)
1031     {
1032
1033         float buffer = max(input.height, input.width);
1034
1035         float derivX, derivY;
1036         float x1, x2, x3;
1037         float y1, y2, y3;
1038
1039         for (int diag = 1; diag < input.width + input.height - 1; diag++)
1040         {
1041
1042         }
1043
1044         for (int y = 1; y < input.height -1 ; y++)
1045             for (int x = 1; x < input.width -1; x++)
1046             {
1047                 x1 = input.get(x-1,y);
1048                 x2 = input.get(x,y);
1049                 x3 = input.get(x+1,y);
1050
1051                 derivX = x1 - 2*x2 + x3;
1052
1053                 y1 = input.get(x,y-1);
1054                 y2 = input.get(x,y);
1055                 y3 = input.get(x,y+1);
1056
1057                 derivY = y1 - 2*y2 + y3;
1058
1059                 //output.set((float)Math.sqrt(Math.pow(derivX,2) + Math.pow(derivY,2)), x, y);
1060             }
1061
1062         return;
1063     }*/
1064
1065     public static FloatArray2D computeLaPlaceFilter5(FloatArray2D input)
1066     {
1067         FloatArray2D output = new FloatArray2D(input.width, input.height);
1068
1069         float derivX, derivY;
1070         float x1, x3, x5;
1071         float y1, y3, y5;
1072
1073         for (int y = 2; y < input.height -2 ; y++)
1074             for (int x = 2; x < input.width -2; x++)
1075             {
1076                 x1 = input.get(x-2,y);
1077                 x3 = input.get(x,y);
1078                 x5 = input.get(x+2,y);
1079
1080                 derivX = x1 - 2*x3 + x5;
1081
1082                 y1 = input.get(x,y-2);
1083                 y3 = input.get(x,y);
1084                 y5 = input.get(x,y+2);
1085
1086                 derivY = y1 - 2*y3 + y5;
1087
1088                 output.set((float)Math.sqrt(Math.pow(derivX,2) + Math.pow(derivY,2)), x, y);
1089             }
1090
1091         return output;
1092     }
1093
1094     public static FloatArray3D computeLaPlaceFilter5(FloatArray3D input)
1095     {
1096         FloatArray3D output = new FloatArray3D(input.width, input.height, input.depth);
1097
1098         float derivX, derivY, derivZ;
1099         float x1, x3, x5;
1100         float y1, y3, y5;
1101         float z1, z3, z5;
1102
1103         for (int z = 2; z < input.depth -2 ; z++)
1104             for (int y = 2; y < input.height -2 ; y++)
1105                 for (int x = 2; x < input.width -2; x++)
1106                 {
1107                     x1 = input.get(x-2,y,z);
1108                     x3 = input.get(x,y,z);
1109                     x5 = input.get(x+2,y,z);
1110
1111                     derivX = x1 - 2*x3 + x5;
1112
1113                     y1 = input.get(x,y-2,z);
1114                     y3 = input.get(x,y,z);
1115                     y5 = input.get(x,y+2,z);
1116
1117                     derivY = y1 - 2*y3 + y5;
1118
1119                     z1 = input.get(x,y,z-2);
1120                     z3 = input.get(x,y,z);
1121                     z5 = input.get(x,y,z+2);
1122
1123                     derivZ = z1 - 2*z3 + z5;
1124
1125                     output.set((float)Math.sqrt(Math.pow(derivX,2) + Math.pow(derivY,2)  + Math.pow(derivZ,2)), x, y, z);
1126                 }
1127
1128         return output;
1129     }
1130
1131     public static FloatArray2D[] createGradients( FloatArray2D array)
1132      {
1133          FloatArray2D[] gradients = new FloatArray2D[2];
1134          gradients[0] = new FloatArray2D(array.width, array.height);
1135          gradients[1] = new FloatArray2D(array.width, array.height);
1136
1137          for (int y = 0; y < array.height; ++y)
1138          {
1139                  int[] ro = new int[3];
1140                      ro[0] = array.width * Math.max(0, y - 1);
1141                      ro[1] = array.width * y;
1142                      ro[2] = array.width * Math.min(y + 1, array.height - 1);
1143                  for (int x = 0; x < array.width; ++x)
1144                  {
1145                          // L(x+1, y) - L(x-1, y)
1146                          float der_x = (
1147                                          array.data[ro[1] + Math.min(x + 1, array.width - 1)] -
1148                                          array.data[ro[1] + Math.max(0, x - 1)]) / 2;
1149
1150                          // L(x, y+1) - L(x, y-1)
1151                          float der_y = (
1152                                  array.data[ro[2] + x] -
1153                                  array.data[ro[0] + x]) / 2;
1154
1155                          //! amplitude
1156                          gradients[0].data[ro[1]+x] = (float)Math.sqrt( Math.pow( der_x, 2 ) + Math.pow( der_y, 2 ) );
1157                          //! orientation
1158                          gradients[1].data[ro[1]+x] = (float)Math.atan2( der_y, der_x );
1159                  }
1160          }
1161          //ImageArrayConverter.FloatArrayToImagePlus( gradients[ 1 ], "gradients", 0, 0 ).show();
1162          return gradients;
1163      }
1164
1165     /**
1166      * in place enhance all values of a FloatArray to fill the given range
1167      */
1168     public static final void enhance( FloatArray2D src, float scale )
1169     {
1170         float min = Float.MAX_VALUE;
1171         float max = Float.MIN_VALUE;
1172         for ( int i = 0; i < src.data.length; ++i )
1173         {
1174                 if ( src.data[ i ] < min ) min = src.data[ i ];
1175                 else if ( src.data[ i ] > max ) max = src.data[ i ];
1176         }
1177         scale /= ( max - min );
1178         for ( int i = 0; i < src.data.length; ++i )
1179                 src.data[ i ] = scale * ( src.data[ i ] - min );
1180     }
1181
1182     /**
1183          * convolve an image with a horizontal and a vertical kernel
1184          * simple straightforward, not optimized---replace this with a trusted better version soon
1185          *
1186          * @param input the input image
1187          * @param h horizontal kernel
1188          * @param v vertical kernel
1189          *
1190          * @return convolved image
1191          */
1192         public static FloatArray2D convolveSeparable( FloatArray2D input, float[] h, float[] v )
1193         {
1194                 FloatArray2D output = new FloatArray2D( input.width, input.height );
1195                 FloatArray2D temp = new FloatArray2D( input.width, input.height );
1196
1197                 int hl = h.length / 2;
1198                 int vl = v.length / 2;
1199
1200                 int xl = input.width - h.length + 1;
1201                 int yl = input.height - v.length + 1;
1202
1203                 // create lookup tables for coordinates outside the image range
1204                 int[] xb = new int[ h.length + hl - 1 ];
1205                 int[] xa = new int[ h.length + hl - 1 ];
1206                 for ( int i = 0; i < xb.length; ++i )
1207                 {
1208                         xb[ i ] = flipInRange( i - hl, input.width );
1209                         xa[ i ] = flipInRange( i + xl, input.width );
1210                 }
1211
1212                 int[] yb = new int[ v.length + vl - 1 ];
1213                 int[] ya = new int[ v.length + vl - 1 ];
1214                 for ( int i = 0; i < yb.length; ++i )
1215                 {
1216                         yb[ i ] = input.width * flipInRange( i - vl, input.height );
1217                         ya[ i ] = input.width * flipInRange( i + yl, input.height );
1218                 }
1219
1220 //              String xa_str = "xa: ";
1221 //              String xb_str = "xb: ";
1222 //              String ya_str = "ya: ";
1223 //              String yb_str = "yb: ";
1224 //              for ( int i = 0; i < xa.length; ++i )
1225 //              {
1226 //                      xa_str = xa_str + xa[ i ] + ", ";
1227 //                      xb_str = xb_str + xb[ i ] + ", ";
1228 //                      ya_str = ya_str + ( ya[ i ] / input.width ) + ", ";
1229 //                      yb_str = yb_str + ( yb[ i ] / input.width ) + ", ";
1230 //              }
1231 //
1232 //              System.out.println( xb_str );
1233 //              System.out.println( xa_str );
1234 //              System.out.println( yb_str );
1235 //              System.out.println( ya_str );
1236
1237
1238                 xl += hl;
1239                 yl += vl;
1240                 // horizontal convolution per row
1241                 int rl = input.height * input.width;
1242                 for ( int r = 0; r < rl; r += input.width )
1243                 {
1244                         for ( int x = hl; x < xl; ++x )
1245                         {
1246                                 int c = x - hl;
1247                                 float val = 0;
1248                                 for ( int xk = 0; xk < h.length; ++xk )
1249                                 {
1250                                         val += h[ xk ] * input.data[ r + c + xk ];
1251                                 }
1252                                 temp.data[ r + x ] = val;
1253                         }
1254                         for ( int x = 0; x < hl; ++x )
1255                         {
1256                                 float valb = 0;
1257                                 float vala = 0;
1258                                 for ( int xk = 0; xk < h.length; ++xk )
1259                                 {
1260                                         valb += h[ xk ] * input.data[ r + xb[ x + xk ] ];
1261                                         vala += h[ xk ] * input.data[ r + xa[ x + xk ] ];
1262                                 }
1263                                 temp.data[ r + x ] = valb;
1264                                 temp.data[ r + x + xl ] = vala;
1265                         }
1266                 }
1267
1268                 // vertical convolution per column
1269                 rl = yl * temp.width;
1270                 int vlc = vl * temp.width;
1271                 for ( int x = 0; x < temp.width; ++x )
1272                 {
1273                         for ( int r = vlc; r < rl; r += temp.width )
1274                         {
1275                                 float val = 0;
1276                                 int c = r - vlc;
1277                                 int rk = 0;
1278                                 for ( int yk = 0; yk < v.length; ++yk )
1279                                 {
1280                                         val += v[ yk ] * temp.data[ c + rk + x ];
1281                                         rk += temp.width;
1282                                 }
1283                                 output.data[ r + x ] = val;
1284                         }
1285                         for ( int y = 0; y < vl; ++y )
1286                         {
1287                                 int r = y * temp.width;
1288                                 float valb = 0;
1289                                 float vala = 0;
1290                                 for ( int yk = 0; yk < v.length; ++yk )
1291                                 {
1292                                         valb += h[ yk ] * temp.data[ yb[ y + yk ] + x ];
1293                                         vala += h[ yk ] * temp.data[ ya[ y + yk ] + x ];
1294                                 }
1295                                 output.data[ r + x ] = valb;
1296                                 output.data[ r + rl + x ] = vala;
1297                         }
1298                 }
1299
1300                 return output;
1301         }
1302
1303 }