IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16404


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).

Location:
branches/pap_branch_080207/ppStack/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r16367 r16404  
    1212
    1313
     14#if 0
    1415// Define an output convolved image file
    1516static pmFPAfile *defineOutputConvolved(const char *name, // FPA file name
    … …  
    6465    return imageFile;
    6566}
    66 
     67#endif
    6768
    6869
    … …  
    196197        }
    197198
     199#if 0
    198200        // Output convolved files
    199201        pmFPAfile *outconvImage  = defineOutputConvolved("PPSTACK.OUTCONV", imageFile->fpa, config,
    … …  
    217219            return false;
    218220        }
    219 
     221#endif
    220222
    221223        psMetadataAddF32(imageFile->fpa->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0,
    … …  
    321323    }
    322324
     325
    323326    // Output PSF
    324327    if (havePSFs) {

  • 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

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

    r16382 r16404  
    1010#include "ppStack.h"
    1111
     12
     13#define ARRAY_BUFFER 16                 // Number to add to array at a time
    1214#define WCS_TOLERANCE 0.001             // Tolerance for WCS
    1315
    … …  
    2628    psMaskType maskBad = psMetadataLookupU8(NULL, config->arguments, "MASK.BAD"); // Value to mask
    2729    psMaskType maskBlank = psMetadataLookupU8(NULL, config->arguments, "MASK.BLANK"); // Mask for blank reg.
    28 //    float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
     30    float threshold = psMetadataLookupF32(NULL, config->arguments, "THRESHOLD.MASK"); // Threshold for mask deconvolution
    2931
    3032    int num = readouts->n;              // Number of inputs
    … …  
    4951            weighting = 1.0;
    5052        }
    51 
    52 #if 0
    53         // I don't think 'scale' is needed, since the convolution matches the scales
    54         float scale = psMetadataLookupF32(&mdok, fpa->analysis, "PPSTACK.SCALE"); // Rel. scale
    55         if (!mdok || !isfinite(scale)) {
    56             psWarning("No SCALE supplied for image %d --- set to unity.", i);
    57             scale = 1.0;
    58         }
    59 #endif
    6053
    6154        float exposure = psMetadataLookupF32(NULL, ro->parent->concepts, "CELL.EXPOSURE"); // Exposure time
    … …  
    135128#endif
    136129
    137 
    138 
    139     // XXX THIS NEEDS WORK
    140 #if 0
    141130    for (int i = 0; i < num; i++) {
    142131        pmStackData *data = stack->data[i]; // Data for this image
    143132        pmReadout *readout = data->readout; // Readout for this image
    144 
    145 
    146         // XXX Need some other mechanism to get the subtraction kernel
    147 
    148133
    149134        // Extract the regions and solutions used in the image matching
    … …  
    151136        {
    152137            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
    153                                                                "^SUBTRACTION.REGION$"); // Iterator
     138                                                               PM_SUBTRACTION_ANALYSIS_REGION); // Iterator
    154139            psMetadataItem *item = NULL;// Item from iteration
    155140            while ((item = psMetadataGetAndIncrement(iter))) {
    … …  
    159144            psFree(iter);
    160145        }
    161         psArray *solutions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of solutions
     146        psArray *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels
    162147        {
    163148            psMetadataIterator *iter = psMetadataIteratorAlloc(readout->analysis, PS_LIST_HEAD,
    164                                                                "^SUBTRACTION.SOLUTION$"); // Iterator
     149                                                               PM_SUBTRACTION_ANALYSIS_KERNEL); // Iterator
    165150            psMetadataItem *item = NULL;// Item from iteration
    166151            while ((item = psMetadataGetAndIncrement(iter))) {
    167152                assert(item->type == PS_DATA_VECTOR);
    168                 solutions = psArrayAdd(solutions, ARRAY_BUFFER, item->data.V);
     153                kernels = psArrayAdd(kernels, ARRAY_BUFFER, item->data.V);
    169154            }
    170155            psFree(iter);
    171156        }
    172         assert(regions->n == solutions->n);
    173 
    174         pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis,
    175                                                             "SUBTRACTION.KERNEL"); // Kernels
    176 
    177         psPixels *reject = pmStackReject(data->pixels, threshold, regions,
    178                                          solutions, kernels); // List of pixels to reject
     157        assert(regions->n == kernels->n);
     158
     159        psPixels *reject = pmStackReject(data->pixels, threshold, regions, kernels); // Pixels to reject
    179160        psFree(data->pixels);
    180161        data->pixels = reject;
    181162
    182         psFree(solutions);
     163        psFree(kernels);
    183164        psFree(regions);
    184165    }
    185 #endif
    186 
    187166
    188167#ifdef REJECTION_FILES
Note: See TracChangeset for help on using the changeset viewer.