TrakEM2_/src/main/java/mpi/fruitfly/math/General.java

   1 package mpi.fruitfly.math;
   2
   3 /**
   4  * <p>Title: </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 java.awt.Point;
  32 import java.util.Vector;
  33
  34 public class General
  35 {
  36         final public static double M_2PI = 2.0 * Math.PI;
  37
  38         /**
  39          * simple square
  40          * @param value
  41          * @return value * value
  42          */
  43         public static double sq( double a )
  44         {
  45                 return a * a;
  46         }
  47
  48         public static float sq( float a )
  49         {
  50                 return a * a;
  51         }
  52
  53         /**
  54          * fast floor ld of unsigned v
  55          */
  56         public static int ldu( int v )
  57         {
  58             int c = 0; // c will be lg( v )
  59             do
  60             {
  61                 v >>= 1;
  62                 c++;
  63             }
  64             while ( v > 1 );
  65             return c;
  66         }
  67
  68         /**
  69          * return a integer that is flipped in the range [0 ... mod - 1]
  70          *
  71          * @param a the value to be flipped
  72          * @param range the size of the range
  73          * @return a flipped in range like a ping pong ball
  74          */
  75         public static int flipInRange( int a, int mod )
  76         {
  77                 int p = 2 * mod;
  78                 if ( a < 0 ) a = p + a % p;
  79                 if ( a >= p ) a = a % p;
  80                 if ( a >= mod ) a = mod - a % mod - 1;
  81                 return a;
  82         }
  83
  84     public static int sign (double value)
  85     {
  86         if (value < 0)
  87             return -1;
  88         else
  89             return 1;
  90     }
  91
  92     public static boolean isApproxEqual(double a, double b, double threshold)
  93     {
  94       if (a==b)
  95         return true;
  96
  97       if (a + threshold > b && a - threshold < b)
  98         return true;
  99       else
 100         return false;
 101     }
 102
 103     public static double computeCorrelationCoefficient(double[] X1, double[] X2)
 104     {
 105         // check for valid input arrays
 106         if (X1 == null || X2 == null)
 107             return Double.NaN;
 108
 109         if (X1.length != X2.length)
 110             return Double.NaN;
 111
 112         int length = X1.length;
 113
 114         if (length == 0)
 115             return Double.NaN;
 116
 117         // compute averages
 118         double avgX1 = 0, avgX2 = 0;
 119
 120         for (int pos = 0; pos < length; pos++)
 121         {
 122             avgX1 += X1[pos];
 123             avgX2 += X2[pos];
 124         }
 125
 126         avgX1 /= (double)length;
 127         avgX2 /= (double)length;
 128
 129         // compute covariance and variances of X1 and X2
 130         double covar =  0, varX1 = 0, varX2 = 0;
 131         double distX1, distX2;
 132
 133         for (int pos = 0; pos < length; pos++)
 134         {
 135             distX1 = X1[pos] - avgX1;
 136             distX2 = X2[pos] - avgX2;
 137
 138             covar += distX1 * distX2;
 139
 140             varX1 += distX1 * distX1;
 141             varX2 += distX2 * distX2;
 142         }
 143
 144         covar /= (double)length;
 145         varX1 /= (double)length;
 146         varX2 /= (double)length;
 147
 148         varX1 = Math.sqrt(varX1);
 149         varX2 = Math.sqrt(varX2);
 150
 151         double R = covar / (varX1 * varX2);
 152
 153         return R;
 154     }
 155
 156     public static int min(int a, int b)
 157     {
 158         if (a > b)
 159             return b;
 160         else
 161             return a;
 162     }
 163
 164
 165     public static double min(double a, double b)
 166     {
 167         if (a > b)
 168             return b;
 169         else
 170             return a;
 171     }
 172
 173     public static float min(float a, float b)
 174     {
 175         if (a > b)
 176             return b;
 177         else
 178             return a;
 179     }
 180
 181     public static int max(int a, int b)
 182     {
 183         if (a > b)
 184             return a;
 185         else
 186             return b;
 187     }
 188
 189     public static float max(float a, float b)
 190     {
 191         if (a > b)
 192             return a;
 193         else
 194             return b;
 195     }
 196     public static double max(double a, double b)
 197     {
 198         if (a > b)
 199             return a;
 200         else
 201             return b;
 202     }
 203
 204     /**
 205      * Berechnet die Differenz(Offset) der ganzzahligen Eingabe zur naechsten Zweierpotenz.
 206      * Bsp: Eingabe 23 -> naechste Potenz 32 -> Ergebnis 9
 207      * @param size ganzzahliger Eingabewert
 208      * @return int der berechnete Offset
 209      */
 210     public static int offsetToPowerOf2(int size)
 211     {
 212         int offset = 0;
 213         int power = 2;
 214
 215         while (power < size)
 216         {
 217             power = power * 2;
 218             if (power == size)
 219             {
 220                 return 0;
 221             }
 222             offset = power - size;
 223         }
 224         return offset;
 225     }
 226
 227     public static double gLog(double z, double c)
 228     {
 229         if (c == 0)
 230             return z;
 231         else
 232             return Log10((z + Math.sqrt(z * z + c * c)) / 2.0);
 233     }
 234
 235     public static double gLogInv(double w, double c)
 236     {
 237         if (c == 0)
 238             return w;
 239         else
 240             return Exp10(w) - (((c * c) * Exp10( -w)) / 4.0);
 241     }
 242
 243     public static double Log10(double x)
 244     {
 245         if (x > 0)
 246             return (Math.log(x) / Math.log(10));
 247         else
 248             return 0;
 249     }
 250
 251     public static double Exp10(double x)
 252     {
 253         return (Math.pow(10, x));
 254     }
 255
 256     public static int pow(int a, int b)
 257     {
 258         if (b == 0)
 259             return 1;
 260
 261         if (b == 1)
 262             return a;
 263
 264         int result = a;
 265
 266         for (int i = 1; i < b; i++)
 267             result *= a;
 268
 269         return result;
 270     }
 271
 272     public static double computeTurkeyBiAverage(double[] conc)
 273     {
 274         return computeTurkeyBiAverage(conc, null, 5);
 275     }
 276
 277     public static double computeTurkeyBiAverage(double[] conc, double c)
 278     {
 279         return computeTurkeyBiAverage(conc, null, c);
 280     }
 281
 282     public static double computeTurkeyBiAverage(double[] conc, double[] W, double c)
 283     {
 284         // weighted or non-weighted
 285         if (W == null)
 286         {
 287             W = new double[conc.length];
 288             for (int i = 0; i < W.length; i++)
 289                 W[i] = 1;
 290         }
 291
 292         // compute average
 293         double avgConc = 0;
 294         double sumW = 0;
 295
 296         for (int i = 0; i < conc.length; i++)
 297         {
 298             avgConc += W[i] * conc[i];
 299             sumW += W[i];
 300         }
 301
 302         avgConc /= sumW;
 303
 304         // compute standard deviation
 305         double SD = 0;
 306
 307         for (int i = 0; i < conc.length; i++)
 308         {
 309             SD += W[i] * Math.pow(conc[i] - avgConc, 2);
 310         }
 311
 312         SD = Math.sqrt(SD / sumW);
 313
 314         // compute outlier
 315         double u[] = new double[conc.length];
 316         double e;
 317
 318         for (int i = 0; i < u.length; i++)
 319         {
 320             e = Math.abs(0.01 * conc[i]); // prevent division by zero
 321             u[i] = (W[i] * (conc[i] - avgConc)) / ((c * SD) + e);
 322         }
 323
 324         // compute turkey-bi weightening values
 325
 326         double tb[] = new double[conc.length];
 327
 328         for (int i = 0; i < conc.length; i++)
 329         {
 330             if (Math.abs(u[i]) > 1)
 331                 tb[i] = 0;
 332             else
 333                 tb[i] = Math.pow(1 - (u[i] * u[i]), 2);
 334         }
 335
 336         // compute result
 337         double concTB = 0;
 338         double temp = 0;
 339
 340         for (int i = 0; i < conc.length; i++)
 341         {
 342             concTB += tb[i] * conc[i];
 343             temp += tb[i];
 344         }
 345
 346         return (concTB / temp);
 347     }
 348
 349     public static double computeTurkeyBiMedian(double[] concinput, double c)
 350     {
 351         double[] conc = concinput.clone();
 352
 353         // compute median
 354         double medianConc = computeMedian(conc);
 355
 356         // compute "standard deviation"
 357         double SD = 0;
 358         double SDValues[] = new double[conc.length];
 359
 360         for (int i = 0; i < conc.length; i++)
 361         {
 362             SDValues[i] = Math.abs(conc[i] - medianConc);
 363         }
 364
 365         SD = computeMedian(SDValues);
 366
 367         // compute "ausreisser"
 368         double u[] = new double[conc.length];
 369         double e;
 370
 371         for (int i = 0; i < u.length; i++)
 372         {
 373             e = Math.abs(0.01 * conc[i]); // prevent division by zero
 374             u[i] = ((conc[i] - medianConc)) / ((c * SD) + e);
 375         }
 376
 377         // compute turkey-bi weightening values
 378
 379         double tb[] = new double[conc.length];
 380
 381         for (int i = 0; i < conc.length; i++)
 382         {
 383             if (Math.abs(u[i]) > 1)
 384                 tb[i] = 0;
 385             else
 386                 tb[i] = Math.pow(1 - (u[i] * u[i]), 2);
 387         }
 388
 389         // compute result
 390         double concTB = 0;
 391         double temp = 0;
 392
 393         for (int i = 0; i < conc.length; i++)
 394         {
 395             concTB += tb[i] * conc[i];
 396             temp += tb[i];
 397         }
 398
 399         return (concTB / temp);
 400     }
 401
 402     public static float computeTurkeyBiMedian(float[] concinput, float c)
 403     {
 404         float[] conc = concinput.clone();
 405
 406         // compute median
 407         float medianConc = computeMedian(conc);
 408
 409         // compute "standard deviation"
 410         float SD = 0;
 411         float SDValues[] = new float[conc.length];
 412
 413         for (int i = 0; i < conc.length; i++)
 414         {
 415             SDValues[i] = Math.abs(conc[i] - medianConc);
 416         }
 417
 418         SD = computeMedian(SDValues);
 419
 420         // compute "ausreisser"
 421         float u[] = new float[conc.length];
 422         float e;
 423
 424         for (int i = 0; i < u.length; i++)
 425         {
 426             e = (float)Math.abs(0.01 * conc[i]); // prevent division by zero
 427             u[i] = ((conc[i] - medianConc)) / ((c * SD) + e);
 428         }
 429
 430         // compute turkey-bi weightening values
 431
 432         float tb[] = new float[conc.length];
 433
 434         for (int i = 0; i < conc.length; i++)
 435         {
 436             if (Math.abs(u[i]) > 1)
 437                 tb[i] = 0;
 438             else
 439                 tb[i] = (float)Math.pow(1 - (u[i] * u[i]), 2);
 440         }
 441
 442         // compute result
 443         float concTB = 0;
 444         float temp = 0;
 445
 446         for (int i = 0; i < conc.length; i++)
 447         {
 448             concTB += tb[i] * conc[i];
 449             temp += tb[i];
 450         }
 451
 452         return (concTB / temp);
 453     }
 454
 455     public static double[] computeTurkeyBiMedianWeights(double[] concinput, double c)
 456     {
 457         double[] conc = concinput.clone();
 458
 459         // compute median
 460         double medianConc = computeMedian(conc);
 461
 462         // compute "standard deviation"
 463         double SD = 0;
 464         double SDValues[] = new double[conc.length];
 465
 466         for (int i = 0; i < conc.length; i++)
 467         {
 468             SDValues[i] = Math.abs(conc[i] - medianConc);
 469         }
 470
 471         SD = computeMedian(SDValues);
 472
 473         // compute "ausreisser"
 474         double u[] = new double[conc.length];
 475         double e;
 476
 477         for (int i = 0; i < u.length; i++)
 478         {
 479             e = Math.abs(0.01 * conc[i]); // prevent division by zero
 480             u[i] = ((conc[i] - medianConc)) / ((c * SD) + e);
 481         }
 482
 483         // compute turkey-bi weightening values
 484
 485         double tb[] = new double[conc.length];
 486
 487         for (int i = 0; i < conc.length; i++)
 488         {
 489             if (Math.abs(u[i]) > 1)
 490                 tb[i] = 0;
 491             else
 492                 tb[i] = Math.pow(1 - (u[i] * u[i]), 2);
 493         }
 494
 495         return tb;
 496     }
 497
 498     public static double computeMedian(double[] values)
 499     {
 500         double temp[] = values.clone();
 501         double median;
 502
 503         int length = temp.length;
 504
 505         quicksort(temp, 0, length - 1);
 506
 507         if (length % 2 == 1) //odd length
 508             median = temp[length / 2];
 509         else //even length
 510             median = (temp[length / 2] + temp[(length / 2) - 1]) / 2;
 511
 512         temp = null;
 513         return median;
 514     }
 515
 516     public static float computeMedian(float[] values)
 517     {
 518         float temp[] = values.clone();
 519         float median;
 520
 521         int length = temp.length;
 522
 523         quicksort(temp, 0, length - 1);
 524
 525         if (length % 2 == 1) //odd length
 526             median = temp[length / 2];
 527         else //even length
 528             median = (temp[length / 2] + temp[(length / 2) - 1]) / 2;
 529
 530         temp = null;
 531         return median;
 532     }
 533
 534     // entropie integral p * log(p)
 535     // entropie der joint distr. - entropie 1 - entropie 2
 536     // wieviel von der gemeinsamen entropie st
 537
 538     // helligkeiten sind ausgang des experimentes
 539
 540
 541     public static void quicksort(double[] data, int left, int right)
 542     {
 543         if (data == null || data.length < 2)return;
 544         int i = left, j = right;
 545         double x = data[(left + right) / 2];
 546         do
 547         {
 548             while (data[i] < x) i++;
 549             while (x < data[j]) j--;
 550             if (i <= j)
 551             {
 552                 double temp = data[i];
 553                 data[i] = data[j];
 554                 data[j] = temp;
 555                 i++;
 556                 j--;
 557             }
 558         }
 559         while (i <= j);
 560         if (left < j) quicksort(data, left, j);
 561         if (i < right) quicksort(data, i, right);
 562     }
 563
 564     public static void quicksort(float[] data, int left, int right)
 565     {
 566         if (data == null || data.length < 2)return;
 567         int i = left, j = right;
 568         float x = data[(left + right) / 2];
 569         do
 570         {
 571             while (data[i] < x) i++;
 572             while (x < data[j]) j--;
 573             if (i <= j)
 574             {
 575                 float temp = data[i];
 576                 data[i] = data[j];
 577                 data[j] = temp;
 578                 i++;
 579                 j--;
 580             }
 581         }
 582         while (i <= j);
 583         if (left < j) quicksort(data, left, j);
 584         if (i < right) quicksort(data, i, right);
 585     }
 586
 587     public static void quicksort(Quicksortable[] data, int left, int right)
 588     {
 589         if (data == null || data.length < 2)return;
 590         int i = left, j = right;
 591
 592         double x = data[(left + right) / 2].getQuicksortValue();
 593
 594         do
 595         {
 596             while (data[i].getQuicksortValue() < x) i++;
 597             while (x < data[j].getQuicksortValue()) j--;
 598             if (i <= j)
 599             {
 600                 Quicksortable temp = data[i];
 601                 data[i] = data[j];
 602                 data[j] = temp;
 603                 i++;
 604                 j--;
 605             }
 606         }
 607         while (i <= j);
 608         if (left < j) quicksort(data, left, j);
 609         if (i < right) quicksort(data, i, right);
 610     }
 611
 612     public static void quicksort(Vector data, int left, int right)
 613     {
 614         if (data == null || data.size() < 2)return;
 615         int i = left, j = right;
 616
 617         double x = ((Quicksortable)data.get((left + right) / 2)).getQuicksortValue();
 618         do
 619         {
 620             while (((Quicksortable)data.get(i)).getQuicksortValue() < x) i++;
 621             while (x < ((Quicksortable)data.get(j)).getQuicksortValue()) j--;
 622             if (i <= j)
 623             {
 624                 Quicksortable temp = (Quicksortable)data.get(i);
 625                 data.set( i, data.get( j ) );
 626                 data.set(j, temp);
 627                 i++;
 628                 j--;
 629             }
 630         }
 631         while (i <= j);
 632         if (left < j) quicksort(data, left, j);
 633         if (i < right) quicksort(data, i, right);
 634     }
 635
 636     public static void quicksort(float[] data, Point[] sortAlso, int left, int right)
 637     {
 638         if (data == null || data.length < 2)return;
 639         int i = left, j = right;
 640         float x = data[(left + right) / 2];
 641         do
 642         {
 643             while (data[i] < x) i++;
 644             while (x < data[j]) j--;
 645             if (i <= j)
 646             {
 647                 float temp = data[i];
 648                 data[i] = data[j];
 649                 data[j] = temp;
 650
 651                 Point temp2 = sortAlso[i];
 652                 sortAlso[i] = sortAlso[j];
 653                 sortAlso[j] = temp2;
 654
 655                 i++;
 656                 j--;
 657             }
 658         }
 659         while (i <= j);
 660         if (left < j) quicksort(data, sortAlso, left, j);
 661         if (i < right) quicksort(data, sortAlso, i, right);
 662     }
 663
 664     public static void quicksort(double[] data, int[] sortAlso, int left, int right)
 665     {
 666         if (data == null || data.length < 2)return;
 667         int i = left, j = right;
 668         double x = data[(left + right) / 2];
 669         do
 670         {
 671             while (data[i] < x) i++;
 672             while (x < data[j]) j--;
 673             if (i <= j)
 674             {
 675                 double temp = data[i];
 676                 data[i] = data[j];
 677                 data[j] = temp;
 678
 679                 int temp2 = sortAlso[i];
 680                 sortAlso[i] = sortAlso[j];
 681                 sortAlso[j] = temp2;
 682
 683                 i++;
 684                 j--;
 685             }
 686         }
 687         while (i <= j);
 688         if (left < j) quicksort(data, sortAlso, left, j);
 689         if (i < right) quicksort(data, sortAlso, i, right);
 690     }
 691
 692 }