IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 29449


Ignore:
Timestamp:
Oct 17, 2010, 9:10:35 AM (16 years ago)
Author:
eugene
Message:

generate a weight window based on the fwhms; pass around fwhm to ensure the window can be calculated; plug some leaks; better test outputs

Location:
branches/eam_branches/ipp-20100823/psModules/src/imcombine
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmPSFEnvelope.c

    r29004 r29449  
    380380
    381381        // measure the source moments: tophat windowing, no pixel S/N cutoff
    382         // XXX probably should be passing the maskVal to this function so we can pass it along here...
    383         if (!pmSourceMoments(source, maxRadius, 0.0, 1.0, maskVal)) {
     382        if (!pmSourceMoments(source, maxRadius, 0.0, 0.0, maskVal)) {
    384383            // Can't do anything about it; limp along as best we can
    385384            psErrorClear();

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c

    r29170 r29449  
    1919
    2020//#define USE_WEIGHT                      // Include weight (1/variance) in equation?
    21 //#define USE_WINDOW                      // Include weight (1/variance) in equation?
     21
     22// XXX TEST:
     23# define USE_WINDOW                      // window to avoid neighbor contamination
    2224
    2325# define PENALTY false
    … …  
    752754    stamps->normValue2 = stats->robustMedian;
    753755
    754     fprintf (stderr, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2);
     756    psLogMsg ("psModules.imcombine", PS_LOG_INFO, "norm (1): %f (%f)\n", stamps->normValue, stamps->normValue2);
    755757
    756758    psFree(stats);
    … …  
    797799    stamp->normSquare2 = normSquare2;
    798800
    799     fprintf (stderr, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
     801    psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "normValue: %f %f %f  (%f %f)\n", normI1, normI2, stamp->norm, normSquare1, normSquare2);
    800802
    801803    return true;
    … …  
    10891091        }
    10901092
     1093#if 1
     1094        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
     1095        psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
     1096        psVectorWriteFile("sumVector.dat", sumVector);
     1097        psFree (save);
     1098#endif
     1099
    10911100        psVector *solution = NULL;                       // Solution to equation!
    10921101        solution = psVectorAlloc(numParams, PS_TYPE_F64);
    10931102        psVectorInit(solution, 0);
    10941103
     1104        // XXX TEST: try some constraint on the svd solution
     1105        // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    10951106        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
    10961107
    … …  
    11491160        }
    11501161
    1151 #if 0
     1162#if 1
    11521163        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
    11531164        psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
    … …  
    12641275    if (!isfinite(dTotal)) return false;
    12651276
    1266     fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum);
     1277    // fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum);
    12671278    psVectorAppend(fResSigma, sigma/sum);
    12681279    psVectorAppend(fResOuter, dOuter/sum);

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionKernels.c

    r29148 r29449  
    258258    }
    259259
    260 #if 1
     260#if 0
    261261    {
    262262        double Sum = 0.0;   // Sum of kernel component

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.c

    r29215 r29449  
    12021202
    12031203bool pmSubtractionParamsScale(int *kernelSize, int *stampSize, psVector *widths,
    1204                               float fwhm1, float fwhm2, float scaleRef, float scaleMin, float scaleMax)
     1204                              float scaleRef, float scaleMin, float scaleMax)
    12051205{
    12061206    PS_ASSERT_PTR_NON_NULL(kernelSize, false);
    … …  
    12081208    PS_ASSERT_VECTOR_NON_NULL(widths, false);
    12091209    PS_ASSERT_VECTOR_TYPE(widths, PS_TYPE_F32, false);
    1210     PS_ASSERT_FLOAT_LARGER_THAN(fwhm1, 0.0, false);
    1211     PS_ASSERT_FLOAT_LARGER_THAN(fwhm2, 0.0, false);
    12121210    PS_ASSERT_FLOAT_LARGER_THAN(scaleRef, 0.0, false);
    12131211    PS_ASSERT_FLOAT_LARGER_THAN(scaleMin, 0.0, false);
    … …  
    12151213    PS_ASSERT_FLOAT_LARGER_THAN(scaleMax, scaleMin, false);
    12161214
     1215    float fwhm1;
     1216    float fwhm2;
     1217
     1218    pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     1219    psAssert(isfinite(fwhm1), "fwhm 1 not set");
     1220    psAssert(isfinite(fwhm2), "fwhm 2 not set");
     1221   
    12171222    // float diff = sqrtf(PS_SQR(PS_MAX(fwhm1, fwhm2)) - PS_SQR(PS_MIN(fwhm1, fwhm2))); // Difference
    12181223    float scale = PS_MAX(fwhm1, fwhm2) / scaleRef;      // Scaling factor
    12191224
    1220     // XXX save these values in a static for later use
    1221     pmSubtractionSetFWHMs(fwhm1, fwhm2);
    1222 
    12231225    if (isfinite(scaleMin) && scale < scaleMin) {
    12241226        scale = scaleMin;

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionMatch.h

    r29004 r29449  
    104104    int *stampSize,                     ///< Half-size of the stamp (footprint)
    105105    psVector *widths,                   ///< ISIS widths
    106     float fwhm1, float fwhm2,           ///< FWHMs for inputs
    107106    float scaleRef,                     ///< Reference width for scaling
    108107    float scaleMin,                     ///< Minimum scaling ratio, or NAN

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.c

    r29274 r29449  
    5151    psFree(list->y);
    5252    psFree(list->flux);
     53    psFree(list->window);
    5354    psFree(list->window1);
    5455    psFree(list->window2);
    … …  
    9697
    9798    // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax);
    98     float xRaw = xMin;
    99     float yRaw = yMin;
     99    // float xRaw = xMin;
     100    // float yRaw = yMin;
    100101
    101102    // Ensure we're not going to go outside the bounds of the image
    … …  
    126127            if (subMask && subMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] &
    127128                (PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ)) {
    128               fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw);
     129                // fprintf (stderr, "%f,%f : masked\n", xRaw, yRaw);
    129130                return false;
    130131            }
    … …  
    142143    }
    143144
    144     if (!found) {
    145       fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
    146     }
     145    // if (!found) {
     146    //  fprintf (stderr, "%f,%f : fails flux test\n", xRaw, yRaw);
     147    // }
     148
    147149    return found;
    148150}
    … …  
    245247    list->normFrac = normFrac;
    246248    list->normValue = NAN;
     249    list->window = NULL;
    247250    list->window1 = NULL;
    248251    list->window2 = NULL;
    … …  
    269272    out->y = NULL;
    270273    out->flux = NULL;
     274    out->window = psMemIncrRefCounter(in->window);
    271275    out->window1 = psMemIncrRefCounter(in->window1);
    272276    out->window2 = psMemIncrRefCounter(in->window2);
    … …  
    487491                                            yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
    488492#endif
    489                     if (0 && goodStamp) {
    490                         fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n",
    491                                  region->x0 + size, xStamp, region->x1 - size,
    492                                  region->y0 + size, yStamp, region->y1 - size);
    493                     }
    494493                }
    495494            } else {
    … …  
    684683    psImageInit(stamps->window2->image, 0.0);
    685684
     685    // Generate a weighting window based on the fwhms (20% larger than the largest)
     686    {
     687        float fwhm1, fwhm2;
     688
     689        // XXX this is annoyingly hack-ish
     690        pmSubtractionGetFWHMs(&fwhm1, &fwhm2);
     691
     692        float sigma = 1.5 * PS_MAX(fwhm1, fwhm2) / 2.35;
     693
     694        psFree (stamps->window);
     695        stamps->window = psKernelAlloc(-size, size, -size, size);
     696        psImageInit(stamps->window->image, 0.0);
     697
     698        for (int y = -size; y <= size; y++) {
     699            for (int x = -size; x <= size; x++) {
     700                stamps->window->kernel[y][x] = exp(-0.5*(x*x + y*y)/(sigma*sigma));
     701            }
     702        }
     703    }
     704
    686705    // generate normalizations for each stamp
    687706    psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
    … …  
    759778    {
    760779        psFits *fits = NULL;
    761         fits = psFitsOpen ("window1.fits", "w");
     780        fits = psFitsOpen ("window1.raw.fits", "w");
    762781        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
    763782        psFitsClose (fits);
    764         fits = psFitsOpen ("window2.fits", "w");
     783        fits = psFitsOpen ("window2.raw.fits", "w");
    765784        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
    766785        psFitsClose (fits);
    … …  
    818837            // interpolate to the radius at which delta2 is normFrac:
    819838            stamps->normWindow1 = radius - (stamps->normFrac - delta1) / (delta1o - delta1);
    820             fprintf (stderr, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
     839            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 1 (%f : %f)\n", stamps->normWindow1, radius, within1, delta1);
    821840            done1 = true;
    822841        }
    … …  
    825844            // interpolate to the radius at which delta2 is normFrac:
    826845            stamps->normWindow2 = radius - (stamps->normFrac - delta2) / (delta2o - delta2);
    827             fprintf (stderr, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
     846            psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "choosing %f (%d) for 2 (%f : %f)\n", stamps->normWindow2, radius, within2, delta2);
    828847            done2 = true;
    829848        }
    … …  
    873892    {
    874893        psFits *fits = NULL;
    875         fits = psFitsOpen ("window1.fits", "w");
     894        fits = psFitsOpen ("window1.norm.fits", "w");
    876895        psFitsWriteImage (fits, NULL, stamps->window1->image, 0, NULL);
    877896        psFitsClose (fits);
    878         fits = psFitsOpen ("window2.fits", "w");
     897        fits = psFitsOpen ("window2.norm.fits", "w");
    879898        psFitsWriteImage (fits, NULL, stamps->window2->image, 0, NULL);
    880899        psFitsClose (fits);
    … …  
    962981    psAssert (isfinite(kernels->penalties2->data.F32[index]), "invalid penalty");
    963982
    964     fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
     983    // fprintf(stderr, "penalty1: %f, penalty2: %f\n", penalty1, penalty2);
    965984
    966985    return true;
    … …  
    9891008    penalty *= 1.0 / sum2;
    9901009
    991     if (1) {
    992         fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
     1010    if (0) {
     1011        // fprintf(stderr, "min: %lf, max: %lf, moment: %lf, flux^2: %lf\n", min, max, penalty, sum2);
    9931012        // psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index, fwhm, uOrder, vOrder, penalty);
    9941013    }

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionStamps.h

    r29170 r29449  
    2929    psKernel *window1;                  ///< window function generated from ensemble of stamps (input 1)
    3030    psKernel *window2;                  ///< window function generated from ensemble of stamps (input 2)
     31    psKernel *window;                   ///< weighting window function (sigma = 1.1 * MAX(fwhm))
    3132    float normWindow1;                  ///< Size of window for measuring normalisation
    3233    float normWindow2;                  ///< Size of window for measuring normalisation

  • branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.c

    r29148 r29449  
    146146    }
    147147    pmVisualScaleImage(kapa1, canvas, "Subtraction_Stamps", 0, true);
     148    psFree(canvas);
    148149
    149150    pmVisualAskUser(&plotStamps);
    … …  
    281282            }
    282283        }
    283         fprintf (stderr, "kernel %d, sum %f\n", i, sum);
     284        // fprintf (stderr, "kernel %d, sum %f\n", i, sum);
    284285    }                                                   
    285286       
    286287    pmVisualScaleImage(kapa1, output, "Image", 0, true);
     288    psFree(output);
    287289    pmVisualAskUser(&plotImage);
    288290    return true;
    … …  
    309311        if (!isfinite(stamp->flux)) continue;
    310312        if (!stamp->convolutions1 && !stamp->convolutions2) continue;
    311         fprintf (stderr, "flux: %f, maxFlux: %f  ", stamp->flux, maxFlux);
     313        // fprintf (stderr, "flux: %f, maxFlux: %f  ", stamp->flux, maxFlux);
    312314        if (!maxStamp) {
    313315            maxFlux = stamp->flux;
    314316            maxStamp = stamp;
    315             fprintf (stderr, "maxStamp %d\n", i);
     317            // fprintf (stderr, "maxStamp %d\n", i);
    316318            continue;
    317319        } else {
    318             fprintf (stderr, "\n");
     320            // fprintf (stderr, "\n");
    319321        }
    320322        if (stamp->flux > maxFlux) {
    … …  
    361363                }
    362364            }
    363             fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
     365            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
    364366        }               
    365367        pmVisualScaleImage(kapa2, output, "Image", 0, true);
    … …  
    403405                }
    404406            }
    405             fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
     407            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
    406408        }               
    407409        pmVisualScaleImage(kapa2, output, "Image", 1, true);
    … …  
    412414            psFitsClose(fits);
    413415        }
     416        psFree(output);
    414417    }                                   
    415418       
    … …  
    708711    psFree (x);
    709712    psFree (y);
     713    psFree (polyValues);
    710714
    711715    pmVisualAskUser(NULL);
Note: See TracChangeset for help on using the changeset viewer.