IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31629


Ignore:
Timestamp:
Jun 15, 2011, 2:42:41 PM (15 years ago)
Author:
watersc1
Message:

Penultimate edit on new stacking code. Want to remove debug statements and do a final cleanup before claiming we've solved the bright artifact bug.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20110406/psModules/src/imcombine/pmStack.c

    r31618 r31629  
    3939# if (1)
    4040#define TESTING                         // Enable test output
    41 #define TEST_X 5745                       // x coordinate to examine
    42 #define TEST_Y 5331                       // y coordinate to examine
     41/* #define TEST_X 5745                       // x coordinate to examine */
     42/* #define TEST_Y 5331                       // y coordinate to examine */
     43#define TEST_X 25
     44#define TEST_Y 25
    4345#define TEST_RADIUS 0.5                 // Radius to examine
    4446# endif
    … …  
    102104    return;
    103105}
     106
     107// KMM functions to do bimodality rejection of pixels
     108
     109float gaussian(float x, float m, float s) {
     110  return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
     111}
     112
     113static void KMMcalculate(const psVector *values,
     114                         float *Punimodal,int *iter,
     115                         float *mU, float *sU,
     116                         float *pi1, float *m1, float *s1,
     117                         float *pi2, float *m2, float *s2) {
     118  double logL_bimodal = 0, logL_unimodal;
     119  psVector *P1 = psVectorAlloc(values->n,PS_TYPE_F32);
     120  psVector *P2 = psVectorAlloc(values->n,PS_TYPE_F32);
     121  int i;
     122/*   int debug = 0; */
     123/*   if (fabs(values->data.F32[0] - -145.624954) < 1e-4) { */
     124/*     debug = 1; */
     125/*   } */
     126 
     127  // Calculate unimodal properties
     128  *mU = 0;
     129  *sU = 0;
     130  logL_unimodal = 0;
     131  for (i = 0; i < values->n; i++) { // Calculate mean
     132    *mU += values->data.F32[i];
     133  }
     134  *mU /= values->n;
     135  for (i = 0; i < values->n; i++) { // Calculate sigma
     136    *sU += pow(values->data.F32[i] - *mU,2);
     137  }
     138  *sU = sqrt(*sU / values->n);
     139  for (i = 0; i < values->n; i++) { // Calculate log likelihood
     140    logL_unimodal += log(gaussian(values->data.F32[i],*mU,*sU));
     141  }
     142
     143  // Do EM loop
     144  float dL = 0;
     145  float oldL = -999;
     146  *iter = 0;
     147  logL_bimodal = logL_unimodal;
     148  *m1 = *mU - 3 * *sU;
     149  *m2 = *mU + 3 * *sU;
     150  *s1 = *sU / 2;
     151  *s2 = *sU / 2;
     152  *pi1 = 0.5;
     153  *pi2 = 0.5;
     154
     155  float g1,g2,norm;
     156  float w1,w2;
     157
     158  float KMM_TOLERANCE = 1e-6;
     159  int KMM_MAX_ITERATIONS = 500;
     160  float KMM_SMALL_NUMBER = 1e-5;
     161  while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
     162    *iter += 1;
     163    dL = fabs(logL_bimodal - oldL);
     164    oldL = logL_bimodal;
     165/*     if (debug == 1) { */
     166/*       fprintf(stderr,"KMM: %d %f %f %f %f (%f %f %f) (%f %f %f)\n", */
     167/*            *iter,logL_unimodal,logL_bimodal,oldL,dL, */
     168/*            *m1,*s1,*pi1, */
     169/*            *m2,*s2,*pi2); */
     170/*     } */
     171    // Expectation/P-stage
     172    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
     173      g1 = gaussian(values->data.F32[i],*m1,*s1);
     174      g2 = gaussian(values->data.F32[i],*m2,*s2);
     175      norm = (*pi1 * g1 + *pi2 * g2);
     176      P1->data.F32[i] = (*pi1 * g1) / norm;
     177      P2->data.F32[i] = (*pi2 * g2) / norm;
     178    }
     179    // Maximization/M-stage
     180    logL_bimodal = 0;
     181    w1 = 0;
     182    w2 = 0;
     183    for (i = 0; i < values->n; i++) { // Calculate log likelihood
     184      if (!((*pi1 == 0)||(*pi2 == 0))) {
     185        logL_bimodal += log(*pi1 * gaussian(values->data.F32[i],*m1,*s1) +
     186                            *pi2 * gaussian(values->data.F32[i],*m2,*s2));
     187      }
     188    }
     189    *m1 = 0;
     190    *m2 = 0;
     191    *s1 = 0;
     192    *s2 = 0;
     193    for (i = 0; i < values->n; i++) { // Calculate new means
     194      *m1 += values->data.F32[i] * P1->data.F32[i];
     195      *m2 += values->data.F32[i] * P2->data.F32[i];
     196
     197      w1 += P1->data.F32[i];
     198      w2 += P2->data.F32[i];
     199    }
     200    *m1 /= w1;
     201    *m2 /= w2;
     202    for (i = 0; i < values->n; i++) { // Calculate new sigmas
     203      *s1 += pow(values->data.F32[i] - *m1,2) * P1->data.F32[i];
     204      *s2 += pow(values->data.F32[i] - *m2,2) * P2->data.F32[i];
     205    }
     206    *s1 = sqrt(*s1 / w1);
     207    *s2 = sqrt(*s2 / w2);
     208
     209    *pi1 = w1 / values->n;
     210    *pi2 = w2 / values->n;
     211
     212    if (!isfinite(*pi1)) { // finite checks
     213      *pi1 = 0.0;
     214    }
     215    if (!isfinite(*pi2)) { // finite checks
     216      *pi2 = 0.0;
     217    }
     218    if (*s1 == 0) { // sigma may not be zero
     219      *s1 = KMM_SMALL_NUMBER * *m1;
     220    }
     221    if (*s2 == 0) { // sigma may not be zero
     222      *s2 = KMM_SMALL_NUMBER * *m2;
     223    }
     224  } // End EM phase
     225
     226  // Calculate Punimodal
     227  double lambda = -2.0 * (logL_unimodal - logL_bimodal);
     228  int    df     = 2 + 2 * 1;
     229  if (lambda > 0) {
     230    *Punimodal = gsl_cdf_chisq_Q(lambda,df);
     231  }
     232  else {
     233    *Punimodal = 1.0;
     234  }
     235  psFree(P1);
     236  psFree(P2);
     237}
     238
     239static void KMMFindPopular(const psVector *values, float *Punimodal, float *mean, float *sigma, float *pi) {
     240  float KMM_MINIMUM_PVALUE = 0.05;
     241  float mU,sU;
     242  float pi1,m1,s1,pi2,m2,s2;
     243  int iter;
     244
     245  KMMcalculate(values,Punimodal,&iter,
     246               &mU,&sU,
     247               &pi1,&m1,&s1,
     248               &pi2,&m2,&s2);
     249
     250  if (*Punimodal > KMM_MINIMUM_PVALUE) {
     251    // Is unimodal
     252    *mean = mU;
     253    *sigma = sU;
     254    *pi = 1.0;
     255  }
     256  else {
     257    if (pi1 >= pi2) {
     258      *mean = m1;
     259      *sigma = s1;
     260      *pi = pi1;
     261    }
     262    else {
     263      *mean = m2;
     264      *sigma = s2;
     265      *pi = pi2;
     266    }
     267  } 
     268}
     269
     270                               
    104271
    105272
    … …  
    287454    return;
    288455}
    289 
     456#if (0)
     457static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
     458  float KMM_MINIMUM_PVALUE = 0.05;
     459  float mU,sU;
     460  float Punimodal,pi1,m1,s1,pi2,m2,s2;
     461  int iter;
     462  int j,k;
     463
     464  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
     465  k = 0;
     466  for (j = 0; j < inputs->n; j++) {
     467    pmStackData *data = inputs->data[j]; // Stack data of interest
     468    if (!data) {
     469      k++;
     470      continue;
     471    }
     472    psImage *image = data->readout->image; // Image of interest
     473    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
     474    values->data.F32[j - k] = image->data.F32[yIn][xIn];
     475  }
     476 
     477  KMMcalculate(values,&Punimodal,&iter,
     478               &mU,&sU,
     479               &pi1,&m1,&s1,
     480               &pi2,&m2,&s2);
     481  //  fprintf(stderr,
     482/*   psTrace("psModules.imcombine",3, */
     483
     484   CHECKPIX(x, y,
     485          "KMM Unpopular Test: %d,%d: Puni: %g in %d",x,y,Punimodal,iter); 
     486  if (Punimodal < KMM_MINIMUM_PVALUE) {
     487    int i;
     488    float g1,g2,norm;
     489    float P1,P2;
     490
     491    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
     492      g1 = gaussian(values->data.F32[i],m1,s1);
     493      g2 = gaussian(values->data.F32[i],m2,s2);
     494      norm = (pi1 * g1 + pi2 * g2);
     495      P1 = (pi1 * g1) / norm;
     496      P2 = (pi2 * g2) / norm;
     497
     498      CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %f(%d): %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n",
     499               x, y,
     500               Punimodal,iter,
     501               i, values->data.F32[i],
     502               P1,m1,s1,pi1,
     503               P2,m2,s2,pi2,
     504               (pi1 > pi2)&&(P1 < P2),
     505               (pi1 < pi2)&&(P1 > P2));
     506      if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2
     507        combineMarkReject(inputs,x,y,i);
     508      }
     509      if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1
     510        combineMarkReject(inputs,x,y,i);
     511      }
     512    }
     513  }
     514  psFree(values);
     515  // else do nothing.
     516}
     517
     518static void KMMRejectBright(const psArray *inputs, int x, int y) {
     519  float KMM_MINIMUM_PVALUE = 0.05;
     520  float mU,sU;
     521  float Punimodal,pi1,m1,s1,pi2,m2,s2;
     522  int iter;
     523  int j;
     524
     525  psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
     526  for (j = 0; j < inputs->n; j++) {
     527    pmStackData *data = inputs->data[j]; // Stack data of interest
     528    psImage *image = data->readout->image; // Image of interest
     529    int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
     530    values->data.F32[j] = image->data.F32[yIn][xIn];
     531  }
     532 
     533  KMMcalculate(values,&Punimodal,&iter,
     534               &mU,&sU,
     535               &pi1,&m1,&s1,
     536               &pi2,&m2,&s2);
     537  if (Punimodal < KMM_MINIMUM_PVALUE) {
     538    int i;
     539    float g1,g2,norm;
     540    float P1,P2;
     541
     542    for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
     543      g1 = gaussian(values->data.F32[i],m1,s1);
     544      g2 = gaussian(values->data.F32[i],m2,s2);
     545      norm = (pi1 * g1 + pi2 * g2);
     546      P1 = (pi1 * g1) / norm;
     547      P2 = (pi2 * g2) / norm;
     548
     549      if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1
     550        combineMarkReject(inputs,x,y,i);
     551      }
     552      if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2
     553        combineMarkReject(inputs,x,y,i);
     554      }
     555    }
     556  }
     557  psFree(values);
     558  // else do nothing.
     559}
     560#endif
    290561
    291562// Extract vectors for simple combination/rejection operations
    … …  
    430701    if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
    431702        for (int i = 0; i < numGood; i++) {
    432             fprintf(stderr, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
     703          fprintf(stderr,"Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f %d %x %x -> %x %x\n",
    433704                    i, x, y, pixelSources->data.U16[i], pixelData->data.F32[i], pixelVariances->data.F32[i],
    434705                    addVariance->data.F32[i], pixelWeights->data.F32[i], pixelExps->data.F32[i],
    … …  
    571842    psVector *pixelLimits = buffer->limits; // Is the pixel suspect?
    572843
     844    // KMM values;
     845    float Punimodal,KMMmean,KMMsigma,KMMpi;
     846   
    573847    // Set up rejection limits
    574848    float rej2 = PS_SQR(rej); // Rejection level squared
    … …  
    577851        // Using squared rejection limit because it's cheaper than sqrts
    578852        double sumWeights = 0.0;
     853
     854        // Determine the systematic error from the most popular population in the sample
     855        KMMFindPopular(pixelData,&Punimodal,&KMMmean,&KMMsigma,&KMMpi);
     856        CHECKPIX(x,y,"KMM Popularity Contest: (%d,%d) Puni: %g Mean: %f Sigma %f Pi: %f\n",
     857                 x,y,Punimodal,KMMmean,KMMsigma,KMMpi);
     858       
    579859        for (int i = 0; i < num; i++) {
    580860            sumWeights += pixelWeights->data.F32[i];
    … …  
    582862        for (int i = 0; i < num; i++) {
    583863            // Systematic error contributes to the rejection level
    584             float sysVar = PS_SQR(sys * pixelData->data.F32[i]);
     864/*             float sysVar = PS_SQR(sys * pixelData->data.F32[i]); */
     865          float sysVar = PS_SQR(KMMsigma);
    585866            CHECKPIX(x, y, "Variance %d (%d), pixel %d,%d: %f %f %f\n",
    586867                     i, pixelSources->data.U16[i], x, y,
    … …  
    7311012          default: {
    7321013              if (useVariance) {
    733                   float median = combinationWeightedOlympic(pixelData, pixelWeights,
    734                                                             olympic, buffer->sort); // Median for stack
     1014                float median;
     1015                if (Punimodal > 0.05) {
     1016                  median = combinationWeightedOlympic(pixelData, pixelWeights,
     1017                                                      olympic, buffer->sort); // Median for stack
     1018                }
     1019                else {
     1020                  median = KMMmean;
     1021                }
    7351022                  CHECKPIX(x, y, "Flag with variance pixel %d,%d: median = %f\n", x, y, median);
    7361023                  float worst = -INFINITY; // Largest deviation
    7371024                  for (int j = 0; j < num; j++) {
     1025                   
    7381026                      float diff = pixelData->data.F32[j] - median; // Difference from expected
    7391027                      CHECKPIX(x, y, "Testing input %d for pixel %d,%d: %f\n", j, x, y, diff);
    … …  
    7421030                      // pixelVariances includes the rejection limit, from above
    7431031                      float diff2 = PS_SQR(diff); // Square difference
     1032                      CHECKPIX(x,y, "Input %d, pixel %d,%d (%" PRIu16 "): %f %f (%f) %f %f :: %f %f %f %f\n",
     1033                              i, x, y, pixelSources->data.U16[j], pixelData->data.F32[j], pixelVariances->data.F32[j],
     1034                              1.0, pixelWeights->data.F32[j], 1.0,
     1035                              pixelLimits->data.F32[j], diff2, diff2 / pixelLimits->data.F32[j],worst);
     1036
    7441037                      if (diff2 > pixelLimits->data.F32[j]) {
    7451038                          float dev = diff2 / pixelLimits->data.F32[j]; // Deviation
    … …  
    9801273}
    9811274
    982 // KMM functions to do bimodality rejection of pixels
    983 
    984 float gaussian(float x, float m, float s) {
    985   return(pow(s * sqrt(2 * M_PI),-1) * exp(-0.5 * pow( (x - m) / s, 2)));
    986 }
    987 
    988 static void KMMcalculate(const psVector *values,
    989                          float *Punimodal,int *iter,
    990                          float *pi1, float *m1, float *s1,
    991                          float *pi2, float *m2, float *s2) {
    992   double logL_bimodal = 0, logL_unimodal;
    993   float mU,sU;
    994   psVector *P1 = psVectorAlloc(values->n,PS_TYPE_F32);
    995   psVector *P2 = psVectorAlloc(values->n,PS_TYPE_F32);
    996   int i;
    997  
    998   // Calculate unimodal properties
    999   mU = 0;
    1000   sU = 0;
    1001   logL_unimodal = 0;
    1002   for (i = 0; i < values->n; i++) { // Calculate mean
    1003     mU += values->data.F32[i];
    1004   }
    1005   mU /= values->n;
    1006   for (i = 0; i < values->n; i++) { // Calculate sigma
    1007     sU += pow(values->data.F32[i] - mU,2);
    1008   }
    1009   sU = sqrt(sU / values->n);
    1010   for (i = 0; i < values->n; i++) { // Calculate log likelihood
    1011     logL_unimodal += log(gaussian(values->data.F32[i],mU,sU));
    1012   }
    1013 
    1014   // Do EM loop
    1015   float dL = 0;
    1016   float oldL = -999;
    1017   *iter = 0;
    1018   logL_bimodal = logL_unimodal;
    1019   *m1 = mU - 3 * sU;
    1020   *m2 = mU + 3 * sU;
    1021   *s1 = sU / 2;
    1022   *s2 = sU / 2;
    1023   *pi1 = 0.5;
    1024   *pi2 = 0.5;
    1025 
    1026   float g1,g2,norm;
    1027   float w1,w2;
    1028 
    1029   float KMM_TOLERANCE = 1e-6;
    1030   int KMM_MAX_ITERATIONS = 500;
    1031   float KMM_SMALL_NUMBER = 1e-5;
    1032   while (((dL > KMM_TOLERANCE)||(*iter < 3))&&(*iter < KMM_MAX_ITERATIONS)) {
    1033     *iter += 1;
    1034     dL = fabs(logL_bimodal - oldL);
    1035     oldL = logL_bimodal;
    1036 
    1037     // Expectation/P-stage
    1038     for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
    1039       g1 = gaussian(values->data.F32[i],*m1,*s1);
    1040       g2 = gaussian(values->data.F32[i],*m2,*s2);
    1041       norm = (*pi1 * g1 + *pi2 * g2);
    1042       P1->data.F32[i] = (*pi1 * g1) / norm;
    1043       P2->data.F32[i] = (*pi2 * g2) / norm;
    1044     }
    1045     // Maximization/M-stage
    1046     logL_bimodal = 0;
    1047     w1 = 0;
    1048     w2 = 0;
    1049     for (i = 0; i < values->n; i++) { // Calculate log likelihood
    1050       if (!((*pi1 == 0)||(*pi2 == 0))) {
    1051         logL_bimodal += log(*pi1 * gaussian(values->data.F32[i],*m1,*s1) +
    1052                             *pi2 * gaussian(values->data.F32[i],*m2,*s2));
    1053       }
    1054     }
    1055     *m1 = 0;
    1056     *m2 = 0;
    1057     *s1 = 0;
    1058     *s2 = 0;
    1059     for (i = 0; i < values->n; i++) { // Calculate new means
    1060       *m1 += values->data.F32[i] * P1->data.F32[i];
    1061       *m2 += values->data.F32[i] * P2->data.F32[i];
    1062 
    1063       w1 += P1->data.F32[i];
    1064       w2 += P2->data.F32[i];
    1065     }
    1066     *m1 /= w1;
    1067     *m2 /= w2;
    1068     for (i = 0; i < values->n; i++) { // Calculate new sigmas
    1069       *s1 += pow(values->data.F32[i] - *m1,2) * P1->data.F32[i];
    1070       *s2 += pow(values->data.F32[i] - *m2,2) * P2->data.F32[i];
    1071     }
    1072     *s1 = sqrt(*s1 / w1);
    1073     *s2 = sqrt(*s2 / w2);
    1074 
    1075     *pi1 = w1 / values->n;
    1076     *pi2 = w2 / values->n;
    1077 
    1078     if (!isfinite(*pi1)) { // finite checks
    1079       *pi1 = 0.0;
    1080     }
    1081     if (!isfinite(*pi2)) { // finite checks
    1082       *pi2 = 0.0;
    1083     }
    1084     if (*s1 == 0) { // sigma may not be zero
    1085       *s1 = KMM_SMALL_NUMBER * *m1;
    1086     }
    1087     if (*s2 == 0) { // sigma may not be zero
    1088       *s2 = KMM_SMALL_NUMBER * *m2;
    1089     }
    1090   } // End EM phase
    1091 
    1092   // Calculate Punimodal
    1093   double lambda = -2.0 * (logL_unimodal - logL_bimodal);
    1094   int    df     = 2 + 2 * 1;
    1095   if (lambda > 0) {
    1096     *Punimodal = gsl_cdf_chisq_Q(lambda,df);
    1097   }
    1098   else {
    1099     *Punimodal = 1.0;
    1100   } 
    1101 }
    1102 
    1103 static void KMMRejectUnpopular(const psArray *inputs, int x, int y) {
    1104   float KMM_MINIMUM_PVALUE = 0.05;
    1105   float Punimodal,pi1,m1,s1,pi2,m2,s2;
    1106   int iter;
    1107   int j;
    1108 
    1109   psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
    1110   for (j = 0; j < inputs->n; j++) {
    1111     pmStackData *data = inputs->data[j]; // Stack data of interest
    1112     psImage *image = data->readout->image; // Image of interest
    1113     int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
    1114     values->data.F32[j] = image->data.F32[yIn][xIn];
    1115   }
    1116  
    1117   KMMcalculate(values,&Punimodal,&iter,
    1118                &pi1,&m1,&s1,
    1119                &pi2,&m2,&s2);
    1120   //  fprintf(stderr,
    1121   psTrace("psModules.imcombine",3,
    1122           "KMM Unpopular Test: %d,%d: Puni: %f in %d",x,y,Punimodal,iter);
    1123   //  CHECKPIX(x, y, "
    1124  
    1125   if (Punimodal < KMM_MINIMUM_PVALUE) {
    1126     int i;
    1127     float g1,g2,norm;
    1128     float P1,P2;
    1129 
    1130     for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
    1131       g1 = gaussian(values->data.F32[i],m1,s1);
    1132       g2 = gaussian(values->data.F32[i],m2,s2);
    1133       norm = (pi1 * g1 + pi2 * g2);
    1134       P1 = (pi1 * g1) / norm;
    1135       P2 = (pi2 * g2) / norm;
    1136 
    1137       CHECKPIX(x, y, "KMM Unpopular Rejection: %d,%d: %d %f %f:(%f %f %f ) %f:(%f %f %f) rejection? %d %d\n",
    1138                x, y, i, values->data.F32[i],
    1139                P1,m1,s1,pi1,
    1140                P2,m2,s2,pi2,
    1141                (pi1 > pi2)&&(P1 < P2),
    1142                (pi1 < pi2)&&(P1 > P2));
    1143       if ((pi1 > pi2)&&(P1 < P2)) { // mode 1 is more popular, but this element belongs to mode 2
    1144         combineMarkReject(inputs,x,y,i);
    1145       }
    1146       if ((pi1 < pi2)&&(P1 > P2)) { // mode 2 is more popular, but this element belongs to mode 1
    1147         combineMarkReject(inputs,x,y,i);
    1148       }
    1149     }
    1150   }
    1151   psFree(values);
    1152   // else do nothing.
    1153 }
    1154 
    1155 static void KMMRejectBright(const psArray *inputs, int x, int y) {
    1156   float KMM_MINIMUM_PVALUE = 0.05;
    1157   float Punimodal,pi1,m1,s1,pi2,m2,s2;
    1158   int iter;
    1159   int j;
    1160 
    1161   psVector *values = psVectorAlloc(inputs->n, PS_TYPE_F32);
    1162   for (j = 0; j < inputs->n; j++) {
    1163     pmStackData *data = inputs->data[j]; // Stack data of interest
    1164     psImage *image = data->readout->image; // Image of interest
    1165     int xIn = x - data->readout->col0, yIn = y - data->readout->row0; // Coordinates on input readout
    1166     values->data.F32[j] = image->data.F32[yIn][xIn];
    1167   }
    1168  
    1169   KMMcalculate(values,&Punimodal,&iter,
    1170                &pi1,&m1,&s1,
    1171                &pi2,&m2,&s2);
    1172   if (Punimodal < KMM_MINIMUM_PVALUE) {
    1173     int i;
    1174     float g1,g2,norm;
    1175     float P1,P2;
    1176 
    1177     for (i = 0; i < values->n; i++) { // Calculate probabilities for each mode
    1178       g1 = gaussian(values->data.F32[i],m1,s1);
    1179       g2 = gaussian(values->data.F32[i],m2,s2);
    1180       norm = (pi1 * g1 + pi2 * g2);
    1181       P1 = (pi1 * g1) / norm;
    1182       P2 = (pi2 * g2) / norm;
    1183 
    1184       if ((m1 > m2)&&(P1 > P2)) { // m1 is larger, and this element belongs to mode 1
    1185         combineMarkReject(inputs,x,y,i);
    1186       }
    1187       if ((m1 < m2)&&(P1 < P2)) { // m2 is larger, and this element belongs to mode 2
    1188         combineMarkReject(inputs,x,y,i);
    1189       }
    1190     }
    1191   }
    1192   // else do nothing.
    1193 }
    1194                                
    1195                                
    1196    
    1197  
    11981275 
    11991276 
    … …  
    13531430    }
    13541431
    1355     // Pre-reject inputs using KMM bimodality test.
    1356     if (1)  {
    1357 /*       KMMRejectUnpopular(input,x,y); */
    1358       rejection = true;
    1359     }
     1432/*     // Pre-reject inputs using KMM bimodality test. */
     1433/*     if (1)  { */
     1434/* /\*       KMMRejectUnpopular(input,x,y); *\/ */
     1435/*       rejection = true; */
     1436/*     } */
    13601437   
    13611438    // Set up rejection list
    … …  
    13691446        for (int x = minInputCols; x < maxInputCols; x++) {
    13701447
    1371             CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
    1372 
    1373             // Pre-reject inputs using KMM bimodality test.
    1374           if (1)  {
    1375             KMMRejectUnpopular(input,x,y);
    1376 /*          rejection = true; */
    1377           }
    1378           else {
    1379             KMMRejectBright(input,x,y);
    1380           }
     1448/*          CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection); */
     1449
     1450/*          // Pre-reject inputs using KMM bimodality test. */
     1451/*        if (1)  { */
     1452/*          KMMRejectUnpopular(input,x,y); */
     1453/* /\*      rejection = true; *\/ */
     1454/*        } */
     1455/*        else { */
     1456/*          KMMRejectBright(input,x,y); */
     1457/*        } */
    13811458
    13821459#ifdef TESTING
Note: See TracChangeset for help on using the changeset viewer.