IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 2, 2007, 4:28:24 PM (19 years ago)
Author:
Paul Price
Message:

Adding function pmSubtractionOrder (choice of name is probably not the best) to determine which of the two images under consideration should be convolved to match the other. Originally was doing this by solving for a RING kernel of width 1, going both ways, and comparing the deviations for each. However, when doing stacks this didn't work (convolving the wider fake Gaussian image gave smaller deviations than convolving the narrower input image), so have settled on measuring the second moments for each stamp, and using the ratio of moments between the two images to determine which is wider; this seems to work for both subtracting and stacking images. In the process, plugged a few memory leaks, and added code to support convolving either way (or both ways; this will be useful when it comes time to code the dual convolution algorithm).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/imcombine/pmSubtractionMatch.c

    r15329 r15443  
    5050
    5151static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
    52                       const pmReadout *reference, // Reference readout
    53                       const pmReadout *input, // Input readout, or NULL to generate fake stamps
     52                      const pmReadout *ro1, // Readout 1
     53                      const pmReadout *ro2, // Readout 2
    5454                      const psImage *subMask, // Mask for subtraction, or NULL
    5555                      psImage *weight,  // Weight map
    … …  
    5757                      float threshold,  // Threshold for stamp finding
    5858                      float stampSpacing, // Spacing between stamps
    59                       float targetWidth,// Target width for fake stamps
    6059                      int size,         // Kernel half-size
    61                       int footprint     // Convolution footprint for stamps
     60                      int footprint,     // Convolution footprint for stamps
     61                      pmSubtractionMode mode // Mode for subtraction
    6262    )
    6363{
    6464    psTrace("psModules.imcombine", 3, "Finding stamps...\n");
    65     *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing);
     65    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode);
    6666    if (!*stamps) {
    6767        psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
    … …  
    7171    memCheck("  find stamps");
    7272
    73     if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) {
    74         psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
    75         return false;
    76     }
    77 
    78     memCheck("   generate stamps");
    79 
    8073    psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
    81     if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL,
    82                                     weight, footprint, size)) {
     74    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) {
    8375        psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
    8476        return false;
    … …  
    9486
    9587
    96 bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input,
     88bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2,
    9789                        int footprint, float regionSize, float stampSpacing, float threshold,
    98                         const psArray *sources, const char *stampsName, float targetWidth,
     90                        const psArray *sources, const char *stampsName,
    9991                        pmSubtractionKernelsType type, int size, int spatialOrder,
    10092                        const psVector *isisWidths, const psVector *isisOrders,
    10193                        int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs,
    10294                        int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
    103                         psMaskType maskBlank, float badFrac)
     95                        psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
    10496{
    10597    PS_ASSERT_PTR_NON_NULL(convolved, false);
    106     PS_ASSERT_PTR_NON_NULL(reference, false);
    107     PS_ASSERT_IMAGE_NON_NULL(reference->image, false);
    108     PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false);
    109     if (reference->mask) {
    110         PS_ASSERT_IMAGE_NON_NULL(reference->mask, false);
    111         PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false);
    112         PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false);
    113     }
    114     if (reference->weight) {
    115         PS_ASSERT_IMAGE_NON_NULL(reference->weight, false);
    116         PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false);
    117         PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false);
    118     }
    119     if (input) {
    120         PS_ASSERT_IMAGE_NON_NULL(input->image, false);
    121         PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false);
    122         PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false);
    123         if (input->mask) {
    124             PS_ASSERT_IMAGE_NON_NULL(input->mask, false);
    125             PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false);
    126             PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false);
    127         }
    128         if (input->weight) {
    129             PS_ASSERT_IMAGE_NON_NULL(input->weight, false);
    130             PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false);
    131             PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false);
     98    PS_ASSERT_PTR_NON_NULL(ro1, false);
     99    PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
     100    PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false);
     101    if (ro1->mask) {
     102        PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false);
     103        PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false);
     104        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false);
     105    }
     106    if (ro1->weight) {
     107        PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false);
     108        PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false);
     109        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false);
     110    }
     111    if (ro2) {
     112        PS_ASSERT_IMAGE_NON_NULL(ro2->image, false);
     113        PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false);
     114        PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false);
     115        if (ro2->mask) {
     116            PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false);
     117            PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false);
     118            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false);
     119        }
     120        if (ro2->weight) {
     121            PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false);
     122            PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false);
     123            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false);
    132124        }
    133125    } else if (!stampsName && !sources) {
    … …  
    144136    }
    145137    // stampsName may be anything
    146     // targetWidth can be just about anything (except maybe negative, but it can be NAN)
    147138    // We'll check kernel type when we allocate the kernels
    148139    PS_ASSERT_INT_POSITIVE(size, false);
    … …  
    196187    }
    197188
    198     psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight
    199     if (input) {
    200         inImage = input->image;
    201         inMask = input->mask;
    202     }
    203 
    204189    // Where does our weight map come from?
    205190    psImage *weight = NULL;             // Weight image to use
    206     if (input && input->weight) {
    207         weight = input->weight;
    208     } else if (reference->weight) {
    209         weight = reference->weight;
    210     } else if (input) {
    211         weight = input->image;
     191    if (ro1->weight && ro2 && ro2->weight) {
     192        weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight);
     193    } else if (ro1->weight) {
     194        weight = psMemIncrRefCounter(ro1->weight);
     195    } else if (ro2) {
     196        if (ro2->weight) {
     197            weight = psMemIncrRefCounter(ro2->weight);
     198        } else {
     199            weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);
     200        }
    212201    } else {
    213         weight = reference->image;
     202        weight = psMemIncrRefCounter(ro1->image);
    214203    }
    215204
    … …  
    221210    psVector *solution = NULL;          // Solution to match PSF
    222211    pmSubtractionKernels *kernels = NULL; // Kernel basis functions
    223 
    224     int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
     212    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
    225213
    226214    memCheck("start");
    227215
    228     subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT);
     216    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint,
     217                                badFrac, useFFT);
    229218    if (!subMask) {
    230219        psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
     220        psFree(weight);
    231221        return false;
    232222    }
    … …  
    260250
    261251            if (sources) {
    262                 stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing,
    263                                                            input ? 0 : 2 * footprint);
     252                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode);
    264253            } else if (stampsName && strlen(stampsName) > 0) {
    265                 // Read stamps from file
    266                 psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
    267                 const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file
    268                 psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat);
    269                 psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
    270                 if (!stampsData) {
    271                     psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
    272                     goto ERROR;
    273                 }
    274                 xStamp = stampsData->data[0];
    275                 yStamp = stampsData->data[1];
    276                 if (!input) {
    277                     fluxStamp = stampsData->data[2];
    278                 }
    279 
    280                 // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
    281                 psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    282                 psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
    283 
    284                 stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask,
    285                                                 region, stampSpacing, input ? 0 : footprint + size);
    286                 psFree(stampsData);
     254                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode);
     255            }
     256
     257            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
     258            // doesn't matter.
     259            if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing,
     260                           size, footprint, mode)) {
     261                goto MATCH_ERROR;
     262            }
     263
     264            if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) {
     265                pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode
     266                switch (newMode) {
     267                  case PM_SUBTRACTION_MODE_1:
     268                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
     269                    break;
     270                  case PM_SUBTRACTION_MODE_2:
     271                    if (mode == PM_SUBTRACTION_MODE_TARGET) {
     272                        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     273                                "Input PSF is larger than target PSF --- can't match image.");
     274                        goto MATCH_ERROR;
     275                    }
     276                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
     277                    break;
     278                  default:
     279                    psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order.");
     280                    goto MATCH_ERROR;
     281                }
     282                mode = newMode;
    287283            }
    288284
    289285            // Define kernel basis functions
    290286            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
    291                 if (!getStamps(&stamps, reference, input, subMask, weight, NULL,
    292                                threshold, stampSpacing, targetWidth, size, footprint)) {
    293                     goto ERROR;
    294                 }
    295287                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
    296                                                           stamps, footprint, optThreshold);
     288                                                          stamps, footprint, optThreshold, mode);
    297289                if (!kernels) {
    298290                    psErrorClear();
    … …  
    314306                psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
    315307
    316                 if (!getStamps(&stamps, reference, input, subMask, weight, region,
    317                                threshold, stampSpacing, targetWidth, size, footprint)) {
    318                     goto ERROR;
     308                if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,
     309                               size, footprint, mode)) {
     310                    goto MATCH_ERROR;
    319311                }
    320312
    321313                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
    322                 if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) {
     314                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
    323315                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    324                     goto ERROR;
     316                    goto MATCH_ERROR;
    325317                }
    326318
    … …  
    331323                if (!solution) {
    332324                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    333                     goto ERROR;
     325                    goto MATCH_ERROR;
    334326                }
    335327
    336328                memCheck("  solve equation");
    337329
     330                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
     331                                                                        mode); // Deviations for each stamp
     332                if (!deviations) {
     333                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     334                    goto MATCH_ERROR;
     335                }
     336
     337                memCheck("   calculate deviations");
     338
    338339                psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
    339                 numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels);
     340                numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint);
    340341                if (numRejected < 0) {
    341342                    psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
    342                     goto ERROR;
    343                 }
     343                    psFree(deviations);
     344                    goto MATCH_ERROR;
     345                }
     346                psFree(deviations);
     347
    344348                memCheck("  reject stamps");
    345349            }
    … …  
    350354                if (!solution) {
    351355                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
    352                     goto ERROR;
    353                 }
    354                 (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels);
     356                    goto MATCH_ERROR;
     357                }
     358                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
     359                                                                        mode); // Deviations for each stamp
     360                if (!deviations) {
     361                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     362                    goto MATCH_ERROR;
     363                }
     364                (void)pmSubtractionRejectStamps(stamps, deviations, subMask, footprint, NAN);
     365                psFree(deviations);
    355366            }
    356367            psFree(stamps);
    … …  
    372383                            psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
    373384                            psFree(convKernels);
    374                             goto ERROR;
     385                            goto MATCH_ERROR;
    375386                        }
    376387
    … …  
    380391                            psFree(kernel);
    381392                            psFree(convKernels);
    382                             goto ERROR;
     393                            goto MATCH_ERROR;
    383394                        }
    384395                        psFree(kernel);
    … …  
    416427
    417428            psTrace("psModules.imcombine", 2, "Convolving...\n");
    418             if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask,
    419                                        reference->image, reference->weight, subMask, maskBlank, region,
    420                                        solution, kernels, useFFT)) {
    421                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
    422                 goto ERROR;
     429            if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region,
     430                                       solution, kernels, mode, useFFT)) {
     431                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
     432                goto MATCH_ERROR;
    423433            }
    424434            psFree(kernels);
    … …  
    461471    psFree(subMask);
    462472    subMask = NULL;
     473    psFree(weight);
     474    weight = NULL;
    463475
    464476    if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
    465477        psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
    466         goto ERROR;
     478        goto MATCH_ERROR;
    467479    }
    468480
    … …  
    472484    return true;
    473485
    474 ERROR:
     486MATCH_ERROR:
    475487    psFree(region);
    476488    psFree(regionString);
    … …  
    479491    psFree(stamps);
    480492    psFree(solution);
     493    psFree(weight);
    481494    return false;
    482495}
     496
     497// Calculate the second order moments for an image
     498static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments
     499                                    int radius   // Maximum radius
     500                                    )
     501{
     502    assert(kernel && kernel->kernel);
     503
     504    int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x
     505    int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y
     506
     507    float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment)
     508    float sum = 0.0;       // Sum (zero-th moment)
     509    for (int y = yMin; y <= yMax; y++) {
     510        for (int x = xMin; x <= xMax; x++) {
     511            xCentroid += kernel->kernel[y][x] * x;
     512            yCentroid += kernel->kernel[y][x] * y;
     513            sum += kernel->kernel[y][x];
     514        }
     515    }
     516    xCentroid /= sum;
     517    yCentroid /= sum;
     518
     519    float eta20 = 0.0, eta02 = 0.0;     // Second moments
     520    for (int y = yMin; y <= yMax; y++) {
     521        float yDiff = y - yCentroid;
     522        for (int x = xMin; x <= xMax; x++) {
     523            float xDiff = x - xCentroid;
     524            eta20 += PS_SQR(xDiff) * kernel->kernel[y][x];
     525            eta02 += PS_SQR(yDiff) * kernel->kernel[y][x];
     526        }
     527    }
     528
     529    // Normalise to calculate the scale-invariant
     530    float sum2 = PS_SQR(sum);
     531    eta20 /= sum2;
     532    eta02 /= sum2;
     533    // eta11 /= sum2;
     534
     535    return eta20 + eta02;
     536}
     537
     538#if 0
     539// Calculate the deviations for a particular subtraction order
     540static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel
     541                                           pmSubtractionStampList *stamps, // Stamps to convolve
     542                                           const pmSubtractionKernels *kernels, // Kernel basis functions
     543                                           int footprint, // Stamp footprint
     544                                           pmSubtractionMode mode // Mode of subtraction
     545                                           )
     546{
     547    assert(stamps);
     548    assert(footprint >= 0);
     549    assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2);
     550
     551    if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
     552        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     553        return NULL;
     554    }
     555
     556    psVector *solution = pmSubtractionSolveEquation(NULL, stamps);
     557    if (!solution) {
     558        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
     559        return NULL;
     560    }
     561
     562    if (sumKernel) {
     563        float sum = 0.0;                // Sum of the kernel
     564        psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel
     565        for (int y = 0; y < image->numRows; y++) {
     566            for (int x = 0; x < image->numCols; x++) {
     567                sum += image->data.F32[y][x];
     568            }
     569        }
     570        psFree(image);
     571        *sumKernel = sum;
     572    }
     573
     574    psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode);
     575    psFree(solution);
     576    if (!deviations) {
     577        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
     578        return NULL;
     579    }
     580
     581    return deviations;
     582}
     583#endif
     584
     585pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius)
     586{
     587    PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR);
     588
     589    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps
     590    psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments
     591    for (int i = 0; i < stamps->num; i++) {
     592        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
     593        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
     594            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
     595            continue;
     596        }
     597        mask->data.PS_TYPE_MASK_DATA[i] = 0;
     598        moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) /
     599            subtractionOrderMoment(stamp->image2, radius);
     600    }
     601
     602    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     603    if (!psVectorStats(stats, moments, NULL, mask, 0xff)) {
     604        psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
     605        psFree(mask);
     606        psFree(moments);
     607        psFree(stats);
     608        return PM_SUBTRACTION_MODE_ERR;
     609    }
     610    psFree(moments);
     611    psFree(mask);
     612
     613    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian);
     614    pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2);
     615    psFree(stats);
     616
     617    return mode;
     618}
     619
     620
     621
Note: See TracChangeset for help on using the changeset viewer.