IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 17, 2009, 12:08:50 PM (17 years ago)
Author:
beaumont
Message:

merged with head

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301

  • branches/cnb_branches/cnb_branch_20090301/psphot

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/Makefile.am

    r21392 r23352  
    11lib_LTLIBRARIES = libpsphot.la
    2 libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     2
     3# PSPHOT_VERSION=`if [ -e ../../VERSION ]; then cat ../../VERSION; else svnversion; fi`
     4# PSPHOT_BRANCH=`if [ -e ../../BRANCH ]; then cat ../../BRANCH; else svn info | sed -n -e '/URL:/ h' -e '/Repository Root:/ { x; H; x; s|Repository Root: \(.*\)\nURL: \1\(.*\)|\2| ; s|^/|| ; s|/[a-zA-Z]*/src.*|| ; p }'; fi`
     5# PSPHOT_SOURCE=`if [ -e ../../SOURCE ]; then cat ../../SOURCE; else svn info | sed -n -e 's/Repository UUID: // p'; fi`
     6#
     7# # Force recompilation of psphotVersion.c, since it gets the version information
     8# psphotVersion.c: FORCE
     9#       touch psphotVersion.c
     10# FORCE: ;
     11
     12libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -DPSPHOT_VERSION=$(PSPHOT_VERSION) -DPSPHOT_BRANCH=$(PSPHOT_BRANCH) -DPSPHOT_SOURCE=$(SVN_SOURCE)
    313libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    414

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphot.c

    r20411 r23352  
    99
    1010    psTimerStart ("complete");
    11     pmErrorRegister();                  // register psModule's error codes/messages
     11    pmErrorRegister();                  // register psModule's error codes/messages
    1212    psphotInit();
    1313
    … …  
    1515    pmConfig *config = psphotArguments (argc, argv);
    1616    if (!config) {
    17         psErrorStackPrint(stderr, "Error reading arguments\n");
    18         usage ();
     17        psErrorStackPrint(stderr, "Error reading arguments\n");
     18        usage ();
    1919    }
     20
     21    psphotVersionPrint();
    2022
    2123    // load input data (config and images (signal, noise, mask)

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphot.h

    r21392 r23352  
    1616const char     *psphotCVSName(void);
    1717psString        psphotVersion(void);
     18psString        psphotSource(void);
    1819psString        psphotVersionLong(void);
     20bool            psphotVersionHeader(psMetadata *header);
     21bool            psphotVersionHeaderFull(psMetadata *header);
     22void            psphotVersionPrint(void);
    1923
    2024bool            psphotModelTest (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
    … …  
    167171
    168172// psphotVisual functions
    169 bool psphotSetVisual (bool mode);
    170173bool psphotVisualShowImage (pmReadout *readout);
    171174bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout);

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotArguments.c

    r21458 r23352  
    8282    if ((N = psArgumentGet (argc, argv, "-visual"))) {
    8383        psArgumentRemove (N, &argc, argv);
    84         psphotSetVisual (true);
    85         // pmSourceSetVisual (true);
     84        pmVisualSetVisual(true);
    8685    }
    8786

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotImageLoop.c

    r21458 r23352  
    2626
    2727    pmFPAview *view = pmFPAviewAlloc (0);
     28    pmHDU *lastHDU = NULL;              // Last HDU updated
    2829
    2930    // files associated with the science image
    … …  
    6465                psLogMsg ("psphot", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    6566                if (! readout->data_exists) { continue; }
     67
     68                // Update the header
     69                {
     70                    pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
     71                    if (hdu && hdu != lastHDU) {
     72                        psphotVersionHeaderFull(hdu->header);
     73                        lastHDU = hdu;
     74                    }
     75                }
    6676
    6777                // run the actual photometry analysis on this chip/cell/readout

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotMagnitudes.c

    r21519 r23352  
    1515    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    1616    if (!status) {
    17         nThreads = 0;
     17        nThreads = 0;
    1818    }
    1919
    … …  
    3636
    3737    // the binning details are saved on the analysis metadata
    38     psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING");
     38    psImageBinning *binning = NULL;
     39    if (backModel) {
     40        binning = psMetadataLookupPtr(&status, backModel->analysis, "PSPHOT.BACKGROUND.BINNING");
     41    }
    3942
    4043    bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
    … …  
    5356    for (int i = 0; i < cellGroups->n; i++) {
    5457
    55         psArray *cells = cellGroups->data[i];
    56 
    57         for (int j = 0; j < cells->n; j++) {
    58 
    59             // allocate a job -- if threads are not defined, this just runs the job
    60             psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES");
    61 
    62             psArrayAdd(job->args, 1, cells->data[j]); // sources
    63             psArrayAdd(job->args, 1, psf);
    64             psArrayAdd(job->args, 1, binning);
    65             psArrayAdd(job->args, 1, backModel);
    66             psArrayAdd(job->args, 1, backStdev);
    67 
    68             PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
    69             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    70             PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
    71 
    72             if (!psThreadJobAddPending(job)) {
    73                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    74                 psFree (job);
    75                 return false;
    76             }
    77             psFree(job);
     58        psArray *cells = cellGroups->data[i];
     59
     60        for (int j = 0; j < cells->n; j++) {
     61
     62            // allocate a job -- if threads are not defined, this just runs the job
     63            psThreadJob *job = psThreadJobAlloc ("PSPHOT_MAGNITUDES");
     64
     65            psArrayAdd(job->args, 1, cells->data[j]); // sources
     66            psArrayAdd(job->args, 1, psf);
     67            psArrayAdd(job->args, 1, binning);
     68            psArrayAdd(job->args, 1, backModel);
     69            psArrayAdd(job->args, 1, backStdev);
     70
     71            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
     72            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     73            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for nAp
     74
     75            if (!psThreadJobAddPending(job)) {
     76                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     77                psFree (job);
     78                return false;
     79            }
     80            psFree(job);
    7881
    7982# if (0)
    80                 int nap = 0;
    81                 if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) {
    82                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    83                     return false;
    84                 }
    85                 Nap += nap;
     83                int nap = 0;
     84                if (!psphotMagnitudes_Unthreaded (&nap, cells->data[j], psf, binning, backModel, backStdev, photMode, maskVal)) {
     85                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     86                    return false;
     87                }
     88                Nap += nap;
    8689# endif
    87         }
    88 
    89         // wait for the threads to finish and manage results
    90         if (!psThreadPoolWait (false)) {
    91             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    92             return false;
    93         }
    94 
    95         // we have only supplied one type of job, so we can assume the types here
    96         psThreadJob *job = NULL;
    97         while ((job = psThreadJobGetDone()) != NULL) {
    98             if (job->args->n < 1) {
    99                 fprintf (stderr, "error with job\n");
    100             } else {
    101                 psScalar *scalar = job->args->data[7];
    102                 Nap += scalar->data.S32;
    103             }
    104             psFree(job);
    105         }
     90        }
     91
     92        // wait for the threads to finish and manage results
     93        if (!psThreadPoolWait (false)) {
     94            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     95            return false;
     96        }
     97
     98        // we have only supplied one type of job, so we can assume the types here
     99        psThreadJob *job = NULL;
     100        while ((job = psThreadJobGetDone()) != NULL) {
     101            if (job->args->n < 1) {
     102                fprintf (stderr, "error with job\n");
     103            } else {
     104                psScalar *scalar = job->args->data[7];
     105                Nap += scalar->data.S32;
     106            }
     107            psFree(job);
     108        }
    106109    }
    107110
    … …  
    130133        if (status && isfinite(source->apMag)) Nap ++;
    131134
    132         if (backModel) {
    133             psAssert (binning, "if backModel is defined, so should binning be");
    134             source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
    135             if (isnan(source->sky) && false) {
    136                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
    137             }
    138         } else {
    139             source->sky = NAN;
    140         }
    141 
    142         if (backStdev) {
    143             psAssert (binning, "if backStdev is defined, so should binning be");
    144             source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
    145             if (isnan(source->skyErr) && false) {
    146                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
    147             }
    148         } else {
    149             source->skyErr = NAN;
    150         }
     135        if (backModel) {
     136            psAssert (binning, "if backModel is defined, so should binning be");
     137            source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
     138            if (isnan(source->sky) && false) {
     139                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
     140            }
     141        } else {
     142            source->sky = NAN;
     143        }
     144
     145        if (backStdev) {
     146            psAssert (binning, "if backStdev is defined, so should binning be");
     147            source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
     148            if (isnan(source->skyErr) && false) {
     149                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
     150            }
     151        } else {
     152            source->skyErr = NAN;
     153        }
    151154    }
    152155
    … …  
    169172        if (status && isfinite(source->apMag)) Nap ++;
    170173
    171         if (backModel) {
    172             psAssert (binning, "if backModel is defined, so should binning be");
    173             source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
    174             if (isnan(source->sky) && false) {
    175                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
    176             }
    177         } else {
    178             source->sky = NAN;
    179         }
    180 
    181         if (backStdev) {
    182             psAssert (binning, "if backStdev is defined, so should binning be");
    183             source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
    184             if (isnan(source->skyErr) && false) {
    185                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
    186             }
    187         } else {
    188             source->skyErr = NAN;
    189         }
     174        if (backModel) {
     175            psAssert (binning, "if backModel is defined, so should binning be");
     176            source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
     177            if (isnan(source->sky) && false) {
     178                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
     179            }
     180        } else {
     181            source->sky = NAN;
     182        }
     183
     184        if (backStdev) {
     185            psAssert (binning, "if backStdev is defined, so should binning be");
     186            source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
     187            if (isnan(source->skyErr) && false) {
     188                psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
     189            }
     190        } else {
     191            source->skyErr = NAN;
     192        }
    190193    }
    191194
    … …  
    210213    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    211214    if (!status) {
    212         nThreads = 0;
     215        nThreads = 0;
    213216    }
    214217    nThreads = 0; // XXX until testing is complete, do not thread this function
    … …  
    233236    for (int i = 0; i < cellGroups->n; i++) {
    234237
    235         psArray *cells = cellGroups->data[i];
    236 
    237         for (int j = 0; j < cells->n; j++) {
    238 
    239             // allocate a job -- if threads are not defined, this just runs the job
    240             psThreadJob *job = psThreadJobAlloc ("PSPHOT_PSF_WEIGHTS");
    241 
    242             psArrayAdd(job->args, 1, cells->data[j]); // sources
    243             PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
    244 
    245             if (!psThreadJobAddPending(job)) {
    246                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    247                 psFree (job);
    248                 return false;
    249             }
    250             psFree(job);
    251 
    252         }
    253 
    254         // wait for the threads to finish and manage results
    255         if (!psThreadPoolWait (false)) {
    256             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    257             return false;
    258         }
    259 
    260         // we have only supplied one type of job, so we can assume the types here
    261         psThreadJob *job = NULL;
    262         while ((job = psThreadJobGetDone()) != NULL) {
    263             if (job->args->n < 1) {
    264                 fprintf (stderr, "error with job\n");
    265             }
    266             psFree(job);
    267         }
     238        psArray *cells = cellGroups->data[i];
     239
     240        for (int j = 0; j < cells->n; j++) {
     241
     242            // allocate a job -- if threads are not defined, this just runs the job
     243            psThreadJob *job = psThreadJobAlloc ("PSPHOT_PSF_WEIGHTS");
     244
     245            psArrayAdd(job->args, 1, cells->data[j]); // sources
     246            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
     247
     248            if (!psThreadJobAddPending(job)) {
     249                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     250                psFree (job);
     251                return false;
     252            }
     253            psFree(job);
     254
     255        }
     256
     257        // wait for the threads to finish and manage results
     258        if (!psThreadPoolWait (false)) {
     259            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     260            return false;
     261        }
     262
     263        // we have only supplied one type of job, so we can assume the types here
     264        psThreadJob *job = NULL;
     265        while ((job = psThreadJobGetDone()) != NULL) {
     266            if (job->args->n < 1) {
     267                fprintf (stderr, "error with job\n");
     268            }
     269            psFree(job);
     270        }
    268271    }
    269272
    … …  
    285288        pmSource *source = (pmSource *) sources->data[i];
    286289
    287         // we must have a valid model
    288         pmModel *model = pmSourceGetModel (&isPSF, source);
    289         if (model == NULL) {
    290           psTrace ("psphot", 3, "fail mag : no valid model");
    291           source->pixWeight = NAN;
    292           continue;
    293         }
     290        // we must have a valid model
     291        pmModel *model = pmSourceGetModel (&isPSF, source);
     292        if (model == NULL) {
     293          psTrace ("psphot", 3, "fail mag : no valid model");
     294          source->pixWeight = NAN;
     295          continue;
     296        }
    294297
    295298        status = pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
    296         if (!status) {
    297           psTrace ("psphot", 3, "fail to measure pixel weight");
    298           source->pixWeight = NAN;
    299           continue;
    300         }
    301 
    302     }
    303 
    304     return true;
    305 }
     299        if (!status) {
     300          psTrace ("psphot", 3, "fail to measure pixel weight");
     301          source->pixWeight = NAN;
     302          continue;
     303        }
     304
     305    }
     306
     307    return true;
     308}

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotMaskReadout.c

    r21366 r23352  
    2929    }
    3030
     31    // make this an option via the recipe
     32    if (0) {
     33      psImage *im = readout->image;
     34      psImage *wt = readout->variance;
     35      psImage *mk = readout->mask;
     36      for (int j = 0; j < im->numRows; j++) {
     37        for (int i = 0; i < im->numCols; i++) {
     38          if (isfinite(im->data.F32[j][i]) && isfinite(wt->data.F32[j][i])) continue;
     39          mk->data.PS_TYPE_IMAGE_MASK_DATA[j][i] |= maskBad;
     40        }
     41      }
     42    }
     43
    3144    // mask the excluded outer pixels
    3245    // these coordinates refer to the parent image

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotModelBackground.c

    r21183 r23352  
    3131static bool backgroundModel(psImage *model,  // Model image
    3232                            psImage *modelStdev, // Model stdev image
    33                             psImage *image, // Image for which to generate a background model
    34                             psImage *mask, // Mask for image
     33                            psMetadata *analysis, // Analysis metadata for outputs
     34                            pmReadout *readout, // Readout for which to generate a background model
    3535                            psImageBinning *binning, // Binning parameters
    3636                            const pmConfig *config // Configuration
    … …  
    4141    bool status = true;
    4242
     43    psImage *image = readout->image, *mask = readout->mask; // Image and mask for readout
     44
    4345    // select the appropriate recipe information
    4446    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    … …  
    4951    assert (maskVal);
    5052
    51     // user supplied seed, if available
    52     unsigned long seed = psMetadataLookupS32 (&status, recipe, "IMSTATS_SEED");
    53     if (!status) {
    54         seed = 0;
    55     }
    56     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, seed);
     53    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    5754
    5855    // subtract this amount extra from the sky
    … …  
    140137
    141138    // we save the binning structure for use in psphotMagnitudes
    142     status = psMetadataAddPtr(recipe, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
    143     PS_ASSERT (status, false);
    144 
     139    psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING",
     140                     PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
    145141
    146142    psF32 **modelData = model->data.F32;
    … …  
    343339    psImage *modelStdev = psImageAlloc(binning->nXruff, binning->nYruff, PS_TYPE_F32); // Standard deviation
    344340
    345     if (!backgroundModel(model, modelStdev, ro->image, ro->mask, binning, config)) {
     341    if (!backgroundModel(model, modelStdev, ro->analysis, ro, binning, config)) {
    346342        psFree(model);
    347343        psFree(modelStdev);
    … …  
    365361    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, filename);
    366362    pmFPA *inFPA = file->fpa;
    367     pmReadout *readout = pmFPAviewThisReadout (view, inFPA);
    368     psImage *image = readout->image;
    369     psImage *mask  = readout->mask;
    370 
    371     psImageBinning *binning = backgroundBinning(image, config); // Image binning parameters
     363    pmReadout *readout = pmFPAviewThisReadout(view, inFPA);
     364
     365    psImageBinning *binning = backgroundBinning(readout->image, config); // Image binning parameters
    372366    pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning);
    373367    pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning);
    374368
    375     if (!backgroundModel(model->image, modelStdev->image, image, mask, binning, config)) {
     369    if (!backgroundModel(model->image, modelStdev->image, model->analysis, readout, binning, config)) {
    376370        psError(PS_ERR_UNKNOWN, false, "Unable to generate background model");
    377371        return false;

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotSetThreads.c

    r21392 r23352  
    2525    psFree(task);
    2626
    27     task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 4);
     27    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 5);
    2828    task->function = &psphotSourceStats_Threaded;
    2929    psThreadTaskAdd(task);

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotSourceStats.c

    r21519 r23352  
    1515    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    1616    if (!status) {
    17         nThreads = 0;
     17        nThreads = 0;
    1818    }
    1919
    … …  
    4141        pmSource *source = pmSourceAlloc();
    4242
    43         // add the peak
     43        // add the peak
    4444        source->peak = psMemIncrRefCounter(peak);
    4545
    46         // allocate space for moments
     46        // allocate space for moments
    4747        source->moments = pmMomentsAlloc();
    4848
    … …  
    5050        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    5151
    52         peak->assigned = true;
    53         psArrayAdd (sources, 100, source);
    54         psFree (source);
     52        peak->assigned = true;
     53        psArrayAdd (sources, 100, source);
     54        psFree (source);
    5555    }
    5656
    5757    if (!strcasecmp (breakPt, "PEAKS")) {
    58         psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
    59         psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
    60         psphotVisualShowMoments (sources);
    61         return sources;
     58        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
     59        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
     60        psphotVisualShowMoments (sources);
     61        return sources;
    6262    }
    6363
    … …  
    6565    int Nfail = 0;
    6666    int Nmoments = 0;
     67    int Nfaint = 0;
    6768
    6869    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    … …  
    7475    for (int i = 0; i < cellGroups->n; i++) {
    7576
    76         psArray *cells = cellGroups->data[i];
    77 
    78         for (int j = 0; j < cells->n; j++) {
    79 
    80             // allocate a job -- if threads are not defined, this just runs the job
    81             psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
    82 
    83             psArrayAdd(job->args, 1, cells->data[j]); // sources
    84             psArrayAdd(job->args, 1, recipe);
    85             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
    86             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    87 
    88             if (!psThreadJobAddPending(job)) {
    89                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    90                 psFree (job);
    91                 return NULL;
    92             }
    93             psFree(job);
     77        psArray *cells = cellGroups->data[i];
     78
     79        for (int j = 0; j < cells->n; j++) {
     80
     81            // allocate a job -- if threads are not defined, this just runs the job
     82            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
     83
     84            psArrayAdd(job->args, 1, cells->data[j]); // sources
     85            psArrayAdd(job->args, 1, recipe);
     86            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
     87            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     88            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
     89
     90            if (!psThreadJobAddPending(job)) {
     91                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     92                psFree (job);
     93                return NULL;
     94            }
     95            psFree(job);
    9496
    9597# if (0)
    96                 int nfail = 0;
    97                 int nmoments = 0;
    98                 if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
    99                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    100                     return NULL;
    101                 }
    102                 Nfail += nfail;
    103                 Nmoments += nmoments;
     98                int nfail = 0;
     99                int nmoments = 0;
     100                if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
     101                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     102                    return NULL;
     103                }
     104                Nfail += nfail;
     105                Nmoments += nmoments;
    104106# endif
    105         }
    106 
    107         // wait for the threads to finish and manage results
    108         if (!psThreadPoolWait (false)) {
    109             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    110             return NULL;
    111         }
    112 
    113         // we have only supplied one type of job, so we can assume the types here
    114         psThreadJob *job = NULL;
    115         while ((job = psThreadJobGetDone()) != NULL) {
    116             if (job->args->n < 1) {
    117                 fprintf (stderr, "error with job\n");
    118             } else {
    119                 psScalar *scalar = NULL;
    120                 scalar = job->args->data[2];
    121                 Nmoments += scalar->data.S32;
    122                 scalar = job->args->data[3];
    123                 Nfail += scalar->data.S32;
    124             }
    125             psFree(job);
    126         }
     107        }
     108
     109        // wait for the threads to finish and manage results
     110        if (!psThreadPoolWait (false)) {
     111            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     112            return NULL;
     113        }
     114
     115        // we have only supplied one type of job, so we can assume the types here
     116        psThreadJob *job = NULL;
     117        while ((job = psThreadJobGetDone()) != NULL) {
     118            if (job->args->n < 1) {
     119                fprintf (stderr, "error with job\n");
     120            } else {
     121                psScalar *scalar = NULL;
     122                scalar = job->args->data[2];
     123                Nmoments += scalar->data.S32;
     124                scalar = job->args->data[3];
     125                Nfail += scalar->data.S32;
     126                scalar = job->args->data[4];
     127                Nfaint += scalar->data.S32;
     128            }
     129            psFree(job);
     130        }
    127131    }
    128132
    129133    psFree (cellGroups);
    130134
    131     psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot.stats"));
     135    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
    132136
    133137    psphotVisualShowMoments (sources);
    … …  
    166170    int Nfail = 0;
    167171    int Nmoments = 0;
     172    int Nfaint = 0;
    168173    for (int i = 0; i < sources->n; i++) {
    169174        pmSource *source = sources->data[i];
    … …  
    171176        // skip faint sources for moments measurement
    172177        if (source->peak->SN < MIN_SN) {
    173             source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN;
     178            source->mode |= PM_SOURCE_MODE_BELOW_MOMENTS_SN;
     179            Nfaint++;
    174180            continue;
    175181        }
    … …  
    179185        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    180186        if (!status) {
    181             source->mode |= PM_SOURCE_MODE_SKY_FAILURE;
    182             psErrorClear(); // XXX re-consider the errors raised here
    183             Nfail ++;
    184             continue;
     187            source->mode |= PM_SOURCE_MODE_SKY_FAILURE;
     188            psErrorClear(); // XXX re-consider the errors raised here
     189            Nfail ++;
     190            continue;
    185191        }
    186192
    … …  
    189195        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    190196        if (!status) {
    191             source->mode |= PM_SOURCE_MODE_SKYVAR_FAILURE;
    192             Nfail ++;
    193             psErrorClear();
    194             continue;
     197            source->mode |= PM_SOURCE_MODE_SKYVAR_FAILURE;
     198            Nfail ++;
     199            psErrorClear();
     200            continue;
    195201        }
    196202
    … …  
    208214        status = pmSourceMoments (source, BIG_RADIUS);
    209215        if (status) {
    210             source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
     216            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    211217            Nmoments ++;
    212218            continue;
    213219        }
    214220
    215         source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
     221        source->mode |= PM_SOURCE_MODE_MOMENTS_FAILURE;
    216222        Nfail ++;
    217223        psErrorClear();
    … …  
    225231    scalar = job->args->data[3];
    226232    scalar->data.S32 = Nfail;
    227    
     233
     234    scalar = job->args->data[4];
     235    scalar->data.S32 = Nfaint;
     236
    228237    return true;
    229238}
    … …  
    268277        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    269278        if (!status) {
    270             psErrorClear(); // XXX re-consider the errors raised here
    271             Nfail ++;
    272             continue;
     279            psErrorClear(); // XXX re-consider the errors raised here
     280            Nfail ++;
     281            continue;
    273282        }
    274283
    … …  
    277286        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
    278287        if (!status) {
    279             Nfail ++;
    280             psErrorClear();
    281             continue;
     288            Nfail ++;
     289            psErrorClear();
     290            continue;
    282291        }
    283292
    … …  
    307316    *nmoments = Nmoments;
    308317    *nfail = Nfail;
    309    
     318
    310319    return true;
    311320}
    312 # endif 
     321# endif

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotSubtractBackground.c

    r21183 r23352  
    44// generate the median in NxN boxes, clipping heavily
    55// linear interpolation to generate full-scale model
    6 bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename) 
     6bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view, const char *filename)
    77{
    88    bool status = true;
    … …  
    3131    assert (maskVal);
    3232
    33     psImageBinning *binning = psMetadataLookupPtr(&status, recipe, "PSPHOT.BACKGROUND.BINNING");
     33    psImageBinning *binning = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.BACKGROUND.BINNING");
    3434    assert (binning);
    3535

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotVersion.c

    r12805 r23352  
    1 # include "psphotInternal.h"
     1#include "psphotInternal.h"
    22
    3 # if (HAVE_KAPA)
    4 # include <kapa.h>
    5 # endif
     3#ifdef HAVE_KAPA
     4#include <kapa.h>
     5#endif
    66
    7 static const char *cvsTag = "$Name: not supported by cvs2svn $";// CVS tag name
     7#ifndef PSPHOT_VERSION
     8#error "PSPHOT_VERSION is not set"
     9#endif
     10#ifndef PSPHOT_BRANCH
     11#error "PSPHOT_BRANCH is not set"
     12#endif
     13#ifndef PSPHOT_SOURCE
     14#error "PSPHOT_SOURCE is not set"
     15#endif
     16
     17#define xstr(s) str(s)
     18#define str(s) #s
    819
    920psString psphotVersion(void)
    1021{
    11     psString version = NULL;            // Version, to return
    12     psStringAppend(&version, "%s-%s",PACKAGE_NAME,PACKAGE_VERSION);
    13     return version;
     22    char *value = NULL;
     23    psStringAppend(&value, "%s@%s", xstr(PSPHOT_BRANCH), xstr(PSPHOT_VERSION));
     24    return value;
     25}
     26
     27psString psphotSource(void)
     28{
     29    return psStringCopy(xstr(PSPHOT_SOURCE));
    1430}
    1531
    1632psString psphotVersionLong(void)
    1733{
    18     psString version = psphotVersion(); // Version, to return
    19     psString tag = psStringStripCVS(cvsTag, "Name"); // CVS tag
    20     psStringAppend(&version, " (cvs tag %s) %s, %s", tag, __DATE__, __TIME__);
     34    psString version = psLibVersion();  // Version, to return
     35    psString source = psLibSource();    // Source
    2136
    22 # if (HAVE_KAPA)
    23     psString ohanaVersion = psStringStripCVS (ohana_version(), "Name");
    24     psString libdvoVersion = psStringStripCVS (libdvo_version(), "Name");
     37    psStringPrepend(&version, "psphot ");
     38    psStringAppend(&version, " from %s, built %s, %s", source, __DATE__, __TIME__);
     39    psFree(source);
    2540
    26     psStringAppend (&version, " with libkapa (ohana %s, libdvo: %s)", ohanaVersion, libdvoVersion);
    27     psFree (ohanaVersion);
    28     psFree (libdvoVersion);
    29 # else
    30     psStringAppend (&version, " WITHOUT libkapa");
    31 # endif
     41#ifdef __OPTIMIZE__
     42    psStringAppend(&version, " optimised");
     43#else
     44    psStringAppend(&version, " unoptimised");
     45#endif
    3246
    33     psFree(tag);
     47#ifdef HAVE_KAPA
     48#if 0
     49    // XXX Need to get ohana and libdvo versions
     50    psString ohanaVersion = psStringStripCVS(ohana_version(), "Name");
     51    psString libdvoVersion = psStringStripCVS(libdvo_version(), "Name");
     52    psStringAppend(&version, " with libkapa (ohana %s, libdvo: %s)", ohanaVersion, libdvoVersion);
     53    psFree(ohanaVersion);
     54    psFree(libdvoVersion);
     55#else
     56    psStringAppend(&version, " with libkapa");
     57#endif
     58
     59#else
     60    psStringAppend (&version, " without libkapa");
     61#endif
     62
    3463    return version;
    3564}
    3665
    37 // Defined by RHL; leaving for backwards compatibility.
    38 const char *psphotCVSName(void) {
    39    return cvsTag;
     66bool psphotVersionHeader(psMetadata *header)
     67{
     68    PS_ASSERT_METADATA_NON_NULL(header, false);
     69
     70    psString version = psphotVersion(); // Software version
     71    psString source = psphotSource();   // Software source
     72
     73    psStringPrepend(&version, "psphot version: ");
     74    psStringPrepend(&source, "psphot source: ");
     75
     76    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, version);
     77    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, source);
     78
     79    psFree(version);
     80    psFree(source);
     81
     82    return true;
    4083}
     84
     85
     86bool psphotVersionHeaderFull(psMetadata *header)
     87{
     88    PS_ASSERT_METADATA_NON_NULL(header, false);
     89
     90    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
     91    psString timeString = psTimeToISO(time); // The time in an ISO string
     92    psFree(time);
     93    psString history = NULL;               // History string
     94    psStringAppend(&history, "psphot at %s", timeString);
     95    psFree(timeString);
     96    psMetadataAddStr(header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, NULL, history);
     97    psFree(history);
     98
     99    psLibVersionHeader(header);
     100    psModulesVersionHeader(header);
     101    psphotVersionHeader(header);
     102
     103    return true;
     104}
     105
     106
     107void psphotVersionPrint(void)
     108{
     109    psTime *time = psTimeGetNow(PS_TIME_TAI); // The time now
     110    psString timeString = psTimeToISO(time); // The time in an ISO string
     111    psFree(time);
     112    psLogMsg("psphot", PS_LOG_INFO, "psphot at %s", timeString);
     113    psFree(timeString);
     114
     115    psString pslib = psLibVersionLong();// psLib version
     116    psString psmodules = psModulesVersionLong(); // psModules version
     117    psString psphot = psphotVersionLong(); // psphot version
     118
     119    psLogMsg("psphot", PS_LOG_INFO, "%s", pslib);
     120    psLogMsg("psphot", PS_LOG_INFO, "%s", psmodules);
     121    psLogMsg("psphot", PS_LOG_INFO, "%s", psphot);
     122
     123    psFree(pslib);
     124    psFree(psmodules);
     125    psFree(psphot);
     126
     127    return;
     128}

  • branches/cnb_branches/cnb_branch_20090301/psphot/src/psphotVisual.c

    r21519 r23352  
    1818// these are invoked by the -visual options
    1919
    20 static bool isVisual = false;
    2120static int kapa = -1;
    2221static int kapa2 = -1;
    2322static int kapa3 = -1;
    2423
    25 bool psphotSetVisual (bool mode) {
    26 
    27     isVisual = mode;
    28     return true;
    29 }
    3024
    3125bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) {
    … …  
    3832
    3933    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    40     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
     34    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    4135    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
    4236        fprintf (stderr, "failed to get background values\n");
    … …  
    8478
    8579    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    86     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
     80    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    8781    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
    8882        fprintf (stderr, "failed to get background values\n");
    … …  
    135129bool psphotVisualShowImage (pmReadout *readout) {
    136130
    137     if (!isVisual) return true;
     131    if (!pmVisualIsVisual()) return true;
    138132
    139133    if (kapa == -1) {
    … …  
    141135        if (kapa == -1) {
    142136            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    143             isVisual = false;
     137            pmVisualSetVisual(false);
    144138            return false;
    145139        }
    … …  
    164158    pmReadout *backgnd;
    165159
    166     if (!isVisual) return true;
     160    if (!pmVisualIsVisual()) return true;
    167161
    168162    if (kapa == -1) {
    … …  
    170164        if (kapa == -1) {
    171165            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    172             isVisual = false;
     166            pmVisualSetVisual(false);
    173167            return false;
    174168        }
    … …  
    199193bool psphotVisualShowSignificance (psImage *image) {
    200194
    201     if (!isVisual) return true;
     195    if (!pmVisualIsVisual()) return true;
    202196
    203197    if (kapa == -1) {
    … …  
    205199        if (kapa == -1) {
    206200            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    207             isVisual = false;
     201            pmVisualSetVisual(false);
    208202            return false;
    209203        }
    … …  
    228222    KiiOverlay *overlay;
    229223
    230     if (!isVisual) return true;
     224    if (!pmVisualIsVisual()) return true;
    231225
    232226    if (kapa == -1) {
    … …  
    306300    KiiOverlay *overlay;
    307301
    308     if (!isVisual) return true;
     302    if (!pmVisualIsVisual()) return true;
    309303
    310304    if (kapa == -1) {
    … …  
    423417    psEllipseAxes axes;
    424418
    425     if (!isVisual) return true;
     419    if (!pmVisualIsVisual()) return true;
    426420
    427421    if (kapa == -1) {
    … …  
    478472    Graphdata graphdata;
    479473
    480     if (!isVisual) return true;
     474    if (!pmVisualIsVisual()) return true;
    481475
    482476    if (kapa3 == -1) {
    … …  
    484478        if (kapa3 == -1) {
    485479            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    486             isVisual = false;
     480            pmVisualSetVisual(false);
    487481            return false;
    488482        }
    … …  
    694688bool psphotVisualShowRoughClass (psArray *sources) {
    695689
    696     if (!isVisual) return true;
     690    if (!pmVisualIsVisual()) return true;
    697691
    698692    if (kapa == -1) {
    … …  
    723717bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf) {
    724718
    725     if (!isVisual) return true;
     719    if (!pmVisualIsVisual()) return true;
    726720
    727721    if (kapa2 == -1) {
    … …  
    729723        if (kapa2 == -1) {
    730724            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    731             isVisual = false;
     725            pmVisualSetVisual(false);
    732726            return false;
    733727        }
    … …  
    801795    bool status;
    802796
    803     if (!isVisual) return true;
     797    if (!pmVisualIsVisual()) return true;
    804798
    805799    if (kapa2 == -1) {
    … …  
    807801        if (kapa2 == -1) {
    808802            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    809             isVisual = false;
     803            pmVisualSetVisual(false);
    810804            return false;
    811805        }
    … …  
    961955    bool status;
    962956
    963     if (!isVisual) return true;
     957    if (!pmVisualIsVisual()) return true;
    964958
    965959    if (kapa2 == -1) {
    … …  
    967961        if (kapa2 == -1) {
    968962            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    969             isVisual = false;
     963            pmVisualSetVisual(false);
    970964            return false;
    971965        }
    … …  
    12141208    KapaSection section;  // put the positive profile in one and the residuals in another?
    12151209
    1216     if (!isVisual) return true;
     1210    if (!pmVisualIsVisual()) return true;
    12171211
    12181212    if (kapa3 == -1) {
    … …  
    12201214        if (kapa3 == -1) {
    12211215            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1222             isVisual = false;
     1216            pmVisualSetVisual(false);
    12231217            return false;
    12241218        }
    … …  
    12861280    psEllipseAxes axes;
    12871281
    1288     if (!isVisual) return true;
     1282    if (!pmVisualIsVisual()) return true;
    12891283
    12901284    if (kapa == -1) {
    … …  
    13871381    KiiOverlay *overlay;
    13881382
    1389     if (!isVisual) return true;
     1383    if (!pmVisualIsVisual()) return true;
    13901384
    13911385    if (kapa == -1) {
    … …  
    14641458    KapaSection section;
    14651459
    1466     if (!isVisual) return true;
     1460    if (!pmVisualIsVisual()) return true;
    14671461
    14681462    if (kapa3 == -1) {
    … …  
    14701464        if (kapa3 == -1) {
    14711465            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1472             isVisual = false;
     1466            pmVisualSetVisual(false);
    14731467            return false;
    14741468        }
    … …  
    16041598bool psphotVisualShowResidualImage (pmReadout *readout) {
    16051599
    1606     if (!isVisual) return true;
     1600    if (!pmVisualIsVisual()) return true;
    16071601
    16081602    if (kapa == -1) {
    … …  
    16101604        if (kapa == -1) {
    16111605            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1612             isVisual = false;
     1606            pmVisualSetVisual(false);
    16131607            return false;
    16141608        }
    … …  
    16311625    Graphdata graphdata;
    16321626
    1633     if (!isVisual) return true;
     1627    if (!pmVisualIsVisual()) return true;
    16341628
    16351629    if (kapa3 == -1) {
    … …  
    16371631        if (kapa3 == -1) {
    16381632            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    1639             isVisual = false;
     1633            pmVisualSetVisual(false);
    16401634            return false;
    16411635        }
    … …  
    17121706# else
    17131707
    1714 bool psphotSetVisual (bool mode){}
    17151708bool psphotVisualShowImage (pmConfig *config, pmReadout *readout) { return true; }
    17161709bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout) { return true; }
Note: See TracChangeset for help on using the changeset viewer.