IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 27625


Ignore:
Timestamp:
Apr 6, 2010, 1:36:51 PM (16 years ago)
Author:
eugene
Message:

working on stackphot

Location:
branches/eam_branches/stackphot.20100406/psphot/src
Files:
1 deleted
9 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/stackphot.20100406/psphot/src/Makefile.am

    r26894 r27625  
    1717# Force recompilation of psphotVersion.c, since it gets the version information
    1818psphotVersion.c: psphotVersionDefinitions.h
    19 psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
    20         -$(RM) psphotVersionDefinitions.h
    21         $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
    22 FORCE: ;
     19# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
     20#       -$(RM) psphotVersionDefinitions.h
     21#       $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
     22# FORCE: ;
    2323
    2424libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotForced psphotMakePSF
     27bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack
    2828# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy
    2929
    … …  
    3939psphotMakePSF_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    4040psphotMakePSF_LDADD = libpsphot.la
     41
     42psphotStack_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     43psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     44psphotStack_LDADD = libpsphot.la
    4145
    4246# psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    8084        psphotMosaicChip.c         \
    8185        psphotCleanup.c
     86
     87# a psphot-variant for stack photometry
     88psphotStack_SOURCES = \
     89        psphotStack.c              \
     90        psphotStackArguments.c     \
     91        psphotStackParseCamera.c   \
     92        psphotStackImageLoop.c     \
     93        psphotStackReadout.c       \
     94        psphotStackChisqImage.c    \
     95        psphotCleanup.c
     96
     97# psphotFitSourceLinearStack.c   
     98
    8299
    83100# # psphot analysis of the detectability of specified positions

  • branches/eam_branches/stackphot.20100406/psphot/src/psphot.h

    r27532 r27625  
    341341bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view);
    342342
     343/**** psphotStack prototypes ****/
     344
     345pmConfig *psphotStackArguments(int argc, char **argv);
     346bool psphotStackParseCamera (pmConfig *config);
     347bool psphotStackImageLoop (pmConfig *config);
     348bool psphotStackReadout (pmConfig *config, const pmFPAview *view);
     349bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view);
     350bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
     351                                     const pmFPAview *view,
     352                                     pmReadout *chiReadout,
     353                                     char *filename,
     354                                     int index);
     355
     356
    343357#endif

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotFitSourcesLinearStack.c

    r27547 r27625  
    11# include "psphotInternal.h"
    2 
    3 // for now, let's store the detections on the readout->analysis for each readout
    4 bool psphotFitSourcesLinearStack (pmConfig *config, const pmFPAview *view, bool final)
    5 {
    6     bool status = true;
    7 
    8     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    9     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    10 
    11     // loop over the available readouts
    12     for (int i = 0; i < num; i++) {
    13         if (!psphotFitSourcesLinearReadoutStack (config, view, "PSPHOT.INPUT", i, final)) {
    14             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
    15             return false;
    16         }
    17     }
    18     return true;
    19 }
    20 
    21 bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) {
     2// XXX need array of covar factors for each image
     3// XXX define the 'good' / 'bad' flags?
     4
     5bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) {
    226
    237    bool status;
    … …  
    259    float y;
    2610    float f;
    27     // float r;
    2811
    2912    // select the appropriate recipe information
    … …  
    3114    assert (recipe);
    3215
    33     // find the currently selected readout
    34     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    35     psAssert (file, "missing file?");
    36 
    37     pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
    38     psAssert (readout, "missing readout?");
    39 
    40     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    41     psAssert (detections, "missing detections?");
    42 
    43     psArray *sources = detections->allSources;
    44     psAssert (sources, "missing sources?");
    45 
    46     if (!sources->n) {
    47         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
    48         return true;
    49     }
    50 
    51     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
    52     psAssert (sources, "missing psf?");
    53 
    5416    psTimerStart ("psphot.linear");
    5517
    … …  
    6527    maskVal |= markVal;
    6628
    67     // source analysis is done in spatial order
    68     sources = psArraySort (sources, pmSourceSortByY);
     29    // analysis is done in spatial order (to speed up overlap search)
     30    // sort by first element in each source list
     31    objects = psArraySort (objects, pmPhotObjSortByY);
    6932
    7033    // storage array for fitSources
    71     psArray *fitSources = psArrayAllocEmpty (sources->n);
    72 
    73     // option to limit analysis to a specific region
    74     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    75     psRegion AnalysisRegion = psRegionFromString (region);
    76     AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
    77     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    78 
    79     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
    80         psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
    81     if (!status) {
    82         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
    83     }
    84     int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
    85     if (!status) {
    86         SKY_FIT_ORDER = 0;
    87     }
    88     bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
    89     if (!status) {
    90         SKY_FIT_LINEAR = false;
    91     }
    92 
    93     // XXX test: choose a larger-than expected radius:
    94     float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
    95     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
    96 
    97     // XXX do not apply covarFactor for the moment...
    98     // covarFactor = 1.0;
     34    psArray *fitSources = psArrayAllocEmpty (objects->n);
     35
     36    bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
     37    psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
     38
     39    // XXX store a local static array of covar factors for each of the images (by image ID)
     40    // float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
     41    // psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
    9942
    10043    // select the sources which will be used for the fitting analysis
    … …  
    10447        if (!object->sources) continue;
    10548
    106         // check an element of the group to see if we should use it
    107         if (!object->flags & PM_PHOT_OBJ_BAD) continue;
     49        // XXX check an element of the group to see if we should use it
     50        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
    10851
    10952        for (int j = 0; j < object->sources->n; j++) {
    … …  
    162105            pmSource *SRCj = fitSources->data[j];
    163106
    164             // XXX I need to know if this source is on the same image as SRCi --
    165             if (!sameImge) { continue; }
     107            // we only need to generate dot terms for source on the same image
     108            if (SRCj->imageID != SRCi->imageID) { continue; }
    166109
    167110            // skip over disjoint source images, break after last possible overlap
    … …  
    190133    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
    191134
    192     // XXXX **** philosophical question:
    193     // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
    194     // should retain the chisq and errors from the intermediate non-linear fit?
    195     // the non-linear fit provides better values for the position errors, and for
    196     // extended sources, the shape errors
    197 
    198135    // adjust I0 for fitSources and subtract
    199136    for (int i = 0; i < fitSources->n; i++) {
    … …  
    208145        model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
    209146        model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
    210         // XXX is the value of 'errors' modified by the sky fit?
    211147
    212148        // subtract object
    … …  
    229165    psFree (fitSources);
    230166    psFree (norm);
    231     psFree (skyfit);
    232167    psFree (errors);
    233     psFree (border);
    234168
    235169    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
    236 
    237     psphotVisualShowResidualImage (readout);
    238     psphotVisualPlotChisq (sources);
    239     // psphotVisualShowFlags (sources);
    240 
    241     // We have to place this visualization here because the models are not realized until
    242     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    243     psphotVisualShowPSFStars (recipe, psf, sources);
    244170
    245171    return true;
    246172}
    247173
    248 // XXX do we need this?
    249 // XXX disallow the simultaneous sky fit and remove this code...
    250 
    251 // Calculate the weight terms for the sky fit component of the matrix.  This function operates
    252 // on the pixels which correspond to all of the sources of interest.  These elements fill in
    253 // the border matrix components in the sparse matrix equation.
    254 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal) {
    255 
    256     // generate the image-wide weight terms
    257     // turn on MARK for all image pixels
    258     psRegion fullArray = psRegionSet (0, 0, 0, 0);
    259     fullArray = psRegionForImage (readout->mask, fullArray);
    260     psImageMaskRegion (readout->mask, fullArray, "OR", markVal);
    261 
    262     // turn off MARK for all object pixels
    263     for (int i = 0; i < sources->n; i++) {
    264         pmSource *source = sources->data[i];
    265         pmModel *model = pmSourceGetModel (NULL, source);
    266         if (model == NULL) continue;
    267         float x = model->params->data.F32[PM_PAR_XPOS];
    268         float y = model->params->data.F32[PM_PAR_YPOS];
    269         psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
    270     }
    271 
    272     // accumulate the image statistics from the masked regions
    273     psF32 **image  = readout->image->data.F32;
    274     psF32 **variance = readout->variance->data.F32;
    275     psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
    276 
    277     double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
    278     w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
    279 
    280     int col0 = readout->image->col0;
    281     int row0 = readout->image->row0;
    282 
    283     for (int j = 0; j < readout->image->numRows; j++) {
    284         for (int i = 0; i < readout->image->numCols; i++) {
    285             if (mask[j][i]) continue;
    286             if (constant_weights) {
    287                 wt = 1.0;
    288             } else {
    289                 wt = variance[j][i];
    290             }
    291             f = image[j][i];
    292             w   += 1/wt;
    293             fo  += f/wt;
    294             if (SKY_FIT_ORDER == 0) continue;
    295 
    296             xc  = i + col0;
    297             yc  = j + row0;
    298             x  +=    xc/wt;
    299             y  +=    yc/wt;
    300             x2 += xc*xc/wt;
    301             xy += xc*yc/wt;
    302             y2 += yc*yc/wt;
    303             fx +=  f*xc/wt;
    304             fy +=  f*yc/wt;
    305         }
    306     }
    307 
    308     // turn off MARK for all image pixels
    309     psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal));
    310 
    311     // set the Border T elements
    312     psSparseBorderElementG (border, 0, fo);
    313     psSparseBorderElementT (border, 0, 0, w);
    314     if (SKY_FIT_ORDER > 0) {
    315         psSparseBorderElementG (border, 0, fx);
    316         psSparseBorderElementG (border, 0, fy);
    317         psSparseBorderElementT (border, 1, 0, x);
    318         psSparseBorderElementT (border, 2, 0, y);
    319         psSparseBorderElementT (border, 0, 1, x);
    320         psSparseBorderElementT (border, 1, 1, x2);
    321         psSparseBorderElementT (border, 2, 1, xy);
    322         psSparseBorderElementT (border, 0, 2, y);
    323         psSparseBorderElementT (border, 1, 2, xy);
    324         psSparseBorderElementT (border, 2, 2, y2);
    325     }
    326 
    327     return true;
     174// sort by Y (ascending)
     175int pmPhotObjSortBySN (const void **a, const void **b)
     176{
     177    pmPhotObj *objA = *(pmPhotObj **)a;
     178    pmPhotObj *objB = *(pmPhotObj **)b;
     179
     180    psAssert (objA->sources, "missing sources?");
     181    psAssert (objB->sources, "missing sources?");
     182
     183    psAssert (objA->sources->n, "missing sources");
     184    psAssert (objB->sources->n, "missing sources");
     185
     186    psAssert (objA->sources->data[0], "missing sources");
     187    psAssert (objB->sources->data[0], "missing sources");
     188
     189    pmSource *A = objA->sources->data[0];
     190    pmSource *B = objB->sources->data[0];
     191
     192    psF32 fA = (A->peak == NULL) ? 0 : A->peak->y;
     193    psF32 fB = (B->peak == NULL) ? 0 : B->peak->y;
     194
     195    psF32 diff = fA - fB;
     196    if (diff > FLT_EPSILON) return (+1);
     197    if (diff < FLT_EPSILON) return (-1);
     198    return (0);
    328199}
     200

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotReadout.c

    r27532 r27625  
    3434
    3535    // Generate the mask and weight images, including the user-defined analysis region of interest
    36     psphotSetMaskAndVariance (config, view);
     36    if (!psphotSetMaskAndVariance (config, view)) {
     37        return psphotReadoutCleanup(config, view);
     38    }
    3739    if (!strcasecmp (breakPt, "NOTHING")) {
    3840        return psphotReadoutCleanup(config, view);

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotStackArguments.c

    r26894 r27625  
    22
    33static void usage(pmConfig *config, int exitCode);
    4 static void writeHelpInfo(pmConfig* config, FILE* ofile);
     4static void writeHelpInfo(FILE* ofile);
    55
    66pmConfig *psphotStackArguments(int argc, char **argv) {
    77
    88    int N;
    9     bool status;
    109
    1110    // print help info
    12     if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(argv[0], config, stdout);
    13     if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(argv[0], config, stdout);
     11    if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout);
     12    if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout);
    1413
    1514    // load config data from default locations
    … …  
    2928    // Number of threads is handled
    3029    PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
    31 
    32     // photcode : used in output to supplement header data (argument or recipe?)
    33     if ((N = psArgumentGet (argc, argv, "-photcode"))) {
    34         if (argc <= N+1) {
    35           psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    36           exit(PS_EXIT_CONFIG_ERROR);
    37         }
    38         psArgumentRemove (N, &argc, argv);
    39         psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]);
    40         psArgumentRemove (N, &argc, argv);
    41     }
    4230
    4331    // visual : interactive display mode
    … …  
    9078}
    9179
    92 static void writeHelpInfo(pmConfig* config, FILE* ofile)
     80static void writeHelpInfo(FILE* ofile)
    9381{
    9482  fprintf(ofile,
    95           "Usage: psphotStack -input (input.mdc) outroot\n"
     83          "Usage: psphotStack -input (INPUTS.mdc) (OUTROOT)\n"
    9684          "\n"
    97           "where:\n"
    98           "  FileNameList is a text file containing filenames, one per line\n"
    99           "  MaskFileNameList is a text file of mask filenames, one per line\n"
    100           "  VarFileNameList is a text file of variance filenames, one per line\n"
    101           "  OutFileBaseName is the 'root name' for output files\n"
     85          "where INPUTS.mdc contains various METADATAs, each with:\n"
     86          "\tIMAGE(STR):     Image filename\n"
     87          "\tMASK(STR):      Mask filename\n"
     88          "\tVARIANCE(STR):  Variance map filename\n"
     89          "OUTROOT is the 'root name' for output files\n"
    10290          "\n"
    10391          "additional options:\n"
    … …  
    136124          "  -trace-levels\n"
    137125          "     print current trace levels\n");
    138     psFree(config);
    139     pmConfigDone();
    140126    psLibFinalize();
    141127    exit(PS_EXIT_SUCCESS);

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotStackChisqImage.c

    r27547 r27625  
    1212    psTimerStart ("psphot.chisq.image");
    1313
     14    pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.OUTPUT", 0);
     15    psAssert (chisqFile, "missing chisq image FPA?");
     16
     17    pmCell *chisqCell = pmFPAviewThisCell(view, chisqFile->fpa);
     18
    1419    // create a holder for the image
    15     pmReadout *chiImage = pmReadoutAlloc();
     20    pmReadout *chiReadout = pmFPAviewThisReadout(view, chisqFile->fpa);
     21    if (!chiReadout) {
     22      chiReadout = pmReadoutAlloc(chisqCell);
     23    } else {
     24      psMemIncrRefCounter(chiReadout);
     25    }
    1626
    1727    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    … …  
    2030    // loop over the available readouts
    2131    for (int i = 0; i < num; i++) {
    22         if (!psphotStackChisqImageAddReadout(config, view, chiImage, "PSPHOT.INPUT", i)) {
     32        if (!psphotStackChisqImageAddReadout(config, view, chiReadout, "PSPHOT.INPUT", i)) {
    2333            psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i);
    2434            return false;
    … …  
    2636    }
    2737
     38    num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     39    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     40
     41    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num);
     42    num++;
     43    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num);
     44
    2845    // need to save the resulting image somewhere (PSPHOT.INPUT?)
     46    if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) {
     47        psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files");
     48        return false;
     49    }
    2950
    3051    psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image"));
    3152
     53    psFree (chiReadout);
    3254    return true;
    3355}
    3456
    3557bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration
    36                                      pmFPAview *view,
     58                                     const pmFPAview *view,
    3759                                     pmReadout *chiReadout,
    3860                                     char *filename,

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotStackImageLoop.c

    r27547 r27625  
    77}
    88
    9 bool psphotImageLoop (pmConfig *config) {
    10 
    11     bool status;
    12     pmChip *chip;
    13     pmCell *cell;
    14     pmReadout *readout;
     9bool psphotStackImageLoop (pmConfig *config) {
    1510
    1611    pmFPAview *view = pmFPAviewAlloc (0);

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotStackParseCamera.c

    r26894 r27625  
    1515    }
    1616
    17     psMetadata *item == NULL;
     17    psMetadataItem *item = NULL;
    1818    int nInputs = 0;
    19     while ((item = psMetadataGetIndex(inputs, nInputs)) != NULL) {
     19    while ((item = psMetadataGet(inputs, nInputs)) != NULL) {
    2020        if (item->type != PS_DATA_METADATA) {
    2121            psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type METADATA", item->name);
    … …  
    3434        pmFPAfile *imageFile = defineFile(config, NULL, "PSPHOT.INPUT", image, PM_FPA_FILE_IMAGE); // File for image
    3535        if (!imageFile) {
    36             psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", i, image);
     36            psError(PS_ERR_UNKNOWN, false, "Unable to define file from image %d (%s)", nInputs, image);
    3737            return false;
    3838        }
    … …  
    4141        if (mask && strlen(mask) > 0) {
    4242            if (!defineFile(config, imageFile, "PSPHOT.INPUT.MASK", mask, PM_FPA_FILE_MASK)) {
    43                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask);
     43                psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", nInputs, mask);
    4444                return false;
    4545            }
    4646        }
    4747
    48         psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map
     48        psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map
    4949        if (variance && strlen(variance) > 0) {
    5050            if (!defineFile(config, imageFile, "PSPHOT.INPUT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) {
    51                 psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance);
     51                psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", nInputs, variance);
    5252                return false;
    5353            }
    … …  
    5858    psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs);
    5959
     60    // generate an pmFPAimage for the chisqImage
     61    pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE");
     62    if (!chisqImage) {
     63        psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE");
     64        return false;
     65    }
     66    pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PPIMAGE.CHISQ.MASK");
     67    if (!chisqMask) {
     68        psError(PS_ERR_IO, false, _("Unable to generate output file from PPIMAGE.CHISQ.MASK"));
     69        return NULL;
     70    }
     71    if (chisqMask->type != PM_FPA_FILE_MASK) {
     72        psError(PS_ERR_IO, true, "PPIMAGE.CHISQ.MASK is not of type MASK");
     73        return NULL;
     74    }
     75    pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PPIMAGE.CHISQ.VARIANCE");
     76    if (!chisqVariance) {
     77        psError(PS_ERR_IO, false, _("Unable to generate output file from PPIMAGE.CHISQ.VARIANCE"));
     78        return NULL;
     79    }
     80    if (chisqVariance->type != PM_FPA_FILE_VARIANCE) {
     81        psError(PS_ERR_IO, true, "PPIMAGE.CHISQ.VARIANCE is not of type VARIANCE");
     82        return NULL;
     83    }
     84
     85# if (0)   
    6086    // define the additional input/output files associated with psphot
    6187    // XXX figure out which files are needed by psphotStack
    … …  
    6490        return false;
    6591    }
     92# endif
    6693
    6794    psTrace("psphot", 1, "Done with psphotStackParseCamera...\n");

  • branches/eam_branches/stackphot.20100406/psphot/src/psphotStackReadout.c

    r26894 r27625  
    1313        return false;
    1414    }
     15    // optional break-point for processing
     16    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
     17    PS_ASSERT_PTR_NON_NULL (breakPt, false);
    1518
    16     // set the photcode for this image (XXX currently goes into recipe, should it go into analysis?)
     19    // set the photcode for each image
    1720    if (!psphotAddPhotcode (config, view)) {
    1821        psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
    … …  
    2023    }
    2124
    22     // generate a background model (median, smoothed image)
    23     if (!psphotSetMaskAndVariance (config, view, recipe)) {
    24         return psphotReadoutCleanup (config, NULL, recipe, NULL, NULL, NULL);
     25    // Generate the mask and weight images
     26    if (!psphotSetMaskAndVariance (config, view)) {
     27        return psphotReadoutCleanup (config, view);
    2528    }
    2629    if (!strcasecmp (breakPt, "NOTHING")) {
    27         return psphotReadoutCleanup(config, NULL, recipe, NULL, NULL, NULL);
     30        return psphotReadoutCleanup (config, view);
    2831    }
    2932
    30     // optional break-point for processing
    31     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    32     PS_ASSERT_PTR_NON_NULL (breakPt, false);
    33 
    3433    // generate a background model (median, smoothed image)
     34    // XXX I think this is not defined correctly for an array of images.
    3535    if (!psphotModelBackground (config, view)) {
    36         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     36        return psphotReadoutCleanup (config, view);
    3737    }
    3838    if (!psphotSubtractBackground (config, view)) {
    39         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     39        return psphotReadoutCleanup (config, view);
    4040    }
    4141    if (!strcasecmp (breakPt, "BACKMDL")) {
    42         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     42        return psphotReadoutCleanup (config, view);
     43    }
     44
     45    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
     46    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
     47    if (!psphotLoadPSF (config, view)) {
     48        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
     49        return psphotReadoutCleanup (config, view);
     50    }
     51
     52    if (!psphotStackChisqImage(config, view)) {
     53        psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
     54        return psphotReadoutCleanup (config, view);
     55    }
     56    if (!strcasecmp (breakPt, "CHISQ")) {
     57        return psphotReadoutCleanup (config, view);
    4358    }
    4459
    4560    // find the detections (by peak and/or footprint) in the image.
    46     pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
    47     if (!detections) {
    48         psLogMsg ("psphot", 3, "unable to find detections in this image");
    49         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
     61    if (!psphotFindDetections (config, view, true)) { // pass 1
     62        // this only happens if we had an error in psphotFindDetections
     63        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
     64        return psphotReadoutCleanup (config, view);
    5065    }
    5166
    52     // construct sources and measure basic stats
    53     psArray *sources = psphotSourceStats (config, readout, detections, true);
    54     if (!sources) return false;
     67    // construct sources and measure basic stats (saved on detections->newSources)
     68    if (!psphotSourceStats (config, view, true)) { // pass 1
     69        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
     70        return psphotReadoutCleanup (config, view);
     71    }
    5572    if (!strcasecmp (breakPt, "PEAKS")) {
    56         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     73        return psphotReadoutCleanup(config, view);
    5774    }
    5875
    59     // find blended neighbors of very saturated stars
    60     // XXX merge this with Basic Deblend?
    61     psphotDeblendSatstars (readout, sources, recipe);
     76    // find blended neighbors of very saturated stars (detections->newSources)
     77    if (!psphotDeblendSatstars (config, view)) {
     78        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
     79        return psphotReadoutCleanup (config, view);
     80    }
    6281
    63     // mark blended peaks PS_SOURCE_BLEND
    64     if (!psphotBasicDeblend (sources, recipe)) {
    65         psLogMsg ("psphot", 3, "failed on deblend analysis");
    66         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     82    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
     83    if (!psphotBasicDeblend (config, view)) {
     84        psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
     85        return psphotReadoutCleanup (config, view);
    6786    }
    6887
    6988    // classify sources based on moments, brightness
    70     if (!psphotRoughClass (readout, sources, recipe, havePSF)) {
    71         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    72         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     89    if (!psphotRoughClass (config, view)) {
     90        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
     91        return psphotReadoutCleanup (config, view);
     92    }
     93    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     94    if (!psphotImageQuality (config, view)) { // pass 1
     95        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
     96        return psphotReadoutCleanup(config, view);
    7397    }
    7498    if (!strcasecmp (breakPt, "MOMENTS")) {
    75         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     99        return psphotReadoutCleanup (config, view);
    76100    }
    77101
    78     if (!havePSF && !psphotImageQuality (recipe, sources)) {
    79         psLogMsg("psphot", 3, "failed to measure image quality");
    80         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     102    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     103    // this step is skipped
     104    if (!psphotChoosePSF (config, view)) { // pass 1
     105        psLogMsg ("psphot", 3, "failure to construct a psf model");
     106        return psphotReadoutCleanup (config, view);
     107    }
     108    if (!strcasecmp (breakPt, "PSFMODEL")) {
     109        return psphotReadoutCleanup (config, view);
    81110    }
    82111
    83     // if we were not supplied a PSF, choose one here
    84     if (psf == NULL) {
    85         // use bright stellar objects to measure PSF
    86         // XXX if we do not have enough stars to generate the PSF, build one
    87         // from the SEEING guess and model class
    88         psf = psphotChoosePSF (readout, sources, recipe);
    89         if (psf == NULL) {
    90             psLogMsg ("psphot", 3, "failure to construct a psf model");
    91             return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    92         }
    93         havePSF = true;
    94     }
    95     if (!strcasecmp (breakPt, "PSFMODEL")) {
    96         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    97     }
    98     psphotVisualShowPSFModel (readout, psf);
    99 
    100112    // include externally-supplied sources
    101     psphotLoadExtSources (config, view, sources);
     113    // XXX fix this in the new multi-input context
     114    // psphotLoadExtSources (config, view); // pass 1
    102115
    103116    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    104     psphotGuessModels (config, readout, sources, psf);
     117    psphotGuessModels (config, view);
    105118
     119    // merge the newly selected sources into the existing list
     120    // NOTE: merge OLD and NEW
     121    psphotMergeSources (config, view);
     122
     123    // *** generate the objects (which unify the sources from the different images)
     124    // pmArray *objects = psphotMatchSources (config, view);
     125   
    106126    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    107     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    108 
    109     // We have to place these visualizations here because the models are not realized until
    110     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    111     psphotVisualShowPSFStars (recipe, psf, sources);
     127    // psphotFitSourcesLinearStack (config, objects, FALSE);
    112128
    113129    // identify CRs and extended sources
    114     psphotSourceSize (config, readout, sources, recipe, psf, 0);
    115     if (!strcasecmp (breakPt, "ENSEMBLE")) {
    116         goto finish;
    117     }
    118     psphotVisualShowSatStars (recipe, psf, sources);
    119 
    120     // non-linear PSF and EXT fit to brighter sources
    121     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    122     psphotBlendFit (config, readout, sources, psf);
    123 
    124     // replace all sources
    125     psphotReplaceAllSources (sources, recipe);
    126 
    127     // linear fit to include all sources (subtract again)
    128     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
    129 
    130     // if we only do one pass, skip to extended source analysis
    131     if (!strcasecmp (breakPt, "PASS1")) {
    132         goto pass1finish;
    133     }
    134     // NOTE: possibly re-measure background model here with objects subtracted
    135 
    136     // add noise for subtracted objects
    137     psphotAddNoise (readout, sources, recipe);
    138 
    139     // find fainter sources (pass 2)
    140     detections = psphotFindDetections (detections, readout, recipe);
    141 
    142     // remove noise for subtracted objects (ie, return to normal noise level)
    143     psphotSubNoise (readout, sources, recipe);
    144 
    145     // define new sources based on only the new peaks
    146     psArray *newSources = psphotSourceStats (config, readout, detections, false);
    147 
    148     // set source type
    149     if (!psphotRoughClass (readout, newSources, recipe, havePSF)) {
    150         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    151         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    152     }
    153 
    154     // create full input models, set the radius to fitRadius, set circular fit mask
    155     psphotGuessModels (config, readout, newSources, psf);
    156 
    157     // replace all sources so fit below applies to all at once
    158     psphotReplaceAllSources (sources, recipe);
    159 
    160     // merge the newly selected sources into the existing list
    161     psphotMergeSources (sources, newSources);
    162     psFree (newSources);
    163 
    164     // linear fit to all sources
    165     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
    166 
    167 pass1finish:
    168 
    169     // measure source size for the remaining sources
    170     psphotSourceSize (config, readout, sources, recipe, psf, 0);
    171 
    172     psphotExtendedSourceAnalysis (readout, sources, recipe);
    173 
    174     psphotExtendedSourceFits (readout, sources, recipe);
    175 
    176 finish:
    177 
    178     // plot positive sources
    179     // psphotSourcePlots (readout, sources, recipe);
     130    psphotSourceSize (config, view, TRUE);
    180131
    181132    // measure aperture photometry corrections
    182     if (!psphotApResid (config, readout, sources, psf)) {
     133    if (!psphotApResid (config, view)) {
    183134        psLogMsg ("psphot", 3, "failed on psphotApResid");
    184         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     135        return psphotReadoutCleanup (config, view);
    185136    }
    186137
    187138    // calculate source magnitudes
    188     psphotMagnitudes(config, readout, view, sources, psf);
     139    psphotMagnitudes(config, view);
    189140
    190     if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
     141    if (!psphotEfficiency(config, view)) {
    191142        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    192143        psErrorClear();
    … …  
    200151
    201152    // drop the references to the image pixels held by each source
    202     psphotSourceFreePixels (sources);
     153    psphotSourceFreePixels (config, view);
    203154
    204155    // create the exported-metadata and free local data
    205     return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     156    return psphotReadoutCleanup (config, view);
    206157}
    207158
Note: See TracChangeset for help on using the changeset viewer.