IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27850


Ignore:
Timestamp:
May 4, 2010, 8:59:10 AM (16 years ago)
Author:
eugene
Message:

working on the psf-matching for stack photometry analysis

Location:
trunk/psphot/src
Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotStackMatchPSFs.c

    r27848 r27850  
    55    bool status = true;
    66
    7     // select the appropriate recipe information
    8     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    9     psAssert (recipe, "missing recipe?");
    10 
    117    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    128    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    … …  
    1410    // loop over the available readouts
    1511    for (int i = 0; i < num; i++) {
    16         if (!psphotStackMatchPSFsReadout (config, view, i, recipe)) {
     12        if (!psphotStackMatchPSFsReadout (config, view, i)) {
    1713            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
    1814            return false;
    … …  
    2319
    2420// convolve the image to match desired PSF
    25 bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     21bool psphotStackMatchPSFsReadout (pmConfig *config, const pmFPAview *view, int index) {
    2622
    2723    bool status;
    … …  
    2925    float NSIGMA_PEAK = 25.0;
    3026    int NMAX = 0;
    31 
    32     Raw;
    33     Cnv;
    3427
    3528    // find the currently selected readout
    … …  
    5043    psAssert (maskVal, "missing mask value?");
    5144
    52     psphotStackMatch();
     45    /***** set up recipe options *****/
    5346
    54     return true;
    55 }
    56 
    57 #define ARRAY_BUFFER 16                 // Number to add to array at a time
    58 #define MAG_IGNORE 50                   // Ignore magnitudes fainter than this --- they're not real!
    59 #define FAKE_SIZE 1                     // Size of fake convolution kernel
    60 #define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
    61                      PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources
    62 #define NOISE_FRACTION 0.01             // Set minimum flux to this fraction of noise
    63 #define COVAR_FRAC 0.01                 // Truncation fraction for covariance matrix
    64 
    65 // #define TESTING                         // Enable debugging output
    66 
    67 #ifdef TESTING
    68 // Read a FITS image
    69 static bool readImage(psImage **target, // Target for image
    70                       const char *name, // Name of FITS file
    71                       const pmConfig *config // Configuration
    72     )
    73 {
    74     psString resolved = pmConfigConvertFilename(name, config, false, false); // Resolved filename
    75     psFits *fits = psFitsOpen(resolved, "r");
    76     psFree(resolved);
    77     if (!fits) {
    78         psError(PPSTACK_ERR_IO, false, "Unable to open previously produced image: %s", name);
    79         return false;
    80     }
    81     psImage *image = psFitsReadImage(fits, psRegionSet(0,0,0,0), 0); // Image of interest
    82     if (!image) {
    83         psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image: %s", name);
    84         psFitsClose(fits);
    85         return false;
    86     }
    87     psFitsClose(fits);
    88 
    89     psFree(*target);
    90     *target = image;
    91 
    92     return true;
    93 }
    94 #endif
    95 
    96 // Get coordinates from a source
    97 static void coordsFromSource(float *x, float *y, const pmSource *source)
    98 {
    99     assert(x && y);
    100     assert(source);
    101 
    102     if (source->modelPSF) {
    103         *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
    104         *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
    105     } else {
    106         *x = source->peak->xf;
    107         *y = source->peak->yf;
    108     }
    109     return;
    110 }
    111 
    112 static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
    113                                    int exclusion // Exclusion zone, pixels
    114     )
    115 {
    116     psAssert(sources && sources->n > 0, "Require array of sources");
    117     if (exclusion <= 0) {
    118         return psMemIncrRefCounter(sources);
    119     }
    120 
    121     int num = sources->n;               // Number of sources
    122     psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates
    123     int numGood = 0;                    // Number of good sources
    124     for (int i = 0; i < num; i++) {
    125         pmSource *source = sources->data[i]; // Source of interest
    126         if (!source) {
    127             continue;
    128         }
    129         coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
    130         numGood++;
    131     }
    132     x->n = y->n = numGood;
    133 
    134     psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree
    135 
    136     psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources
    137     psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source
    138     int numFiltered = 0;                // Number of filtered sources
    139     for (int i = 0; i < num; i++) {
    140         pmSource *source = sources->data[i]; // Source of interest
    141         if (!source) {
    142             continue;
    143         }
    144         float xSource, ySource;         // Coordinates of source
    145         coordsFromSource(&xSource, &ySource, source);
    146 
    147         coords->data.F64[0] = xSource;
    148         coords->data.F64[1] = ySource;
    149 
    150         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
    151         psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone",
    152                 coords->data.F64[0], coords->data.F64[1], numWithin);
    153         if (numWithin == 1) {
    154             // Only itself inside the exclusion zone
    155             filtered = psArrayAdd(filtered, filtered->n, source);
    156         } else {
    157             numFiltered++;
    158         }
    159     }
    160     psFree(coords);
    161     psFree(tree);
    162     psFree(x);
    163     psFree(y);
    164 
    165     psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood);
    166 
    167     return filtered;
    168 }
    169 
    170 // Add background into the fake image
    171 // Based on ppSubBackground()
    172 static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
    173                                      const pmConfig *config // Configuration
    174     )
    175 {
    176     psAssert(ro && ro->image, "Need readout image");
    177     psAssert(config, "Need configuration");
    178 
    179     psImage *image = ro->image;         // Image of interest
    180     int numCols = image->numCols, numRows = image->numRows; // Size of image
    181 
    182     psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE);
    183     psAssert(ppStackRecipe, "Need PPSTACK recipe");
    184     psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
    185     psAssert(psphotRecipe, "Need PSPHOT recipe");
    186 
    187     psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad
    188     psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    189 
    190     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    191     psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
    192 
    193     psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model
    194     psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis,
    195                                                   "PSPHOT.BACKGROUND.BINNING"); // Binning for model
    196     psAssert(binning, "Need binning parameters");
    197     psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
    198     if (!psImageUnbin(unbinned, binned, binning)) {
    199         psError(PPSTACK_ERR_DATA, false, "Unable to unbin background model");
    200         psFree(binned);
    201         psFree(unbinned);
    202         return NULL;
    203     }
    204     psFree(binned);
    205 
    206     return unbinned;
    207 }
    208 
    209 // Renormalise a readout's variance map
    210 bool stackRenormaliseReadout(const pmConfig *config, // Configuration
    211                              pmReadout *readout      // Readout to renormalise
    212     )
    213 {
    214 #if 1
    215     bool mdok; // Status of metadata lookups
    216 
    217     psMetadata *recipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE); // Recipe for ppStack
    218     psAssert(recipe, "Need PPSTACK recipe");
    219 
    220     if (!psMetadataLookupBool(&mdok, recipe, "RENORM")) return true;
    221 
    222     int num = psMetadataLookupS32(&mdok, recipe, "RENORM.NUM");
    223     if (!mdok) {
    224         psError(PPSTACK_ERR_CONFIG, true, "RENORM.NUM is not set in the recipe");
    225         return false;
    226     }
    227     float minValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MIN");
    228     if (!mdok) {
    229         psError(PPSTACK_ERR_CONFIG, true, "RENORM.MIN is not set in the recipe");
    230         return false;
    231     }
    232     float maxValid = psMetadataLookupF32(&mdok, recipe, "RENORM.MAX");
    233     if (!mdok) {
    234         psError(PPSTACK_ERR_CONFIG, true, "RENORM.MAX is not set in the recipe");
    235         return false;
    236     }
    237 
    238     psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask
    239 
    240     psImageCovarianceTransfer(readout->variance, readout->covariance);
    241     return pmReadoutVarianceRenormalise(readout, maskBad, num, minValid, maxValid);
    242 #else
    243     return true;
    244 #endif
    245 }
    246 
    247 bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
    248 {
    249     assert(readout);
    250     assert(options);
    251     assert(config);
    252 
    253     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    254     psAssert(recipe, "We've thrown an error on this before.");
    255 
    256     float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
     47    psMetadata *stackRecipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
     48    psAssert(stackRecipe, "We've thrown an error on this before.");
    25749
    25850    // Look up appropriate values from the ppSub recipe
    259     psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
    260     int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size
     51    psMetadata *subRecipe = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe
     52    psAssert(subRecipe, "recipe missing");
    26153
    262     psString maskValStr = psMetadataLookupStr(NULL, ppsub, "MASK.VAL"); // Name of bits to mask going in
     54    int size = psMetadataLookupS32(NULL, subRecipe, "KERNEL.SIZE"); // Kernel half-size
     55
     56    float deconvLimit = psMetadataLookupF32(NULL, stackRecipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
     57
     58    psString maskValStr = psMetadataLookupStr(NULL, subRecipe, "MASK.VAL"); // Name of bits to mask going in
    26359    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch
    264     psString maskPoorStr = psMetadataLookupStr(NULL, recipe, "MASK.POOR"); // Name of bits to mask for poor
     60    psString maskPoorStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.POOR"); // Name of bits to mask for poor
    26561    psImageMaskType maskPoor = pmConfigMaskGet(maskPoorStr, config); // Bits to mask for poor pixels
    266     psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad
     62    psString maskBadStr = psMetadataLookupStr(NULL, stackRecipe, "MASK.BAD"); // Name of bits to mask for bad
    26763    psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
    26864
    26965    bool mdok;                          // Status of MD lookup
    270     float penalty = psMetadataLookupF32(NULL, ppsub, "PENALTY"); // Penalty for wideness
     66    float penalty = psMetadataLookupF32(NULL, subRecipe, "PENALTY"); // Penalty for wideness
    27167    int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads
    27268
    … …  
    27672    }
    27773
    278     // Match the PSF
     74    // Image Matching (PSFs or just flux)
    27975    if (options->convolve) {
     76      // Full match of PSFs
    28077        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
    281 #ifdef TESTING
    282         // This is a hack to use the temporary convolved images and kernel generated previously.
    283         // This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
    284         // It implicitly assumes the output root name is the same between invocations.
    28578
    28679        // Read previously produced kernel
    28780        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
    288             pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
    289             psAssert(file, "Require file");
     81          loadKernel();
     82        } else {
     83          matchKernel();
     84        } // !DEBUG.STACK
    29085
    291             pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
    292             view->chip = view->cell = view->readout = 0;
    293             psString filename = pmFPAfileNameFromRule(file->filerule, file, view); // Filename of interest
     86        saveMatchData();
    29487
    295             // Read convolution kernel
    296             psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
    297             psFree(filename);
    298             psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
    299             psFree(resolved);
    300             if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
    301                 psError(PPSTACK_ERR_IO, false, "Unable to read previously produced kernel");
    302                 psFitsClose(fits);
    303                 return false;
    304             }
    305             psFitsClose(fits);
     88        saveChiSquare();
    30689
    307             if (!readImage(&readout->image, options->convImages->data[index], config) ||
    308                 !readImage(&readout->mask, options->convMasks->data[index], config) ||
    309                 !readImage(&readout->variance, options->convVariances->data[index], config)) {
    310                 psError(PPSTACK_ERR_IO, false, "Unable to read previously produced image.");
    311                 return false;
    312             }
    313 
    314             psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
    315                                                    PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
    316             pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
    317                                                                 PM_SUBTRACTION_ANALYSIS_KERNEL);
    318 
    319             pmSubtractionAnalysis(conv->analysis, NULL, kernels, region,
    320                                   readout->image->numCols, readout->image->numRows);
    321 
    322             psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    323             bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
    324             psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
    325             psImageCovarianceSetThreads(oldThreads);
    326             psFree(readout->covariance);
    327             readout->covariance = covar;
    328             psFree(kernel);
    329         } else {
    330 #endif
    331 
    332             // Normal operations here
    333             psAssert(options->psf, "Require target PSF");
    334             psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
    335 
    336             int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
    337             float regionSize = psMetadataLookupF32(NULL, ppsub, "REGION.SIZE"); // Size of iso-kernel regs
    338             float spacing = psMetadataLookupF32(NULL, ppsub, "STAMP.SPACING"); // Typical stamp spacing
    339             int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size
    340             float threshold = psMetadataLookupF32(NULL, ppsub, "STAMP.THRESHOLD"); // Threshold for stmps
    341             int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches
    342             int iter = psMetadataLookupS32(NULL, ppsub, "ITER"); // Rejection iterations
    343             float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
    344             float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel
    345             float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw
    346             float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images
    347             float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky
    348             float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation
    349 
    350             const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
    351             pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
    352             psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
    353             psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
    354             int inner = psMetadataLookupS32(NULL, ppsub, "INNER"); // Inner radius
    355             int ringsOrder = psMetadataLookupS32(NULL, ppsub, "RINGS.ORDER"); // RINGS polynomial order
    356             int binning = psMetadataLookupS32(NULL, ppsub, "SPAM.BINNING"); // Binning for SPAM kernel
    357             float badFrac = psMetadataLookupF32(NULL, ppsub, "BADFRAC"); // Maximum bad fraction
    358             bool optimum = psMetadataLookupBool(&mdok, ppsub, "OPTIMUM"); // Derive optimum parameters?
    359             float optMin = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MIN"); // Minimum width for search
    360             float optMax = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.MAX"); // Maximum width for search
    361             float optStep = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.STEP"); // Step for search
    362             float optThresh = psMetadataLookupF32(&mdok, ppsub, "OPTIMUM.TOL"); // Tolerance for search
    363             int optOrder = psMetadataLookupS32(&mdok, ppsub, "OPTIMUM.ORDER"); // Order for search
    364             float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor"
    365 
    366             bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE");        // Scale kernel parameters?
    367             float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling
    368             float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling
    369             float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling
    370             if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) {
    371                 psError(PPSTACK_ERR_CONFIG, false,
    372                         "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.",
    373                         scaleRef, scaleMin, scaleMax);
    374                 return false;
    375             }
    376 
    377 
    378             // These values are specified specifically for stacking
    379             const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename
    380 
    381             psVector *optWidths = NULL;         // Vector with FWHMs for optimum search
    382             if (optimum) {
    383                 optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32);
    384             }
    385 
    386             pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
    387 
    388             psStats *bg = psStatsAlloc(PS_STAT_ROBUST_STDEV); // Statistics for background
    389             psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    390             if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    391                 psError(PPSTACK_ERR_DATA, false, "Can't measure background for image.");
    392                 psFree(fake);
    393                 psFree(optWidths);
    394                 psFree(conv);
    395                 psFree(bg);
    396                 psFree(rng);
    397                 return false;
    398             }
    399             float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image
    400             psFree(rng);
    401             psFree(bg);
    402 
    403             // For the sake of stamps, remove nearby sources
    404             psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
    405                                                        footprint); // Filtered list of sources
    406 
    407             bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
    408             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
    409                                           stampSources, SOURCE_MASK, NULL, NULL, options->psf,
    410                                           minFlux, footprint + size, false, true)) {
    411                 psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF.");
    412                 psFree(fake);
    413                 psFree(optWidths);
    414                 psFree(conv);
    415                 return false;
    416             }
    417             pmReadoutFakeThreads(oldThreads);
    418 
    419             fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);
    420 
    421 #if 1
    422             // Add the background into the target image
    423             psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
    424             psBinaryOp(fake->image, fake->image, "+", bgImage);
    425             psFree(bgImage);
    426 #endif
    427 
    428 #ifdef TESTING
    429             {
    430                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    431                 psString name = NULL;
    432                 psStringAppend(&name, "fake_%03d.fits", index);
    433                 pmStackVisualPlotTestImage(fake->image, name);
    434                 psFits *fits = psFitsOpen(name, "w");
    435                 psFree(name);
    436                 psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    437                 psFitsClose(fits);
    438             }
    439             {
    440                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    441                 psString name = NULL;
    442                 psStringAppend(&name, "real_%03d.fits", index);
    443                 pmStackVisualPlotTestImage(readout->image, name);
    444                 psFits *fits = psFitsOpen(name, "w");
    445                 psFree(name);
    446                 psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    447                 psFitsClose(fits);
    448             }
    449 #endif
    450 
    451             if (threads > 0) {
    452                 pmSubtractionThreadsInit();
    453             }
    454 
    455             // Do the image matching
    456             pmSubtractionKernels *kernel = psMetadataLookupPtr(&mdok, readout->analysis,
    457                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
    458             if (kernel) {
    459                 if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis,
    460                                                stride, kernelError, covarFrac, maskVal, maskBad, maskPoor,
    461                                                poorFrac, badFrac)) {
    462                     psError(psErrorCodeLast(), false, "Unable to convolve images.");
    463                     psFree(fake);
    464                     psFree(optWidths);
    465                     psFree(stampSources);
    466                     psFree(conv);
    467                     if (threads > 0) {
    468                         pmSubtractionThreadsFinalize();
    469                     }
    470                     return false;
    471                 }
    472             } else {
    473                 // Scale the input parameters
    474                 psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths
    475                 if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy,
    476                                                        options->inputSeeing->data.F32[index],
    477                                                        options->targetSeeing, scaleRef, scaleMin, scaleMax)) {
    478                     psError(psErrorCodeLast(), false, "Unable to scale kernel parameters");
    479                     psFree(fake);
    480                     psFree(optWidths);
    481                     psFree(stampSources);
    482                     psFree(conv);
    483                     psFree(widthsCopy);
    484                     if (threads > 0) {
    485                         pmSubtractionThreadsFinalize();
    486                     }
    487                     return false;
    488                 }
    489 
    490                 if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing,
    491                                         threshold, stampSources, stampsName, type, size, order, widthsCopy,
    492                                         orders, inner, ringsOrder, binning, penalty,
    493                                         optimum, optWidths, optOrder, optThresh, iter, rej, normFrac,
    494                                         sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor,
    495                                         poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) {
    496                     psError(psErrorCodeLast(), false, "Unable to match images.");
    497                     psFree(fake);
    498                     psFree(optWidths);
    499                     psFree(stampSources);
    500                     psFree(conv);
    501                     psFree(widthsCopy);
    502                     if (threads > 0) {
    503                         pmSubtractionThreadsFinalize();
    504                     }
    505                     return false;
    506                 }
    507                 psFree(widthsCopy);
    508             }
    509 
    510 
    511 #ifdef TESTING
    512             {
    513                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    514                 psString name = NULL;
    515                 psStringAppend(&name, "conv_%03d.fits", index);
    516                 pmStackVisualPlotTestImage(conv->image, name);
    517                 psFits *fits = psFitsOpen(name, "w");
    518                 psFree(name);
    519                 psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL);
    520                 psFitsClose(fits);
    521             }
    522             {
    523                 pmHDU *hdu = pmHDUFromCell(readout->parent);
    524                 psString name = NULL;
    525                 psStringAppend(&name, "diff_%03d.fits", index);
    526                 pmStackVisualPlotTestImage(fake->image, name);
    527                 psFits *fits = psFitsOpen(name, "w");
    528                 psFree(name);
    529                 psBinaryOp(fake->image, conv->image, "-", fake->image);
    530                 psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    531                 psFitsClose(fits);
    532             }
    533 #endif
    534 
    535             psFree(fake);
    536             psFree(optWidths);
    537             psFree(stampSources);
    538 
    539             if (threads > 0) {
    540                 pmSubtractionThreadsFinalize();
    541             }
    542 
    543             // Replace original images with convolved
    544             psFree(readout->image);
    545             psFree(readout->mask);
    546             psFree(readout->variance);
    547             psFree(readout->covariance);
    548             readout->image  = psMemIncrRefCounter(conv->image);
    549             readout->mask   = psMemIncrRefCounter(conv->mask);
    550             readout->variance = psMemIncrRefCounter(conv->variance);
    551             readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC);
    552 #ifdef TESTING
    553         }
    554 #endif
    555 
    556         // Extract the regions and solutions used in the image matching
    557         // This stops them from being freed when we iterate back up the FPA
    558         psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
    559         {
    560             psString regex = NULL;          // Regular expression
    561             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
    562             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    563             psFree(regex);
    564             psMetadataItem *item = NULL;// Item from iteration
    565             while ((item = psMetadataGetAndIncrement(iter))) {
    566                 assert(item->type == PS_DATA_REGION);
    567                 regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
    568             }
    569             psFree(iter);
    570         }
    571         psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
    572         {
    573             psString regex = NULL;          // Regular expression
    574             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    575             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    576             psFree(regex);
    577             psMetadataItem *item = NULL;// Item from iteration
    578             while ((item = psMetadataGetAndIncrement(iter))) {
    579                 assert(item->type == PS_DATA_UNKNOWN);
    580                 pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
    581                 kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
    582             }
    583             psFree(iter);
    584         }
    585         psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
    586 
    587         // Record chi^2
    588         {
    589             double sum = 0.0;           // Sum of chi^2
    590             int num = 0;                // Number of measurements of chi^2
    591             psString regex = NULL;      // Regular expression
    592             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    593             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    594             psFree(regex);
    595             psMetadataItem *item = NULL;// Item from iteration
    596             while ((item = psMetadataGetAndIncrement(iter))) {
    597                 assert(item->type == PS_DATA_UNKNOWN);
    598                 pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
    599                 sum += kernels->mean;
    600                 num++;
    601             }
    602             psFree(iter);
    603             options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
    604         }
    605 
    606         // Kernel normalisation
    607         {
    608             double sum = 0.0;           // Sum of chi^2
    609             int num = 0;                // Number of measurements of chi^2
    610             psString regex = NULL;      // Regular expression
    611             psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_NORM);
    612             psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
    613             psFree(regex);
    614             psMetadataItem *item = NULL;// Item from iteration
    615             while ((item = psMetadataGetAndIncrement(iter))) {
    616                 assert(item->type == PS_TYPE_F32);
    617                 float norm = item->data.F32; // Normalisation
    618                 sum += norm;
    619                 num++;
    620             }
    621             psFree(iter);
    622             float conv = sum/num;       // Mean normalisation from convolution
    623             float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars
    624             float renorm =  stars / conv; // Renormalisation to apply
    625             psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n",
    626                      index, renorm, conv, stars);
    627             psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32));
    628             psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32));
    629         }
     90        renormKernel();
    63091
    63192        // Reject image completely if the maximum deconvolution fraction exceeds the limit
    632         float deconv = psMetadataLookupF32(NULL, conv->analysis,
    633                                            PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
     93        float deconv = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
    63494        if (deconv > deconvLimit) {
    635             psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n",
    636                       deconv, deconvLimit, index);
     95            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", deconv, deconvLimit, index);
    63796            psFree(conv);
    63897            return NULL;
    … …  
    643102        psFree(conv);
    644103    } else {
    645         // Match the normalisation
     104        // only match the flux
    646105        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
    647106        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    … …  
    671130
    672131    // Measure the variance level for the weighting
    673     if (psMetadataLookupBool(NULL, recipe, "WEIGHTS")) {
     132    if (psMetadataLookupBool(NULL, stackRecipe, "WEIGHTS")) {
    674133        if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
    675134            psError(PPSTACK_ERR_DATA, false, "Can't measure mean variance for image.");
    … …  
    678137            return false;
    679138        }
    680         options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
    681                                                       psImageCovarianceFactor(readout->covariance));
     139        options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * psImageCovarianceFactor(readout->covariance));
    682140    } else {
    683141        options->weightings->data.F32[index] = 1.0;
    … …  
    689147    psFree(bg);
    690148
    691 #ifdef TESTING
    692     {
    693         pmHDU *hdu = pmHDUFromCell(readout->parent);
    694         psString name = NULL;
    695         psStringAppend(&name, "convolved_%03d.fits", index);
    696         pmStackVisualPlotTestImage(readout->image, name);
    697         psFits *fits = psFitsOpen(name, "w");
    698         psFree(name);
    699         psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    700         psFitsClose(fits);
    701     }
    702 #endif
     149    dumpImage3();
    703150
    704151    return true;
Note: See TracChangeset for help on using the changeset viewer.