IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 35553


Ignore:
Timestamp:
May 9, 2013, 11:33:40 AM (13 years ago)
Author:
eugene
Message:

merge changes from trunk

Location:
branches/eam_branches/ipp-20130419/ppStack/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130419/ppStack/src/ppStackConvolve.c

    r34372 r35553  
    403403                    continue;
    404404                }
    405                 if ((options->matchChi2->data.F32[i] > thresh) || ! isfinite(options->matchChi2->data.F32[i])) {
     405                if ((options->matchChi2->data.F32[i] > thresh) ||
     406                    ! isfinite(options->matchChi2->data.F32[i])) {
    406407                    numRej++;
    407408                    options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2;

  • branches/eam_branches/ipp-20130419/ppStack/src/ppStackMatch.c

    r35383 r35553  
    467467                sum += kernels->rms;
    468468                num++;
     469                //              psLogMsg("ppStack", PS_LOG_INFO, "KERNMATCHDUMP %d %d %g %g %g\n",
     470                //                       index,num,kernels->mean,kernels->rms,sum);
    469471            }
    470472            psFree(iter);
    471473            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     474
    472475        }
    473476

  • branches/eam_branches/ipp-20130419/ppStack/src/ppStackPrepare.c

    r35400 r35553  
    257257
    258258    bool mdok = false;
     259    bool  simpleClip = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.SIMPLE");
    259260    float maxFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.MAX"); // max allowed input fwhm
    260261    float clipFWHMnSig = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.CLIP.NSIGMA"); // sigma clipping of inputs
    261262    float threshFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.THRESH"); // FWHM size that we refuse to clip below
    262263    float asymmetryFWHM = psMetadataLookupF32(&mdok, recipe, "PSF.INPUT.ASYMMETRY"); // max bimodal asymmetry
     264   
    263265    psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message
    264266    bool havePSFs = false;                                // Do we have any PSFs?
    … …  
    309311
    310312#ifdef NEW_REJECTION
    311     // Do GMM test with two modes to decide where to put the break.
    312     double Punimodal;
    313     int m = 2;
    314     psVector *modes = psVectorAlloc(num,PS_TYPE_F32);
    315     psVector *means = psVectorAlloc(m,PS_TYPE_F32);
    316     psVector *S     = psVectorAlloc(m,PS_TYPE_F32);
    317     psVector *pi    = psVectorAlloc(m,PS_TYPE_F32);
    318     psImage *P      = psImageAlloc(m,num,PS_TYPE_F32);
    319     float limit;
    320     if (!psMM1DClass(options->inputSeeing,
    321                      options->inputSeeing->n,
    322                      modes,means,
    323                      S,pi,P,
    324                      2,
    325                      &Punimodal)) {
    326       // Handle error here
    327       psFree(modes);
    328       psFree(means);
    329       psFree(S);
    330       psFree(pi);
    331       psFree(P);
    332       return(false);
    333     }
    334     //    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
    335     //    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
    336 
    337     //    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
    338    
    339     //    unwrittenGMMfunction(options->inputSeeing, options->inputMask, 0xff,
    340     //           &Punimodal,
    341     //           &m1,&s1,&pi1,
    342     //           &m2,&s2,&pi2);
    343 
    344     // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
    345     if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode
    346         (num <= 4)) {       // Or we have a small number of inputs, making this statistic poor.
    347       limit = maxFWHM;
    348     }
    349     else {                  // This is a bimodal distribution
    350       if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)||  // However, both modes are equally populated
    351           (pi->data.F32[1] > pi->data.F32[0])) {                       // Or the larger FWHM mode is more populated
    352         limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1];
     313    // Begin new rejection code.
     314    float limit = INFINITY; // Default to a big value.
     315    if (!isfinite(threshFWHM)) { // If this isn't set, initialize it to zero.  No inputs may be clipped if smaller than this
     316      threshFWHM = 0.0;
     317    }
     318    if (!isfinite(maxFWHM)) { // No maxFWHM was set, accept all inputs
     319      maxFWHM = INFINITY;
     320    }
     321    limit = maxFWHM; // Initialize the limit to be the maxFWHM
     322    if (limit < threshFWHM) {
     323      limit = threshFWHM;
     324    }
     325    if (!isfinite(clipFWHMnSig)) { // We cannot do a sigma clip with nSig, so use the maxFWHM to set the limit.
     326      psLogMsg("ppStack",PS_LOG_INFO,
     327               "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)",
     328               limit,
     329               clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM);
     330    }
     331    if (simpleClip) { // Do a sigma clip like the old rejection code.
     332      psStats *fwhmStats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     333      psVectorStats (fwhmStats, options->inputSeeing, NULL, options->inputMask, 0xff);
     334      psLogMsg("ppStack", PS_LOG_INFO, "Input FWHMs : %f +/- %f", fwhmStats->clippedMean, fwhmStats->clippedStdev);
     335     
     336      options->clippedMean =  fwhmStats->clippedMean;
     337      options->clippedStdev = fwhmStats->clippedStdev;
     338      limit = options->clippedMean + clipFWHMnSig * options->clippedStdev;
     339      if (limit < threshFWHM) {
     340        limit = threshFWHM;
    353341      }
    354       else {                                   // The smaller FWHM mode is more populated
    355         limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0];
     342      if (limit > maxFWHM) {
     343        limit = maxFWHM;
    356344      }
    357     }
    358     if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
    359     if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
    360 
     345      psLogMsg("ppStack",PS_LOG_INFO,
     346               "PSF FWHM distribution: Used simple clip method.  Limit set to %f.",
     347               limit);
     348      psFree(fwhmStats);
     349    }
     350    else { // Do the mixture model rejection.
     351      if(!isfinite(asymmetryFWHM)) { // Or not, because the parameters aren't set.
     352        psLogMsg("ppStack",PS_LOG_INFO,
     353                 "PSF FWHM distribution: Unable to determine limit based on recipe. Setting to %f. (clipFWHMnSig: %f maxFWHM: %f threshFWHM: %f asymmetryFWHM: %f)",
     354                 limit,
     355                 clipFWHMnSig,maxFWHM,threshFWHM,asymmetryFWHM);
     356      }
     357      else {
     358        // Do GMM test with two modes to decide where to put the break.
     359        double Punimodal;
     360        int m = 2;
     361        psVector *modes = psVectorAlloc(num,PS_TYPE_F32);
     362        psVector *means = psVectorAlloc(m,PS_TYPE_F32);
     363        psVector *S     = psVectorAlloc(m,PS_TYPE_F32);
     364        psVector *pi    = psVectorAlloc(m,PS_TYPE_F32);
     365        psImage *P      = psImageAlloc(m,num,PS_TYPE_F32);
     366       
     367        if (!psMM1DClass(options->inputSeeing,
     368                         options->inputSeeing->n,
     369                         modes,means,
     370                         S,pi,P,
     371                         2,
     372                         &Punimodal)) {
     373          // Handle error here
     374          psFree(modes);
     375          psFree(means);
     376          psFree(S);
     377          psFree(pi);
     378          psFree(P);
     379          return(false);
     380        }
     381        //    fprintf(stderr,"means: %g %g\n",means->data.F32[0],means->data.F32[1]);
     382        //    fprintf(stderr,"sigma: %g %g \n",S->data.F32[0],S->data.F32[1]);
     383        //    fprintf(stderr,"pi:    %g %g\n",pi->data.F32[0],pi->data.F32[1]);
     384       
     385        // Use Gaussian mixture model analysis of the FWHM distribution to decide an optional FWHM limit
     386        if ((Punimodal > 0.5)|| // This distribution is best represented by a single mode
     387            (num <= 4)) {       // Or we have a small number of inputs, making this statistic poor.
     388          limit = maxFWHM;
     389        }
     390        else {                  // This is a bimodal distribution
     391          if ((fabs(pi->data.F32[0] - pi->data.F32[1]) < asymmetryFWHM)||  // However, both modes are equally populated
     392              (pi->data.F32[1] > pi->data.F32[0])) {                       // Or the larger FWHM mode is more populated
     393            limit = means->data.F32[1] + clipFWHMnSig * S->data.F32[1];
     394          }
     395          else {                                   // The smaller FWHM mode is more populated
     396            limit = means->data.F32[0] + clipFWHMnSig * S->data.F32[0];
     397          }
     398        }
     399        if (limit > maxFWHM)    { limit = maxFWHM; }    // We should not be larger than our max
     400        if (limit < threshFWHM) { limit = threshFWHM; } // Nor smaller than our min
     401        psLogMsg("ppStack",PS_LOG_INFO,
     402                 "PSF FWHM distribution: limit: %f (%f %f %f) (%f %f %f) %f",
     403                 limit,means->data.F32[0],S->data.F32[0],pi->data.F32[0],
     404                 means->data.F32[1],S->data.F32[1],pi->data.F32[1],
     405                 Punimodal);
     406        psFree(means);
     407        psFree(modes);
     408        psFree(S);
     409        psFree(pi);
     410        psFree(P);
     411      } // End mixture model case
     412    }
     413    // Perform rejection using the limit set
    361414    for (int i = 0; i < num; i++) {
    362415      if (options->inputSeeing->data.F32[i] > limit) {
    … …  
    367420      }
    368421    }
    369     psFree(means);
    370     psFree(modes);
    371     psFree(S);
    372     psFree(pi);
    373     psFree(P);
    374422    // End new rejection code
    375423#endif
    … …  
    404452    if (options->convolve) {
    405453        options->psf = ppStackPSF(config, numCols, numRows, psfs, options);
    406         psFree(psfs);
     454        psFree(psfs);
    407455        if (!options->psf) {
    408456#if 1
Note: See TracChangeset for help on using the changeset viewer.