IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 11, 2008, 6:12:19 PM (18 years ago)
Author:
Paul Price
Message:

Now reads temporary convolved images bit by bit. Updating the stacking process to use modern APIs (the kernel solution is now a part of the pmSubtractionKernels).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_branch_080207/ppStack/src/ppStackLoop.c

    r16382 r16404  
    1616// All files in the system
    1717static char *allFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT",
    18                             "PPSTACK.OUTCONV", "PPSTACK.OUTCONV.MASK", "PPSTACK.OUTCONV.WEIGHT",
    19                             "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT",
    2018                            "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
    2119                            "PSPHOT.PSF.SAVE", "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES",
    … …  
    2927static char *prepareFiles[] = { "PSPHOT.PSF.LOAD", "PPSTACK.SOURCES", 0 };
    3028
    31 // Input files for the combination
    32 static char *inCombineFiles[] = { "PPSTACK.INCONV", "PPSTACK.INCONV.MASK", "PPSTACK.INCONV.WEIGHT", 0 };
    33 
    34 // Output files for the combination and photometry
    35 static char *outCombineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT",
    36                                    "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
    37                                    "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
    38                                    "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
    39                                    "PSPHOT.INPUT.CMF", 0 };
     29// Files required for the convolution
     30static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.WEIGHT", 0 };
     31
     32// Output files for the combination
     33static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.WEIGHT", 0 };
     34
     35// Files for photometry
     36static char *photFiles[] = { "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL",
     37                             "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
     38                             "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID",
     39                             "PSPHOT.INPUT.CMF", 0 };
     40
     41#define CONVOLVED_ALREADY               // Already have the convolution products --- testing
    4042
    4143
    … …  
    7880}
    7981
     82#if 0
    8083// Set the data level for a list of files
    8184static void fileSetDataLevel(pmConfig *config, // Configuration
    … …  
    98101    return;
    99102}
    100 
     103#endif
    101104
    102105// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells
    … …  
    150153    return true;
    151154}
     155
     156#ifndef CONVOLVED_ALREADY
     157// Write an image to a FITS file
     158static bool writeImage(const char *name, // Name of image
     159                       psMetadata *header, // Header
     160                       const psImage *image // Image
     161                       )
     162{
     163    assert(name);
     164    assert(image);
     165
     166    psFits *fits = psFitsOpen(name, "w");
     167    if (!fits) {
     168        psError(PS_ERR_IO, false, "Unable to open FITS file to write image.");
     169        return false;
     170    }
     171    if (!psFitsWriteImage(fits, header, image, 0, NULL)) {
     172        psError(PS_ERR_IO, false, "Unable to write FITS image.");
     173        return false;
     174    }
     175    psFitsClose(fits);
     176    return true;
     177}
     178#endif
    152179
    153180
    … …  
    263290    }
    264291
     292    const char *suffix = "conv.fits";   // Suffix for convolved images; ultimately this will be from recipe
     293    const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root
     294    assert(outName);
     295    psArray *imageNames = psArrayAlloc(num);
     296    psArray *maskNames = psArrayAlloc(num);
     297    psArray *weightNames = psArrayAlloc(num);
     298    for (int i = 0; i < num; i++) {
     299        psString imageName = NULL, maskName = NULL, weightName = NULL; // Names for convolved images
     300        psStringAppend(&imageName, "%s.im-%d.%s", outName, i, suffix);
     301        psStringAppend(&maskName, "%s.mk-%d.%s", outName, i, suffix);
     302        psStringAppend(&weightName, "%s.wt-%d.%s", outName, i, suffix);
     303        imageNames->data[i] = imageName;
     304        maskNames->data[i] = maskName;
     305        weightNames->data[i] = weightName;
     306    }
    265307
    266308    // Generate convolutions and write them to disk
    267309    psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again
     310    psArray *kernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking
    268311    for (int i = 0; i < num; i++) {
    269312        pmFPAfileActivate(config->files, false, NULL);
    270         psArray *files = fileActivationSingle(config, prepareFiles, true, i);
     313        psArray *files = fileActivationSingle(config, convolveFiles, true, i);
    271314        pmFPAview *view = filesIterateDown(config);
    272315        if (!view) {
    … …  
    276319        }
    277320
    278         pmReadout *readout = pmFPAviewThisReadout(view, files->data[0]); // Input readout
     321        pmReadout *readout = pmFPAviewThisReadout(view, ((pmFPAfile*)files->data[0])->fpa); // Input readout
    279322        psFree(view);
    280323
     324#ifndef CONVOLVED_ALREADY
    281325        // Background subtraction, scaling and normalisation is performed automatically by the image matching
    282326        if (!ppStackMatch(readout, sources, targetPSF, config)) {
    … …  
    287331        }
    288332
    289 #ifdef CONVOLUTION_FILES
    290         if (readout->image) {
    291             psString name = NULL;           // Name of image
    292             psStringAppend(&name, "convolved%03d_image.fits", i);
    293             psFits *fits = psFitsOpen(name, "w");
    294             psFree(name);
    295             psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
    296             psFitsClose(fits);
    297         }
    298         if (readout->mask) {
    299             psString name = NULL;           // Name of image
    300             psStringAppend(&name, "convolved%03d_mask.fits", i);
    301             psFits *fits = psFitsOpen(name, "w");
    302             psFree(name);
    303             psFitsWriteImage(fits, NULL, readout->mask, 0, NULL);
    304             psFitsClose(fits);
    305         }
    306         if (readout->weight) {
    307             psString name = NULL;           // Name of image
    308             psStringAppend(&name, "convolved%03d_weight.fits", i);
    309             psFits *fits = psFitsOpen(name, "w");
    310             psFree(name);
    311             psFitsWriteImage(fits, NULL, readout->weight, 0, NULL);
    312             psFitsClose(fits);
    313         }
    314 #endif // CONVOLUTION_FILES
     333        // Write the temporary convolved files
     334        pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image
     335        assert(hdu);
     336        writeImage(imageNames->data[i],  hdu->header, readout->image);
     337        writeImage(maskNames->data[i],   hdu->header, readout->mask);
     338        writeImage(weightNames->data[i], hdu->header, readout->weight);
     339#endif
    315340
    316341        cells->data[i] = psMemIncrRefCounter(readout->parent);
     342        kernels->data[i] = psMemIncrRefCounter(psMetadataLookupPtr(NULL, readout->analysis,
     343                                                                   PM_SUBTRACTION_ANALYSIS_KERNEL));
    317344        filesIterateUp(config);
    318345    }
    … …  
    320347    psFree(targetPSF);
    321348
    322     // Set files to read at the readout level --- then when we iterate down, the pmFPAfileIO stuff should take
    323     // care of the basic things, and we can just worry about grabbing the readouts by chunks
    324     fileSetDataLevel(config, inCombineFiles, PM_FPA_LEVEL_CHUNK);
    325     pmFPAfileActivate(config->files, false, NULL);
    326     fileActivation(config, inCombineFiles, true);
    327     fileActivation(config, outCombineFiles, true);
    328 
     349    // Stack the convolved files
    329350    {
     351        pmFPAfileActivate(config->files, false, NULL);
     352        fileActivation(config, combineFiles, true);
     353        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY")) {
     354            fileActivation(config, photFiles, true);
     355        }
    330356        pmFPAview *view = filesIterateDown(config);
    331357        if (!view) {
    332358            psFree(cells);
     359            psFree(kernels);
    333360            return false;
    334361        }
    … …  
    341368        }
    342369        psFree(cells);
     370
     371        // FITS files
     372        psArray *imageFits  = psArrayAlloc(num);
     373        psArray *maskFits   = psArrayAlloc(num);
     374        psArray *weightFits = psArrayAlloc(num);
     375        for (int i = 0; i < num; i++) {
     376            imageFits->data[i] = psFitsOpen(imageNames->data[i], "r");
     377            maskFits->data[i] = psFitsOpen(maskNames->data[i], "r");
     378            weightFits->data[i] = psFitsOpen(weightNames->data[i], "r");
     379        }
    343380
    344381        // Read convolutions by chunks
    … …  
    349386                assert(readout);
    350387
    351                 // Files, containing the FITS handle
    352                 pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i);
    353                 pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i);
    354                 pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i);
    355 
    356                 // FITS handles
    357                 psFits *fitsImage = imageFile->fits;
    358                 psFits *fitsMask = maskFile->fits;
    359                 psFits *fitsWeight = weightFile->fits;
    360 
    361                 if (!pmReadoutReadChunk(readout, fitsImage, 0, numScans, overlap) ||
    362                     !pmReadoutReadChunkMask(readout, fitsMask, 0, numScans, overlap) ||
    363                     !pmReadoutReadChunkWeight(readout, fitsWeight, 0, numScans, overlap)) {
     388                if (!pmReadoutReadChunk(readout, imageFits->data[i], 0, numScans, overlap) ||
     389                    !pmReadoutReadChunkMask(readout, maskFits->data[i], 0, numScans, overlap) ||
     390                    !pmReadoutReadChunkWeight(readout, weightFits->data[i], 0, numScans, overlap)) {
    364391                    psError(PS_ERR_IO, false, "Unable to read chunk %d for file %d", numChunk, i);
    365392                    psFree(readouts);
     393                    psFree(kernels);
    366394                    psFree(outRO);
    367395                    psFree(view);
    … …  
    370398            }
    371399
    372 #if 0
     400#ifdef TESTING
     401            {
     402                pmReadout *ro = readouts->data[0];
     403                psTrace("ppStack", 1, "Stack: [%d:%d,%d:%d]\n", ro->col0, ro->col0 + ro->image->numCols,
     404                        ro->row0, ro->row0 + ro->image->numRows);
     405            }
     406#endif
     407
    373408            if (!ppStackReadout(config, outRO, readouts)) {
    374409                psError(PS_ERR_UNKNOWN, false, "Unable to stack images.\n");
    375410                psFree(readouts);
     411                psFree(kernels);
    376412                psFree(outRO);
    377413                psFree(view);
    378414                return false;
    379415            }
    380 #else
    381             printf("Stack...\n");
    382 #endif
    383416
    384417            for (int i = 0; i < num && more; i++) {
    385418                pmReadout *readout = readouts->data[i];
    386419                assert(readout);
    387 
    388                 // Files, containing the FITS handle
    389                 pmFPAfile *imageFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV", i);
    390                 pmFPAfile *maskFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.MASK", i);
    391                 pmFPAfile *weightFile = pmFPAfileSelectSingle(config->files, "PPSTACK.INCONV.WEIGHT", i);
    392 
    393                 // FITS handles
    394                 psFits *fitsImage = imageFile->fits;
    395                 psFits *fitsMask = maskFile->fits;
    396                 psFits *fitsWeight = weightFile->fits;
    397 
    398                 more &= pmReadoutMore(readout, fitsImage, 0, numScans, overlap);
    399                 more &= pmReadoutMoreMask(readout, fitsMask, 0, numScans, overlap);
    400                 more &= pmReadoutMoreWeight(readout, fitsWeight, 0, numScans, overlap);
    401             }
    402         }
     420                more &= pmReadoutMore(readout, imageFits->data[i], 0, numScans, overlap);
     421                more &= pmReadoutMoreMask(readout, maskFits->data[i], 0, numScans, overlap);
     422                more &= pmReadoutMoreWeight(readout, weightFits->data[i], 0, numScans, overlap);
     423            }
     424        }
     425
    403426        psFree(readouts);
    404 
    405         // Get rid of the input files
    406         fileActivation(config, outCombineFiles, false);
    407         filesIterateUp(config);
     427        psFree(kernels);
     428        for (int i = 0; i < num; i++) {
     429            psFitsClose(imageFits->data[i]);
     430            psFitsClose(maskFits->data[i]);
     431            psFitsClose(weightFits->data[i]);
     432        }
    408433
    409434        if (psMetadataLookupBool(&mdok, config->arguments, "PHOTOMETRY") &&
    … …  
    432457
    433458        // Write out the output files
    434         fileActivation(config, outCombineFiles, true);
     459        fileActivation(config, combineFiles, true);
    435460        filesIterateUp(config);
    436461
Note: See TracChangeset for help on using the changeset viewer.