IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 27, 2009, 11:49:17 AM (17 years ago)
Author:
Paul Price
Message:

Allow convolution of images to be optional when stacking.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStackMatch.c

    r23379 r23573  
    163163
    164164
    165 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
    166                   psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config)
     165bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config)
    167166{
    168167    assert(readout);
    169     assert(regions && !*regions);
    170     assert(kernels && !*kernels);
     168    assert(options);
    171169    assert(config);
    172     *weighting = 0.0;
    173170
    174171    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
    … …  
    197194    }
    198195
    199     pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily
    200 
    201     static int numInput = -1;            // Index of input file
    202     numInput++;
     196    // Match the PSF
     197    if (options->convolve) {
     198        pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily
    203199#ifdef TESTING
    204     // Read previously produced kernel
    205     if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
    206         const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
    207         assert(outName);
    208         // Read convolution kernel
    209         psString filename = NULL;   // Output filename
    210         psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
    211         psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
    212         psFree(filename);
    213         psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
    214         psFree(resolved);
    215         if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
    216             psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
     200        // This is a hack to use the temporary convolved images and kernel generated previously.
     201        // This makes the 'matching' operation much faster, allowing debugging of the stack process easier.
     202        // It implicitly assumes the output root name is the same between invocations.
     203
     204        // Read previously produced kernel
     205        if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) {
     206            pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index);
     207            psAssert(file, "Require file");
     208
     209            pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest
     210            view->chip = view->cell = view->readout = 0;
     211            psString filename = pmFPAfileNameFromRule(filerule->rule, file, view); // Filename of interest
     212
     213            // Read convolution kernel
     214            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
     215            psFree(filename);
     216            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
     217            psFree(resolved);
     218            if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) {
     219                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
     220                psFitsClose(fits);
     221                return false;
     222            }
    217223            psFitsClose(fits);
    218             return false;
    219         }
    220         psFitsClose(fits);
    221 
    222         // Add in variance factor
    223         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
    224                                                             PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
    225         float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
    226         psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
    227         if (!isfinite(vf)) {
    228             vf = 1.0;
    229         }
    230         if (isfinite(vfItem->data.F32)) {
    231             vfItem->data.F32 *= vf;
    232         } else {
    233             vfItem->data.F32 = vf;
    234         }
    235 
    236         // Read image, mask, variance
    237         const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image
    238         const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask
    239         const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for variance map
    240         psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images
    241         psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage);
    242         psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask);
    243         psStringAppend(&varianceName, "%s.%d.%s", outName, numInput, tempVariance);
    244 
    245         if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) ||
    246             !readImage(&readout->variance, varianceName, config)) {
    247             psError(PS_ERR_IO, false, "Unable to read previously produced image.");
     224
     225            // Add in variance factor
     226            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis,
     227                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
     228            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
     229            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
     230            if (!isfinite(vf)) {
     231                vf = 1.0;
     232            }
     233            if (isfinite(vfItem->data.F32)) {
     234                vfItem->data.F32 *= vf;
     235            } else {
     236                vfItem->data.F32 = vf;
     237            }
     238
     239            if (!readImage(&readout->image, options->imageNames->data[index], config) ||
     240                !readImage(&readout->mask, options->maskNames->data[index], config) ||
     241                !readImage(&readout->variance, options->varianceNames->data[index], config)) {
     242                psError(PS_ERR_IO, false, "Unable to read previously produced image.");
     243                psFree(imageName);
     244                psFree(maskName);
     245                psFree(varianceName);
     246                return false;
     247            }
    248248            psFree(imageName);
    249249            psFree(maskName);
    250250            psFree(varianceName);
    251             return false;
    252         }
    253         psFree(imageName);
    254         psFree(maskName);
    255         psFree(varianceName);
    256 
    257         psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
    258                                                PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
    259 
    260         pmSubtractionAnalysis(readout->analysis, kernels, region,
    261                               readout->image->numCols, readout->image->numRows);
    262 
    263         psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    264         psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix
    265         psFree(readout->covariance);
    266         readout->covariance = covar;
    267         psFree(kernel);
    268 
    269     } else {
    270 #endif
    271 
    272         // Normal operations here
    273         if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) {
    274             assert(psf);
    275             assert(sources);
     251
     252            psRegion *region = psMetadataLookupPtr(NULL, conv->analysis,
     253                                                   PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
     254
     255            pmSubtractionAnalysis(readout->analysis, kernels, region,
     256                                  readout->image->numCols, readout->image->numRows);
     257
     258            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     259            psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
     260            psFree(readout->covariance);
     261            readout->covariance = covar;
     262            psFree(kernel);
     263        } else {
     264#endif
     265
     266            // Normal operations here
     267            psAssert(options->psf, "Require target PSF");
     268            psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list");
    276269
    277270            int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order
    … …  
    284277            float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold
    285278            float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel
    286             pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(
    287                 psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type
     279            const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type
     280            pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type
    288281            psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths
    289282            psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders
    … …  
    308301            }
    309302
    310 #if 0
    311             // Testing the normalisation of the fake image
    312             {
    313                 pmReadout *test = pmReadoutAlloc(NULL); // Test readout
    314                 psArray *sources = psArrayAlloc(1);     // Array of sources
    315                 pmSource *source = pmSourceAlloc();     // Source
    316                 sources->data[0] = source;
    317                 source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
    318                 source->psfMag = -13.0;
    319                 pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
    320                 float sum = 0.0;
    321                 for (int y = 0; y < test->image->numRows; y++) {
    322                     for (int x = 0; x < test->image->numCols; x++) {
    323                         sum += test->image->data.F32[y][x];
    324                     }
    325                 }
    326                 fprintf(stderr, "Photometry: %f --> %f = -13.0 ???\n", sum, -2.5*log10(sum));
    327 
    328                 psFits *fits = psFitsOpen("testphot.fits", "w");
    329                 psFitsWriteImage(fits, NULL, test->image, 0, NULL);
    330                 psFitsClose(fits);
    331                 exit(0);
    332             }
    333 #endif
    334 
    335303            pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF
    336304
    337305            // For the sake of stamps, remove nearby sources
    338             psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources
     306            psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index],
     307                                                       footprint); // Filtered list of sources
    339308
    340309            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
    341                                           stampSources, NULL, NULL, psf, NAN, footprint + size,
     310                                          stampSources, NULL, NULL, options->psf, NAN, footprint + size,
    342311                                          false, true)) {
    343312                psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
    344313                psFree(fake);
    345314                psFree(optWidths);
    346                 psFree(output);
     315                psFree(conv);
    347316                return false;
    348317            }
    … …  
    357326                pmHDU *hdu = pmHDUFromCell(readout->parent);
    358327                psString name = NULL;
    359                 psStringAppend(&name, "fake_%03d.fits", numInput);
     328                psStringAppend(&name, "fake_%03d.fits", index);
    360329                pmStackVisualPlotTestImage(fake->image, name);
    361330                psFits *fits = psFitsOpen(name, "w");
    … …  
    367336                pmHDU *hdu = pmHDUFromCell(readout->parent);
    368337                psString name = NULL;
    369                 psStringAppend(&name, "real_%03d.fits", numInput);
     338                psStringAppend(&name, "real_%03d.fits", index);
    370339                pmStackVisualPlotTestImage(readout->image, name);
    371340                psFits *fits = psFitsOpen(name, "w");
    … …  
    384353                                                               PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel
    385354            if (kernel) {
    386                 if (!pmSubtractionMatchPrecalc(output, NULL, readout, fake, readout->analysis,
     355                if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis,
    387356                                               stride, sysError, maskVal, maskBad, maskPoor,
    388357                                               poorFrac, badFrac)) {
    … …  
    391360                    psFree(optWidths);
    392361                    psFree(stampSources);
    393                     psFree(output);
     362                    psFree(conv);
    394363                    return false;
    395364                }
    396365            } else {
    397                 if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing,
     366                if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing,
    398367                                        threshold, stampSources, stampsName, type, size, order, widths,
    399368                                        orders, inner, ringsOrder, binning, penalty,
    … …  
    405374                    psFree(optWidths);
    406375                    psFree(stampSources);
    407                     psFree(output);
     376                    psFree(conv);
    408377                    return false;
    409378                }
    … …  
    414383                pmHDU *hdu = pmHDUFromCell(readout->parent);
    415384                psString name = NULL;
    416                 psStringAppend(&name, "conv_%03d.fits", numInput);
    417                 pmStackVisualPlotTestImage(output->image, name);
     385                psStringAppend(&name, "conv_%03d.fits", index);
     386                pmStackVisualPlotTestImage(conv->image, name);
    418387                psFits *fits = psFitsOpen(name, "w");
    419388                psFree(name);
    420                 psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
     389                psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL);
    421390                psFitsClose(fits);
    422391            }
    … …  
    424393                pmHDU *hdu = pmHDUFromCell(readout->parent);
    425394                psString name = NULL;
    426                 psStringAppend(&name, "diff_%03d.fits", numInput);
     395                psStringAppend(&name, "diff_%03d.fits", index);
    427396                pmStackVisualPlotTestImage(fake->image, name);
    428397                psFits *fits = psFitsOpen(name, "w");
    429398                psFree(name);
    430                 psBinaryOp(fake->image, output->image, "-", fake->image);
     399                psBinaryOp(fake->image, conv->image, "-", fake->image);
    431400                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
    432401                psFitsClose(fits);
    … …  
    444413            // Set the variance factor
    445414            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
    446             float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
     415            float vf = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1);
    447416            if (!isfinite(vf)) {
    448417                vf = 1.0;
    … …  
    459428            psFree(readout->variance);
    460429            psFree(readout->covariance);
    461             readout->image  = psMemIncrRefCounter(output->image);
    462             readout->mask   = psMemIncrRefCounter(output->mask);
    463             readout->variance = psMemIncrRefCounter(output->variance);
    464             readout->covariance = psImageCovarianceTruncate(output->covariance, COVAR_FRAC);
    465         } else {
    466             // Fake the convolution
    467             psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1);
    468             psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION,
    469                              PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region);
    470             psFree(region);
    471             pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty,
    472                                                                      PM_SUBTRACTION_MODE_1);
    473             // Set solution to delta function
    474             kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64);
    475             psVectorInit(kernels->solution1, 0.0);
    476             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
    477             kernels->solution1->data.F64[normIndex] = 1.0;
    478             psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL,
    479                              PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels);
    480             psFree(kernels);
    481         }
    482 
     430            readout->image  = psMemIncrRefCounter(conv->image);
     431            readout->mask   = psMemIncrRefCounter(conv->mask);
     432            readout->variance = psMemIncrRefCounter(conv->variance);
     433            readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC);
    483434#ifdef TESTING
    484         // Write convolution kernel
     435        }
     436#endif
     437
     438        // Extract the regions and solutions used in the image matching
     439        // This stops them from being freed when we iterate back up the FPA
     440        psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions
    485441        {
    486             const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
    487             assert(outName);
    488 
    489             psString filename = NULL;   // Output filename
    490             psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
    491             psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename
    492             psFree(filename);
    493             psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel
    494             psFree(resolved);
    495             pmReadoutWriteSubtractionKernels(output, fits);
    496             psFitsClose(fits);
    497         }
    498     }
    499 #endif
    500 
    501     readout->analysis = psMetadataCopy(readout->analysis, output->analysis);
    502 
    503 // Extract the regions and solutions used in the image matching
    504 // This stops them from being freed when we iterate back up the FPA
    505     *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions
    506     {
    507         psString regex = NULL;          // Regular expression
    508         psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
    509         psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
    510         psFree(regex);
    511         psMetadataItem *item = NULL;// Item from iteration
    512         while ((item = psMetadataGetAndIncrement(iter))) {
    513             assert(item->type == PS_DATA_REGION);
    514             *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V);
    515         }
    516         psFree(iter);
    517     }
    518     *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
    519     {
    520         psString regex = NULL;          // Regular expression
    521         psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    522         psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
    523         psFree(regex);
    524         psMetadataItem *item = NULL;// Item from iteration
    525         while ((item = psMetadataGetAndIncrement(iter))) {
    526             assert(item->type == PS_DATA_UNKNOWN);
    527             // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
    528             // images by scans.
    529             pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
    530             kernel->numCols = readout->image->numCols;
    531             kernel->numRows = readout->image->numRows;
    532 
    533             *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel);
    534         }
    535         psFree(iter);
    536     }
    537     assert((*regions)->n == (*kernels)->n);
    538 
    539     // Record chi^2
    540     {
    541         *chi2 = 0.0;
    542         int num = 0;                    // Number of measurements of chi^2
    543         psString regex = NULL;          // Regular expression
    544         psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
    545         psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator
    546         psFree(regex);
    547         psMetadataItem *item = NULL;// Item from iteration
    548         while ((item = psMetadataGetAndIncrement(iter))) {
    549             assert(item->type == PS_DATA_UNKNOWN);
    550             pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
    551             *chi2 += kernels->mean;
    552             num++;
    553         }
    554         psFree(iter);
    555         *chi2 /= psImageCovarianceFactor(readout->covariance) * num;
    556     }
    557 
    558     // Reject image completely if the maximum deconvolution fraction exceeds the limit
    559     float deconv = psMetadataLookupF32(NULL, output->analysis,
    560                                        PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
    561     if (deconv > deconvLimit) {
    562         psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
    563                   deconv, deconvLimit);
    564         psFree(output);
    565         return NULL;
     442            psString regex = NULL;          // Regular expression
     443            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION);
     444            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     445            psFree(regex);
     446            psMetadataItem *item = NULL;// Item from iteration
     447            while ((item = psMetadataGetAndIncrement(iter))) {
     448                assert(item->type == PS_DATA_REGION);
     449                regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V);
     450            }
     451            psFree(iter);
     452        }
     453        psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels
     454        {
     455            psString regex = NULL;          // Regular expression
     456            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     457            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     458            psFree(regex);
     459            psMetadataItem *item = NULL;// Item from iteration
     460            while ((item = psMetadataGetAndIncrement(iter))) {
     461                assert(item->type == PS_DATA_UNKNOWN);
     462                // Set the normalisation dimensions, since these will be otherwise unavailable when reading
     463                // the images by scans.
     464                pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
     465                kernel->numCols = readout->image->numCols;
     466                kernel->numRows = readout->image->numRows;
     467
     468                kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
     469            }
     470            psFree(iter);
     471        }
     472        psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match");
     473
     474        // Record chi^2
     475        {
     476            double sum = 0.0;           // Sum of chi^2
     477            int num = 0;                // Number of measurements of chi^2
     478            psString regex = NULL;      // Regular expression
     479            psStringAppend(&regex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL);
     480            psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex);
     481            psFree(regex);
     482            psMetadataItem *item = NULL;// Item from iteration
     483            while ((item = psMetadataGetAndIncrement(iter))) {
     484                assert(item->type == PS_DATA_UNKNOWN);
     485                pmSubtractionKernels *kernels = item->data.V; // Convolution kernels
     486                sum += kernels->mean;
     487                num++;
     488            }
     489            psFree(iter);
     490            options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num);
     491        }
     492
     493        // Reject image completely if the maximum deconvolution fraction exceeds the limit
     494        float deconv = psMetadataLookupF32(NULL, conv->analysis,
     495                                           PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction
     496        if (deconv > deconvLimit) {
     497            psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
     498                      deconv, deconvLimit);
     499            psFree(conv);
     500            return NULL;
     501        }
     502
     503        readout->analysis = psMetadataCopy(readout->analysis, conv->analysis);
     504
     505        psFree(conv);
     506    } else {
     507        // Match the normalisation
     508        float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation
     509        psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
     510        psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32));
    566511    }
    567512
    568513    // Ensure the background value is zero
    569514    psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
     515    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
    570516    if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) {
    571517        psWarning("Can't measure background for image.");
    … …  
    581527    if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) {
    582528        psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image.");
    583         psFree(output);
     529        psFree(rng);
     530        psFree(bg);
    584531        return false;
    585532    }
    586     *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
    587                         psImageCovarianceFactor(readout->covariance));
     533    options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) *
     534                                                  psImageCovarianceFactor(readout->covariance));
    588535    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
    589                      "Weighting by 1/noise^2 for stack", *weighting);
    590 
     536                     "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]);
     537
     538    psFree(rng);
    591539    psFree(bg);
    592 
    593 #if 0
    594 #define RADIUS 10                       // Radius of photometry
    595 #define MIN_ERR 0.05                    // Minimum photometric error, mag
    596 #define MAX_MAG -13                     // Maximum magnitude for source
    597 
    598     // Ensure the normalisation is correct
    599     // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
    600     int numSources = sources->n;        // Number of sources
    601     psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
    602     psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
    603     psVectorInit(ratioMask, 0xFF);
    604     psImage *image = readout->image;    // Image of interest
    605     psImage *mask = readout->mask;      // Mask of interest
    606     int numCols = image->numCols, numRows = image->numRows; // Size of image
    607     for (int i = 0; i < numSources; i++) {
    608         pmSource *source = sources->data[i]; // Source of interest
    609         if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
    610             source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
    611             continue;
    612         }
    613 
    614         float xSrc, ySrc;              // Source coordinates
    615         coordsFromSource(&xSrc, &ySrc, source);
    616         int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
    617         int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
    618         int numPix = 0;                 // Number of pixels
    619         float sum = 0.0;                // Sum of pixels
    620         for (int y = yMin; y <= yMax; y++) {
    621             for (int x = xMin; x <= xMax; x++) {
    622                 if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
    623                     continue;
    624                 }
    625                 sum += image->data.F32[y][x];
    626                 numPix++;
    627             }
    628         }
    629         if (sum >= 0 && numPix > 0) {
    630             float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
    631             ratio->data.F32[i] = mag - source->psfMag;
    632             ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
    633         }
    634     }
    635 
    636     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
    637     if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
    638         psWarning("Unable to measure normalisation --- assuming correct.");
    639     } else {
    640         psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
    641                  stats->robustMedian, stats->robustStdev);
    642         float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
    643         psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
    644     }
    645     psFree(stats);
    646     psFree(ratio);
    647     psFree(ratioMask);
    648 #endif
    649540
    650541#ifdef TESTING
    … …  
    652543        pmHDU *hdu = pmHDUFromCell(readout->parent);
    653544        psString name = NULL;
    654         psStringAppend(&name, "convolved_%03d.fits", numInput);
    655         pmStackVisualPlotTestImage(output->image, name);
     545        psStringAppend(&name, "convolved_%03d.fits", index);
     546        pmStackVisualPlotTestImage(readout->image, name);
    656547        psFits *fits = psFitsOpen(name, "w");
    657548        psFree(name);
    658         psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
     549        psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
    659550        psFitsClose(fits);
    660551    }
    661552#endif
    662 
    663     psFree(output);
    664553
    665554    return true;
Note: See TracChangeset for help on using the changeset viewer.