IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 10, 2008, 10:54:53 AM (18 years ago)
Author:
Paul Price
Message:

Turn off convolution if PSFs not provided.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStackMatch.c

    r17268 r17426  
    1313
    1414//#define TESTING
    15 //#define NO_CONVOLUTION
    1615
    1716bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels,
    … …  
    2221    assert(kernels && !*kernels);
    2322    assert(sourcesRO);
    24     assert(psf);
    2523    assert(config);
    2624
    … …  
    3735    int renormWidth = psMetadataLookupS32(&mdok, config->arguments, "RENORM.WIDTH"); // Width for renormalisation box
    3836
    39 #ifndef NO_CONVOLUTION
    40     int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
    41     float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
    42     float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
    43     int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
    44     float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
    45     int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
    46     float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
    47     pmSubtractionKernelsType type =
    48         pmSubtractionKernelsTypeFromString(psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE")); // Kernel type
    49     psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
    50     psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
    51     int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
    52     int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
    53     int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
    54     psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BLANK"),
    55                                         config); // Mask for blank reg.
    56     float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
    57     bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
    58     float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
    59     float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
    60     float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
    61     float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
    62     int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
    63 
    64     // These values are specified specifically for stacking
    65     const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps
    66 
    67     psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
    68     if (optimum) {
    69         optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    70     }
    71 
    72     psArray *sources = NULL;            // Sources in image
    73     if (sourcesRO) {
    74         sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
    75     }
    76 
    77     // Generate a fake image to match to
    78     float maxMag = -INFINITY;           // Maximum magnitude of sources
    79     for (int i = 0; i < sources->n; i++) {
    80         pmSource *source = sources->data[i]; // Source of interest
    81         if (source->psfMag > maxMag && source->psfMag != MAG_IGNORE) {
    82             maxMag = source->psfMag;
    83         }
    84     }
    85 
    86     pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
    87 
    88     if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, NULL,
    89                                   psf, powf(10.0, -0.4 * maxMag), 0, false)) {
    90         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
     37    if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
     38        assert(psf);
     39        int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order
     40        float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs
     41        float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing
     42        int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size
     43        float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps
     44        int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations
     45        float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold
     46        pmSubtractionKernelsType type =
     47            pmSubtractionKernelsTypeFromString(psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE")); // Kernel type
     48        psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths
     49        psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders
     50        int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius
     51        int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order
     52        int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel
     53        psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BLANK"),
     54                                            config); // Mask for blank reg.
     55        float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction
     56        bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters?
     57        float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search
     58        float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search
     59        float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search
     60        float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search
     61        int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search
     62
     63        // These values are specified specifically for stacking
     64        const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps
     65
     66        psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
     67        if (optimum) {
     68            optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
     69        }
     70
     71        psArray *sources = NULL;            // Sources in image
     72        if (sourcesRO) {
     73            sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES");
     74        }
     75
     76        // Generate a fake image to match to
     77        float maxMag = -INFINITY;           // Maximum magnitude of sources
     78        for (int i = 0; i < sources->n; i++) {
     79            pmSource *source = sources->data[i]; // Source of interest
     80            if (source->psfMag > maxMag && source->psfMag != MAG_IGNORE) {
     81                maxMag = source->psfMag;
     82            }
     83        }
     84
     85        pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
     86
     87        if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL,
     88                                      NULL, psf, powf(10.0, -0.4 * maxMag), 0, false)) {
     89            psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
     90            psFree(fake);
     91            psFree(optWidths);
     92            psFree(output);
     93            return false;
     94        }
     95
     96#ifdef TESTING
     97        {
     98            psFits *fits = psFitsOpen("fake.fits", "w");
     99            psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
     100            psFitsClose(fits);
     101        }
     102#endif
     103
     104        // Do the image matching
     105        if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
     106                                sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
     107                                binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
     108                                maskBlank, badFrac, PM_SUBTRACTION_MODE_1)) {
     109            psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
     110            psFree(fake);
     111            psFree(optWidths);
     112            psFree(output);
     113            return false;
     114        }
    91115        psFree(fake);
    92116        psFree(optWidths);
    93         psFree(output);
    94         return false;
    95     }
    96 
    97 #ifdef TESTING
    98     {
    99         psFits *fits = psFitsOpen("fake.fits", "w");
    100         psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
    101         psFitsClose(fits);
    102     }
    103 #endif
    104 
    105     // Do the image matching
    106     if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold,
    107                             sources, stampsName, type, size, order, widths, orders, inner, ringsOrder,
    108                             binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad,
    109                             maskBlank, badFrac, PM_SUBTRACTION_MODE_1)) {
    110         psError(PS_ERR_UNKNOWN, false, "Unable to match images.");
    111         psFree(fake);
    112         psFree(optWidths);
    113         psFree(output);
    114         return false;
    115     }
    116     psFree(fake);
    117     psFree(optWidths);
    118 
    119     // Replace original images with convolved
    120     psFree(readout->image);
    121     psFree(readout->mask);
    122     psFree(readout->weight);
    123     readout->image  = psMemIncrRefCounter(output->image);
    124     readout->mask   = psMemIncrRefCounter(output->mask);
    125     readout->weight = psMemIncrRefCounter(output->weight);
    126 
    127 #else  // NO_CONVOLUTION
    128     // Fake the convolution
    129     {
     117
     118        // Replace original images with convolved
     119        psFree(readout->image);
     120        psFree(readout->mask);
     121        psFree(readout->weight);
     122        readout->image  = psMemIncrRefCounter(output->image);
     123        readout->mask   = psMemIncrRefCounter(output->mask);
     124        readout->weight = psMemIncrRefCounter(output->weight);
     125    } else {
     126        // Fake the convolution
    130127        psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
    131128        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    132129                         PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
    133130        psFree(region);
    134         pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(size, 0, PM_SUBTRACTION_MODE_1);
     131        pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(0, 0, PM_SUBTRACTION_MODE_1);
    135132        // Set solution to delta function
    136133        kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
    … …  
    142139        psFree(kernels);
    143140    }
    144 #endif
    145141
    146142    // Extract the regions and solutions used in the image matching
Note: See TracChangeset for help on using the changeset viewer.