IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 7, 2009, 4:08:25 PM (17 years ago)
Author:
Paul Price
Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:
branches/pap
Files:
30 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/pap

  • branches/pap/psphot

    • Property svn:mergeinfo deleted

  • branches/pap/psphot/src

    • Property svn:ignore

      •  

        old new  
        1717growthTest
        1818psphotVersionDefinitions.h
         19psphotMomentsStudy

  • branches/pap/psphot/src/Makefile.am

    r23808 r25027  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotTest
     27bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
    2828
    2929psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    3434psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    3535psphotTest_LDADD = libpsphot.la
     36
     37psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     38psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     39psphotMomentsStudy_LDADD = libpsphot.la
    3640
    3741psphot_SOURCES = \
    … …  
    5458psphotTest_SOURCES = \
    5559        psphotTest.c
     60
     61psphotMomentsStudy_SOURCES = \
     62        psphotMomentsStudy.c
    5663
    5764libpsphot_la_SOURCES = \

  • branches/pap/psphot/src/psphot.c

    r23195 r25027  
    11# include "psphotStandAlone.h"
    2 
    3 static void usage (void) {
    4     fprintf (stderr, "USAGE: psphot [-file image(s)] [-list imagelist] (output)\n");
    5     exit (PS_EXIT_CONFIG_ERROR);
    6 }
    72
    83int main (int argc, char **argv) {
    … …  
    149    // load command-line arguments, options, and system config data
    1510    pmConfig *config = psphotArguments (argc, argv);
    16     if (!config) {
    17         psErrorStackPrint(stderr, "Error reading arguments\n");
    18         usage ();
    19     }
     11    assert(config);
    2012
    2113    psphotVersionPrint();

  • branches/pap/psphot/src/psphot.h

    r23688 r25027  
    157157bool            psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal);
    158158bool            psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal);
    159 bool            psphotDeblendSatstars (psArray *sources, psMetadata *recipe);
     159bool            psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe);
    160160bool            psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first);
    161161

  • branches/pap/psphot/src/psphotApResid.c

    r21519 r25027  
    168168        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
    169169
     170        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
     171        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
     172           
    170173        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
    171174            continue;
    … …  
    317320*/
    318321
     322// XXX this still sucks...  need a better way to estimate the error floor...
    319323bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
    320324
    … …  
    355359
    356360        if (j > 2) {
    357             bool status = true;
    358             status &= psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff);
    359             status &= psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff);
    360             if (!status) { psErrorClear (); }
     361            if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) {
     362                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     363                return false;
     364            }
     365            if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) {
     366                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     367                return false;
     368            }
    361369            dSo->data.F32[i] = statsS->robustStdev;
    362370            dMo->data.F32[i] = statsM->sampleMean;
    363371            dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
    364             //      fprintf (stderr, "%d (%d) : sys: %f, phot: %f, rat: %f\n", i, j, dSo->data.F32[i], dMo->data.F32[i], dRo->data.F32[i]);
    365372        } else {
    366373            dSo->data.F32[i] = NAN;
    … …  
    375382    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
    376383    if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
    377         // XXX better testing of raised error
    378         psErrorClear();
     384        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
     385        return false;
    379386    }
    380387
    … …  
    382389    for (int i = 0; i < dSo->n; i++) {
    383390        *errorFloor = dSo->data.F32[i];
     391        if (fabs(*errorFloor) <= FLT_EPSILON) continue;
    384392        if (isfinite(*errorFloor)) break;
    385393    }

  • branches/pap/psphot/src/psphotArguments.c

    r23242 r25027  
    11# include "psphotStandAlone.h"
     2
     3static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile)
     4{
     5  fprintf(ofile,
     6          "Usage: one of the following\n"
     7          "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n"
     8          "     -variance varfile1[,varfile2,...] OutFileBaseName\n"
     9          "\n"
     10          "%s -list FileNameList [-masklist MaskFileNameList] \n"
     11          "     -variancelist VarFileNameList OutFileBaseName\n"
     12          "\n"
     13          "%s -help\n"
     14          "\n"
     15          "%s -version\n"
     16          "\n"
     17          "where:\n"
     18          "  FileNameList is a text file containing filenames, one per line\n"
     19          "  MaskFileNameList is a text file of mask filenames, one per line\n"
     20          "  VarFileNameList is a text file of variance filenames, one per line\n"
     21          "  OutFileBaseName is the 'root name' for output files\n"
     22          "\n"
     23          "additional options:\n"
     24          "  -modeltest xObj yObj [-model DEFAULT|ModelName] \n"
     25          "        [-fitmode DEFAULT|PSF|CONV] [-fitset FitFileName]\n"
     26          "     Test fit for object at the given coordinates.  ModelName\n"
     27          "     is one of PS_MODEL_GAUSS, PS_MODEL_PGAUSS, PS_MODEL_QGAUSS,\n"
     28          "     PS_MODEL_RGAUSS, PS_MODEL_PS1_V1, or PS_MODEL_SERSIC.\n"
     29          "     FitFileName is a file of x,y,Io triples.\n"
     30          "  -psf PsfFile1[,PsfFile2,...] or -psflist PsfFileNameList\n"
     31          "     specify PSF rather than letting %s estimate it\n"
     32          "  -src SrcFile1[,SrcFile2,...] or -srclist SrcFileNameList\n"
     33          "     specify additional sources for PSF generation\n"
     34          "  -chip nn[,nn,...]\n"
     35          "     select detector chips to process; default is all.\n"
     36          "     Indices correspond to zero-based offset in the FPA metadata table.\n"
     37          "  -photcode PhotoCodeName\n"
     38          "     specify photocode\n"
     39          "  -region RegionString\n"
     40          "     specify analysis region.  String is of form '[x0:x1,y0:y1]'\n"
     41          "     To use this option you must define a default in psphot.config\n"
     42          "  -visual\n"
     43          "     turns on interactive display mode\n"
     44          "  -dumpconfig CfgFileName\n"
     45          "     causes config info to be dumped to the named file.\n"
     46          "  -break NOTHING|BACKMDL|PEAKS|MOMENTS|PSFMODEL|ENSEMBLE|PASS1\n"
     47          "     choose a point at which to exit processing early\n"
     48          "  -nthreads n\n"
     49          "     set number of parallel threads of execution\n"
     50          "  -F OldFileRule ReplacementFileRule\n"
     51          "     change file naming rule; e.g. '-F PSPHOT.OUTPUT PSPHOT.OUT.CMF.MEF'\n"
     52          "  -D name stringval\n"
     53          "     set a string-valued config parameter\n"
     54          "  -Di name intval\n"
     55          "     set an integer-valued config parameter\n"
     56          "  -Df name fval\n"
     57          "     set a float-valued config parameter\n"
     58          "  -Db name boolval\n"
     59          "     set a boolean-valued config parameter\n"
     60          "  -v, -vv, -vvv\n"
     61          "     set increasing levels of verbosity\n"
     62          "  -logfmt FormatString\n"
     63          "     set format string used for log messages\n"
     64          "  -trace Fac Lvl\n"
     65          "     set tracing for facility Fac to integer Lvl, e.g. '-trace err 10'\n"
     66          "  -trace-levels\n"
     67          "     print current trace levels\n",
     68          program,program,program,program,program);
     69    psFree(config);
     70    pmConfigDone();
     71    psLibFinalize();
     72    exit(PS_EXIT_SUCCESS);
     73}
     74
     75static void usage(const char *program,  // Name of the program
     76                  psMetadata *arguments, // Command-line arguments
     77                  pmConfig *config,      // Configuration
     78                  int exitCode
     79                  )
     80{
     81  fprintf(stderr,
     82          "Usage: one of the following\n"
     83          "%s -file fname1[,fname2,...] -mask maskfile1[,maskfile2,...]\n"
     84          "     -variance varfile1[,varfile2,...] OutFileBaseName\n"
     85          "\n"
     86          "%s -list FileNameList [-masklist MaskFileNameList] \n"
     87          "     -variancelist VarFileNameList OutFileBaseName\n"
     88          "\n"
     89          "Try '%s -help' for more options and explanation\n",
     90          program,program,program);
     91    if (exitCode != PS_EXIT_SUCCESS)
     92      psErrorStackPrint(stderr, "Error reading arguments\n");
     93    psFree(config);
     94    pmConfigDone();
     95    psLibFinalize();
     96    exit(exitCode);
     97}
    298
    399pmConfig *psphotArguments(int argc, char **argv) {
    … …  
    5101    int N;
    6102    bool status;
    7 
    8     if (argc == 1) {
    9         psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc);
    10         return NULL;
    11     }
    12 
    13     if ((N = psArgumentGet (argc, argv, "-version"))) {
    14         psString version;
    15         version = psphotVersionLong();    fprintf (stdout, "%s\n", version); psFree (version);
    16         version = psModulesVersionLong(); fprintf (stdout, "%s\n", version); psFree (version);
    17         version = psLibVersionLong();     fprintf (stdout, "%s\n", version); psFree (version);
    18         exit (0);
    19     }
    20103
    21104    // load config data from default locations
    22105    pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE);
    23106    if (config == NULL) {
    24         psError(PSPHOT_ERR_CONFIG, false, "Can't read site configuration");
    25         return NULL;
    26     }
     107      psErrorStackPrint(stderr, "Can't read site configuration");
     108        exit(PS_EXIT_CONFIG_ERROR);
     109    }
     110
     111    PSARGUMENTS_INSTANTIATE_GENERICS( psphot, config, argc, argv );
    27112
    28113    // save the following additional recipe values based on command-line options
    … …  
    30115    psMetadata *options = pmConfigRecipeOptions (config, PSPHOT_RECIPE);
    31116
    32     // Number of threads
    33     if ((N = psArgumentGet(argc, argv, "-threads"))) {
    34         psArgumentRemove(N, &argc, argv);
    35         int nThreads = atoi(argv[N]);
    36         psMetadataAddS32(config->arguments, PS_LIST_TAIL, "NTHREADS", 0, "number of psphot threads", nThreads);
    37         psArgumentRemove(N, &argc, argv);
    38 
    39         // create the thread pool with number of desired threads, supplying our thread launcher function
    40         psThreadPoolInit (nThreads);
    41     }
    42     psphotSetThreads();
     117    // Number of threads is handled
     118    PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv )
    43119
    44120    // run the test model (requires X,Y coordinate)
    45121    if ((N = psArgumentGet (argc, argv, "-modeltest"))) {
     122        if (argc<=N+2) {
     123          psErrorStackPrint(stderr, "Expected to see 2 more arguments; saw %d", argc - 1);
     124          exit(PS_EXIT_CONFIG_ERROR);
     125        }
     126
    46127        psMetadataAddBool (options, PS_LIST_TAIL, "TEST_FIT",   0, "", true);
    47128        psMetadataAddF32  (options, PS_LIST_TAIL, "TEST_FIT_X", 0, "", atof(argv[N+1]));
    … …  
    54135        // specify the modeltest model
    55136        if ((N = psArgumentGet (argc, argv, "-model"))) {
     137            if (argc<=N+1) {
     138              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     139              exit(PS_EXIT_CONFIG_ERROR);
     140            }
    56141            psArgumentRemove (N, &argc, argv);
    57142            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODEL", 0, "", argv[N]);
    … …  
    61146        // specify the test fit mode
    62147        if ((N = psArgumentGet (argc, argv, "-fitmode"))) {
     148            if (argc<=N+1) {
     149              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     150              exit(PS_EXIT_CONFIG_ERROR);
     151            }
    63152            psArgumentRemove (N, &argc, argv);
    64153            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_MODE", 0, "", argv[N]);
    … …  
    66155        }
    67156        if ((N = psArgumentGet (argc, argv, "-fitset"))) {
     157            if (argc<=N+1) {
     158              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     159              exit(PS_EXIT_CONFIG_ERROR);
     160            }
    68161            psArgumentRemove (N, &argc, argv);
    69162            psMetadataAddStr (options, PS_LIST_TAIL, "TEST_FIT_SET", 0, "", argv[N]);
    … …  
    74167    // photcode : used in output to supplement header data (argument or recipe?)
    75168    if ((N = psArgumentGet (argc, argv, "-photcode"))) {
     169        if (argc<=N+1) {
     170          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     171          exit(PS_EXIT_CONFIG_ERROR);
     172        }
    76173        psArgumentRemove (N, &argc, argv);
    77174        psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]);
    … …  
    87184    // break : used from recipe throughout psphotReadout
    88185    if ((N = psArgumentGet (argc, argv, "-break"))) {
     186        if (argc<=N+1) {
     187          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     188          exit(PS_EXIT_CONFIG_ERROR);
     189        }
    89190        psArgumentRemove (N, &argc, argv);
    90191        psMetadataAddStr (options, PS_LIST_TAIL, "BREAK_POINT", PS_META_REPLACE, "", argv[N]);
    91         psArgumentRemove (N, &argc, argv);
    92     }
    93 
    94     // fitmode : used from recipe throughout psphotReadout
    95     if ((N = psArgumentGet (argc, argv, "-fitmode"))) {
    96         psArgumentRemove (N, &argc, argv);
    97         psMetadataAddStr (options, PS_LIST_TAIL, "FITMODE", PS_META_REPLACE, "", argv[N]);
    98192        psArgumentRemove (N, &argc, argv);
    99193    }
    … …  
    101195    // analysis region : overrides recipe value, used in psphotReadout/psphotEnsemblePSF
    102196    if ((N = psArgumentGet (argc, argv, "-region"))) {
     197        if (argc<=N+1) {
     198          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     199          exit(PS_EXIT_CONFIG_ERROR);
     200        }
    103201        psArgumentRemove (N, &argc, argv);
    104202        psMetadataAddStr (options, PS_LIST_TAIL, "ANALYSIS_REGION", 0, "", argv[N]);
    … …  
    108206    // chip selection is used to limit chips to be processed
    109207    if ((N = psArgumentGet (argc, argv, "-chip"))) {
     208        if (argc<=N+1) {
     209          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
     210          exit(PS_EXIT_CONFIG_ERROR);
     211        }
    110212        psArgumentRemove (N, &argc, argv);
    111213        psMetadataAddStr (config->arguments, PS_LIST_TAIL, "CHIP_SELECTIONS", PS_DATA_STRING, "", argv[N]);
    … …  
    121223    pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC",    "-src",    "-srclist");
    122224
     225    if (argc == 1) {
     226        psError(PSPHOT_ERR_ARGUMENTS, true, "Too few arguments: %d", argc);
     227        usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR);
     228    }
     229
     230    if (psArgumentGet(argc, argv, "-help")
     231        || psArgumentGet(argc, argv, "-h"))
     232      writeHelpInfo(argv[0], config, stdout);
     233     
    123234    // the input file is a required argument; if not found, we will exit
    124235    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
    125236    if (!status) {
    126237        psError(PSPHOT_ERR_ARGUMENTS, false, "pmConfigFileSetsMD failed to parse arguments");
    127         return NULL;
     238        usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR);
    128239    }
    129240
    130241    if (argc != 2) {
    131242        psError(PSPHOT_ERR_ARGUMENTS, true, "Expected to see one more argument; saw %d", argc - 1);
    132         return NULL;
     243        usage(argv[0], config->arguments, config, PS_EXIT_CONFIG_ERROR);
    133244    }
    134245

  • branches/pap/psphot/src/psphotBlendFit.c

    r21519 r25027  
    197197        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    198198
    199         // skip DBL second sources (ie, added by psphotFitBlob)
     199        // skip DBL second sources (ie, added by psphotFitBlob
    200200        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
    201201

  • branches/pap/psphot/src/psphotChoosePSF.c

    r21183 r25027  
    361361
    362362    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
    363     psVectorStats (stats, fwhmMajor, NULL, NULL, 0);
     363    if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) {
     364        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR");
     365        return false;
     366    }
     367
    364368    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
    365369    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
    … …  
    367371    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
    368372
    369     psVectorStats (stats, fwhmMinor, NULL, NULL, 0);
     373    if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
     374        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR");
     375        return false;
     376    }
    370377    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
    371378    psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);

  • branches/pap/psphot/src/psphotCleanup.c

    r23487 r25027  
    22
    33void psphotCleanup (pmConfig *config) {
     4
     5    // Dump configuration if requested
     6    bool status;
     7    psString dump_file = psMetadataLookupStr(&status, config->arguments, "DUMP_CONFIG");
     8    if (dump_file) {
     9        (void)pmConfigCamerasCull(config,NULL);
     10        (void)pmConfigRecipesCull(config,NULL);
     11        pmConfigDump(config, dump_file);
     12    }
    413
    514    psFree (config);

  • branches/pap/psphot/src/psphotCullPeaks.c

    r17516 r25027  
    1818        nsigma_min = 0;
    1919    }
     20    float fPadding = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_PAD");
     21    if (!status) {
     22        fPadding = 0;
     23    }
    2024    const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");
    21 
    22     return pmFootprintArrayCullPeaks(image, weight, footprints,
    23                                      nsigma_delta, nsigma_min*skyStdev);
    24 }
    25 
    26 
    27 /*
    28  * Cull an entire psArray of pmFootprints
    29  * XXX drop this intermediate level function?
    30  */
    31 psErrorCode
    32 pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint
    33                           const psImage *weight,        // corresponding variance image
    34                           psArray *footprints, // array of pmFootprints
    35                           const float nsigma_delta, // how many sigma above local background a peak
    36                                         // needs to be to survive
    37                           const float min_threshold) { // minimum permitted coll height
     25    const float min_threshold = nsigma_min*skyStdev;
     26   
    3827    for (int i = 0; i < footprints->n; i++) {
    3928        pmFootprint *fp = footprints->data[i];
    40         if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) {
     29        if (pmFootprintCullPeaks(image, weight, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) {
    4130            return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
    4231        }
    … …  
    4534    return PS_ERR_NONE;
    4635}
    47 
    48  /*
    49   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    50   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    51   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
    52   * starting point, discard the peak.
    53   */
    54 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
    55                                  const psImage *weight, // corresponding variance image
    56                                  pmFootprint *fp, // Footprint containing mortal peaks
    57                                  const float nsigma_delta, // how many sigma above local background a peak
    58                                         // needs to be to survive
    59                                  const float min_threshold) { // minimum permitted coll height
    60     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    61     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    62     assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    63     assert (fp != NULL);
    64 
    65     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    66         return PS_ERR_NONE;
    67     }
    68 
    69     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    70     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    71     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    72     const psImage *subImg = psImageSubset((psImage *)img, subRegion);
    73     const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    74     assert (subImg != NULL && subWt != NULL);
    75     //
    76     // We need a psArray of peaks brighter than the current peak.  We'll fake this
    77     // by reusing the fp->peaks but lying about n.
    78     //
    79     // We do this for efficiency (otherwise I'd need two peaks lists), and we are
    80     // rather too chummy with psArray in consequence.  But it works.
    81     //
    82     psArray *brightPeaks = psArrayAlloc(0);
    83     psFree(brightPeaks->data);
    84     brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
    85     //
    86     // The brightest peak is always safe; go through other peaks trying to cull them
    87     //
    88     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    89         const pmPeak *peak = fp->peaks->data[i];
    90         int x = peak->x - subImg->col0;
    91         int y = peak->y - subImg->row0;
    92         //
    93         // Find the level nsigma below the peak that must separate the peak
    94         // from any of its friends
    95         //
    96         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
    97         const float stdev = sqrt(subWt->data.F32[y][x]);
    98         float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
    99         if (isnan(threshold) || threshold < min_threshold) {
    100 #if 1     // min_threshold is assumed to be below the detection threshold,
    101           // so all the peaks are pmFootprint, and this isn't the brightest
    102             // XXX mark peak to be dropped
    103             (void)psArrayRemoveIndex(fp->peaks, i);
    104             i--;                        // we moved everything down one
    105             continue;
    106 #else
    107 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
    108             threshold = min_threshold;
    109 #endif
    110         }
    111 
    112         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    113         if (threshold > subImg->data.F32[y][x]) {
    114             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    115         }
    116 
    117         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    118         // perhaps this should alloc a single ID image above and pass it in to be set.
    119 
    120         const int peak_id = 1;          // the ID for the peak of interest
    121         brightPeaks->n = i;             // only stop at a peak brighter than we are
    122 
    123         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    124 
    125         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
    126         brightPeaks->n = 0;             // don't double free
    127         psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
    128         psFree(peakFootprint);
    129 
    130         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    131         // If so, the current peak is bogus; drop it.
    132         int j;
    133         for (j = 0; j < i; j++) {
    134             const pmPeak *peak2 = fp->peaks->data[j];
    135             int x2 = peak2->x - subImg->col0;
    136             int y2 = peak2->y - subImg->row0;
    137             const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
    138 
    139             if (peak2_id == peak_id) {  // There's a brighter peak within the footprint above
    140                 ;                       // threshold; so cull our initial peak
    141                 (void)psArrayRemoveIndex(fp->peaks, i);
    142                 i--;                    // we moved everything down one
    143                 break;
    144             }
    145         }
    146         if (j == i) {
    147             j++;
    148         }
    149 
    150         psFree(idImg);
    151     }
    152 
    153     brightPeaks->n = 0; psFree(brightPeaks);
    154     psFree((psImage *)subImg);
    155     psFree((psImage *)subWt);
    156 
    157     return PS_ERR_NONE;
    158 }
    159 
    160  /*
    161   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    162   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    163   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
    164   * starting point, discard the peak.
    165   */
    166 
    167 # define IN_PEAK 1
    168 psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
    169                                  const psImage *weight, // corresponding variance image
    170                                  pmFootprint *fp, // Footprint containing mortal peaks
    171                                  const float nsigma_delta, // how many sigma above local background a peak
    172                                  // needs to be to survive
    173                                  const float min_threshold) { // minimum permitted coll height
    174     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    175     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    176     assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    177     assert (fp != NULL);
    178 
    179     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    180         return PS_ERR_NONE;
    181     }
    182 
    183     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    184     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    185     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    186 
    187     psImage *subImg = psImageSubset((psImage *)img, subRegion);
    188     psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    189     assert (subImg != NULL && subWt != NULL);
    190 
    191     psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
    192 
    193     // We need a psArray of peaks brighter than the current peak. 
    194     // We reject peaks which either:
    195     // 1) are below the local threshold
    196     // 2) have a brighter peak within their threshold
    197 
    198     // allocate the full-sized array.  if the final array is much smaller, we can realloc
    199     // at that point.
    200     psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n);
    201     psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);
    202 
    203     // The brightest peak is always safe; go through other peaks trying to cull them
    204     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    205         const pmPeak *peak = fp->peaks->data[i];
    206         int x = peak->x - subImg->col0;
    207         int y = peak->y - subImg->row0;
    208         //
    209         // Find the level nsigma below the peak that must separate the peak
    210         // from any of its friends
    211         //
    212         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
    213         const float stdev = sqrt(subWt->data.F32[y][x]);
    214         float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
    215         if (isnan(threshold) || threshold < min_threshold) {
    216             // min_threshold is assumed to be below the detection threshold,
    217             // so all the peaks are pmFootprint, and this isn't the brightest
    218             continue;
    219         }
    220 
    221         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    222         if (threshold > subImg->data.F32[y][x]) {
    223             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    224         }
    225 
    226         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    227         // perhaps this should alloc a single ID image above and pass it in to be set.
    228 
    229         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    230 
    231         // at this point brightPeaks only has the peaks brighter than the current
    232         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
    233 
    234         // XXX need to supply the image here
    235         // we set the IDs to either 1 (in peak) or 0 (not in peak)
    236         pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
    237         psFree(peakFootprint);
    238 
    239         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    240         // If so, the current peak is bogus; drop it.
    241         bool keep = true;
    242         for (int j = 0; keep && (j < brightPeaks->n); j++) {
    243             const pmPeak *peak2 = fp->peaks->data[j];
    244             int x2 = peak2->x - subImg->col0;
    245             int y2 = peak2->y - subImg->row0;
    246             if (idImg->data.S32[y2][x2] == IN_PEAK)
    247                 // There's a brighter peak within the footprint above threshold; so cull our initial peak
    248                 keep = false;
    249         }
    250         if (!keep) continue;
    251 
    252         psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
    253     }
    254 
    255     psFree (fp->peaks);
    256     fp->peaks = brightPeaks;
    257 
    258     psFree(idImg);
    259     psFree(subImg);
    260     psFree(subWt);
    261 
    262     return PS_ERR_NONE;
    263 }
    264 

  • branches/pap/psphot/src/psphotDeblendSatstars.c

    r21519 r25027  
    11# include "psphotInternal.h"
    22
    3 bool psphotDeblendSatstars (psArray *sources, psMetadata *recipe) {
     3bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe) {
    44
    55    int N;
    … …  
    99
    1010    int Nblend = 0;
    11     float SAT_TEST_LEVEL = 50000;
    1211    float SAT_MIN_RADIUS = 5.0;
     12
     13    bool status;
     14    pmCell *cell = readout->parent;
     15    float SATURATION = 0.75*psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION");
     16    float SAT_TEST_LEVEL = 0.5*SATURATION;
    1317
    1418    // we need sources spatially-sorted to find overlaps
    … …  
    3236        // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
    3337        if (source->mode & PM_SOURCE_MODE_BLEND) continue;
    34         if (source->peak->flux < SAT_TEST_LEVEL) continue;
     38        if (source->peak->flux < SATURATION) continue;
    3539
    3640        // save these for reference below
    … …  
    4549        psVector *xVec = contour->data[0];
    4650        psVector *yVec = contour->data[1];
    47 
    48         // XXX should we filter based on the number of pixels in the contour?
     51        if (xVec->n < 5) continue;
    4952
    5053        // find the center of the contour (let's just use mid[x,y])
    … …  
    6265        int yCenter = 0.5*(yMin + yMax);
    6366        psFree (contour);
     67
     68        psAssert (xCenter >= source->pixels->col0, "invalid shift in object center");
     69        psAssert (xCenter <  source->pixels->col0 + source->pixels->numCols, "invalid shift in object center");
     70        psAssert (yCenter >= source->pixels->row0, "invalid shift in object center");
     71        psAssert (yCenter <  source->pixels->row0 + source->pixels->numRows, "invalid shift in object center");
    6472
    6573        // reset the peak for this source to the value of the center pixel

  • branches/pap/psphot/src/psphotDiagnosticPlots.c

    r21108 r25027  
    4040
    4141    psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MIN);
    42     psVectorStats (stats, values, NULL, NULL, 0);
     42    if (!psVectorStats (stats, values, NULL, NULL, 0)) {
     43        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for histogram");
     44        return false;
     45    }
    4346
    4447    psHistogram *histogram = psHistogramAlloc (stats->min, stats->max, 1000);
    … …  
    6366    psStatsInit (stats);
    6467    stats->options = PS_STAT_MAX | PS_STAT_MIN;
    65     psVectorStats (stats, histogram->nums, NULL, NULL, 0);
     68    if (!psVectorStats (stats, histogram->nums, NULL, NULL, 0)) {
     69        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for histogram");
     70        return false;
     71    }
    6672
    6773    // scale the plot to hold the histogram
    … …  
    100106    psStatsInit (stats);
    101107    stats->options = PS_STAT_MAX | PS_STAT_MIN;
    102     psVectorStats (stats, histogram->nums, NULL, NULL, 0);
     108    if (!psVectorStats (stats, histogram->nums, NULL, NULL, 0)) {
     109        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for histogram");
     110        return false;
     111    }
    103112
    104113    // scale the plot to hold the histogram

  • branches/pap/psphot/src/psphotFindFootprints.c

    r21366 r25027  
    5959
    6060        psArray *mergedFootprints = pmFootprintArraysMerge(detections->footprints, footprints, includePeaks);
    61         psLogMsg ("psphot", PS_LOG_MINUTIA, "merged %ld new footprints with %ld existing ones: %f sec\n", footprints->n, detections->footprints->n, psTimerMark ("psphot.footprints"));
     61        psLogMsg ("psphot", PS_LOG_MINUTIA, "merged %ld new footprints with %ld existing ones: %f sec\n", footprints->n, (detections->footprints ? detections->footprints->n : 0), psTimerMark ("psphot.footprints"));
    6262
    6363        psFree(footprints);

  • branches/pap/psphot/src/psphotGuessModels.c

    r21519 r25027  
    66// 2) loop over the sources once and associate them with their cell
    77// 3) define the threaded function to work with sources for a given cell
    8 
    9 // A guess for when the moments aren't available
    10 static pmModel *wildGuess(pmSource *source, // Source for which to guess
    11                           pmPSF *psf    // The point-spread function
    12     )
    13 {
    14     pmModel *model = pmModelAlloc(psf->type);
    15     psF32 *PAR = model->params->data.F32;
    16     PAR[PM_PAR_SKY]  = 0;
    17     // XXX get this from the image pixels
    18     PAR[PM_PAR_I0]   = source->peak->flux;
    19     PAR[PM_PAR_XPOS] = source->peak->xf;
    20     PAR[PM_PAR_YPOS] = source->peak->yf;
    21     return model;
    22 }
    238
    249// construct an initial PSF model for each object
    … …  
    8166            }
    8267            psFree(job);
    83 
    84 # if (0)               
    85                 if (!psphotGuessModel_Unthreaded (readout, cells->data[j], psf, maskVal, markVal)) {
    86                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    87                     return false;
    88                 }
    89 # endif
    9068        }
    9169
    … …  
    152130        nSrc ++;
    153131       
    154         // XXX if a source is faint, it will not have moments measured.
    155         // it must be modelled as a PSF.  In this case, we need to use
    156         // the peak centroid to get the coordinates and get the peak flux
    157         // from the image?
    158         pmModel *modelEXT;
    159         if (!source->moments) {
    160             modelEXT = wildGuess(source, psf);
     132        // the guess central intensity comes from the peak:
     133        float Io = source->peak->flux;
     134
     135        // We have two options to get a guess for the object position: the position from the
     136        // peak and the position from the moments.  Use the peak position if (a) there are no
     137        // moments and (b) the sources is not saturated
     138
     139        bool useMoments = false;
     140        useMoments = (source->mode & PM_SOURCE_MODE_SATSTAR);  // we only want to try if SATSTAR is set, but..
     141        useMoments = (useMoments && source->moments);          // can't if there are no moments
     142        useMoments = (useMoments && source->moments->nPixels); // can't if the moments were not measured
     143        useMoments = (useMoments && !(source->mode && PM_SOURCE_MODE_MOMENTS_FAILURE)); // can't if the moments failed...
     144
     145        float Xo, Yo;
     146        if (useMoments) {
     147            Xo = source->moments->Mx;
     148            Yo = source->moments->My;
    161149        } else {
    162             // use the source moments, etc to guess basic model parameters
    163             modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC X5
    164             if (!modelEXT) {
    165                 modelEXT = wildGuess(source, psf);
    166             }
    167             // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
    168             if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
    169                 modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
    170                 modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
    171             } else {
    172                 modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;
    173                 modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf;
    174             }
     150            Xo = source->peak->xf;
     151            Yo = source->peak->yf;
    175152        }
    176153
    177         // set PSF parameters for this model (apply 2D shape model)
    178         pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC X5
     154        // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
     155        pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
     156
    179157        if (modelPSF == NULL) {
    180             psWarning ("Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
    181                     source->peak->y, source->peak->x);
     158            psWarning ("Failed to determine PSF model at (%f,%f); trying image center", Xo, Yo);
    182159
    183             // Try the center of the image
    184             modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
    185             modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
    186             modelPSF = pmModelFromPSF (modelEXT, psf);
     160            float Xc = 0.5*readout->image->numCols;
     161            float Yc = 0.5*readout->image->numRows;
     162            pmModel *modelPSF = pmModelFromPSFforXY(psf, Xc, Yc, Io);
    187163            if (modelPSF == NULL) {
    188164                psError(PSPHOT_ERR_PSF, false, "Failed to determine PSF model at center of image");
    189                 psFree(modelEXT);
    190165                return false;
    191166            }
     167
     168            // Now set the object position at the expected location:
     169            modelPSF->params->data.F32[PM_PAR_XPOS] = Xo;
     170            modelPSF->params->data.F32[PM_PAR_YPOS] = Yo;
    192171            source->mode |= PM_SOURCE_MODE_BADPSF;
    193172        }
    194         psFree (modelEXT); // FREE (x3)
    195173
    196         // XXX need to define the guess flux?
    197174        // set the fit radius based on the object flux limit and the model
    198175        // this function affects the mask pixels
    … …  
    209186    return true;
    210187}
    211 
    212 # if (0)
    213 // construct models only for sources in the specified region
    214 bool psphotGuessModel_Unthreaded (pmReadout *readout, psArray *sources, pmPSF *psf, psImageMaskType maskVal, psImageMaskType markVal) {
    215 
    216     int nSrc = 0;
    217 
    218     for (int i = 0; i < sources->n; i++) {
    219         pmSource *source = sources->data[i];
    220 
    221         // XXXX this is just for a test: use this to mark sources for which the model is measured
    222         // check later that all are used.
    223         source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
    224 
    225         // skip non-astronomical objects (very likely defects)
    226         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    227         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    228         if (!source->peak) continue;
    229 
    230         nSrc ++;
    231        
    232         // XXX if a source is faint, it will not have moments measured.
    233         // it must be modelled as a PSF.  In this case, we need to use
    234         // the peak centroid to get the coordinates and get the peak flux
    235         // from the image?
    236         pmModel *modelEXT;
    237         if (!source->moments) {
    238             modelEXT = wildGuess(source, psf);
    239         } else {
    240             // use the source moments, etc to guess basic model parameters
    241             modelEXT = pmSourceModelGuess (source, psf->type); // ALLOC
    242             if (!modelEXT) {
    243                 modelEXT = wildGuess(source, psf);
    244             }
    245             // these valuse are set in pmSourceModelGuess, should this rule be in there as well?
    246             if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
    247                 modelEXT->params->data.F32[PM_PAR_XPOS] = source->moments->Mx;
    248                 modelEXT->params->data.F32[PM_PAR_YPOS] = source->moments->My;
    249             } else {
    250                 modelEXT->params->data.F32[PM_PAR_XPOS] = source->peak->xf;
    251                 modelEXT->params->data.F32[PM_PAR_YPOS] = source->peak->yf;
    252             }
    253         }
    254 
    255         // set PSF parameters for this model (apply 2D shape model)
    256         pmModel *modelPSF = pmModelFromPSF (modelEXT, psf); // ALLOC
    257         if (modelPSF == NULL) {
    258             psError(PSPHOT_ERR_PSF, false,
    259                     "Failed to determine PSF model at r,c = (%d,%d); trying centre of image",
    260                     source->peak->y, source->peak->x);
    261             //
    262             // Try the centre of the image
    263             //
    264             modelEXT->params->data.F32[PM_PAR_XPOS] = 0.5*readout->image->numCols;
    265             modelEXT->params->data.F32[PM_PAR_YPOS] = 0.5*readout->image->numRows;
    266             modelPSF = pmModelFromPSF (modelEXT, psf);
    267             if (modelPSF == NULL) {
    268                 psError(PSPHOT_ERR_PSF, false,
    269                         "Failed to determine PSF model at centre of image");
    270                 psFree(modelEXT);
    271                 return false;
    272             }
    273 
    274             source->mode |= PM_SOURCE_MODE_BADPSF;
    275         }
    276         psFree (modelEXT);
    277 
    278         // XXX need to define the guess flux?
    279         // set the fit radius based on the object flux limit and the model
    280         // this function affects the mask pixels
    281         psphotCheckRadiusPSF (readout, source, modelPSF, markVal);
    282 
    283         // set the source PSF model
    284         source->modelPSF = modelPSF;
    285         source->modelPSF->residuals = psf->residuals;
    286 
    287         pmSourceCacheModel (source, maskVal);
    288 
    289     }
    290 
    291     return true;
    292 }
    293 # endif

  • branches/pap/psphot/src/psphotImageLoop.c

    r23167 r25027  
    7575                }
    7676
     77                psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
     78                if (!pmReadoutMaskNonfinite(readout, maskSat)) {
     79                    psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
     80                    psFree(view);
     81                    return false;
     82                }
     83
    7784                // run the actual photometry analysis on this chip/cell/readout
    7885                if (!psphotReadout (config, view)) {

  • branches/pap/psphot/src/psphotImageQuality.c

    r20471 r25027  
    8888                     "Number of stars used for IQ measurements", M2->n);
    8989
     90// XXX make this a recipe option
    9091#if (USE_SAMPLE)
    9192    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);

  • branches/pap/psphot/src/psphotMakeResiduals.c

    r21366 r25027  
    200200            // measure the robust median to determine a baseline reference value
    201201            *fluxClip = *fluxClipDef;
    202             psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal);
    203             psErrorClear();             // clear (ignore) any outstanding errors
     202            if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) {
     203                psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
     204                return false;
     205            }
     206            if (isnan(fluxClip->robustMedian)) {
     207                resid->Ro->data.F32[oy][ox] = 0.0;
     208                resid->Rx->data.F32[oy][ox] = 0.0;
     209                resid->Ry->data.F32[oy][ox] = 0.0;
     210                resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
     211                continue;
     212            }
    204213
    205214            // mark input pixels which are more than N sigma from the median
    … …  
    220229                // measure the desired statistic on the unclipped pixels
    221230                *fluxStats = *fluxStatsDef;
    222                 psVectorStats (fluxStats, fluxes, NULL, fmasks, fmaskVal);
    223                 psErrorClear();         // clear (ignore) any outstanding errors
     231                if (!psVectorStats (fluxStats, fluxes, NULL, fmasks, fmaskVal)) {
     232                    psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats");
     233                    return false;
     234                }
    224235
    225236                resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption);
    226237                resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0;
    227                 //resid->variance->data.F32[oy][ox] = fluxStats->sampleStdev;
     238
     239                if (isnan(resid->Ro->data.F32[oy][ox])) {
     240                    resid->Ro->data.F32[oy][ox] = 0.0;
     241                    resid->Rx->data.F32[oy][ox] = 0.0;
     242                    resid->Ry->data.F32[oy][ox] = 0.0;
     243                    resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask;
     244                    continue;
     245                }
    228246
    229247                if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*fluxStats->sampleStdev/sqrt(nKeep)) {

  • branches/pap/psphot/src/psphotMaskReadout.c

    r23688 r25027  
    2727    }
    2828    if (!readout->variance) {
    29         if (!pmReadoutGenerateVariance(readout, true)) {
     29        if (!pmReadoutGenerateVariance(readout, NULL, true)) {
    3030            psError (PSPHOT_ERR_CONFIG, false, "trouble creating variance");
    3131            return false;

  • branches/pap/psphot/src/psphotModelBackground.c

    r23287 r25027  
    189189
    190190            if (gotX && gotY) {
    191                 psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
    192                 psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
     191                (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 6);
     192                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
    193193            } else {
    194                 psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
    195                 psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
     194                (void) psTraceSetLevel ("psLib.math.vectorFittedStats_v4", 0);
     195                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
    196196            }
    197197            # endif

  • branches/pap/psphot/src/psphotModelTest.c

    r21183 r25027  
    139139
    140140    // get the source moments
    141     status = pmSourceMoments (source, mRADIUS);
     141    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0);
    142142    if (!status) psAbort("psSourceMoments error");
    143143    source->peak->value = source->moments->Peak;

  • branches/pap/psphot/src/psphotReadout.c

    r23486 r25027  
    6161    // display the backsub and backgnd images
    6262    psphotVisualShowBackground (config, view, readout);
    63 
     63   
    6464    // run a single-model test if desired (exits from here if test is run)
    6565    psphotModelTest (config, view, recipe);
    … …  
    8787    // find blended neighbors of very saturated stars
    8888    // XXX merge this with Basic Deblend?
    89     psphotDeblendSatstars (sources, recipe);
     89    psphotDeblendSatstars (readout, sources, recipe);
    9090
    9191    // mark blended peaks PS_SOURCE_BLEND
    … …  
    236236    return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    237237}
     238

  • branches/pap/psphot/src/psphotReadoutCleanup.c

    r20769 r25027  
    99    if (psErrorCodeLast() == PSPHOT_ERR_DATA) {
    1010        psErrorStackPrint(stderr, "Error in the psphot readout analysis");
    11         psErrorClear();
    12     } 
    13     if (psErrorCodeLast() != PS_ERR_NONE) { 
    14         psFree (psf);
    15         psFree (sources);
    16         psFree (detections);
    17         return false;
     11        psErrorClear();
     12    }
     13    if (psErrorCodeLast() != PS_ERR_NONE) {
     14        psFree (psf);
     15        psFree (sources);
     16        psFree (detections);
     17        return false;
    1818    }
    1919
    … …  
    2222        if (!psphotPSFstats (readout, recipe, psf)) {
    2323            psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters");
    24             psFree (psf);
    25             psFree (sources);
    26             psFree (detections);
     24            psFree (psf);
     25            psFree (sources);
     26            psFree (detections);
    2727            return false;
    2828        }
    … …  
    3232        if (!psphotMomentsStats (readout, recipe, sources)) {
    3333            psError(PSPHOT_ERR_PROG, false, "Failed to measure Moment shape parameters");
    34             psFree (psf);
    35             psFree (sources);
    36             psFree (detections);
     34            psFree (psf);
     35            psFree (sources);
     36            psFree (detections);
    3737            return false;
    3838        }
    3939    }
     40
     41    // Check to see if the image quality was measured
     42    if (!psf) {
     43        bool mdok;                      // Status of MD lookup
     44        int nIQ = psMetadataLookupS32(&mdok, recipe, "IQ_NSTAR"); // Number of stars for IQ measurement
     45        if (!mdok || nIQ <= 0) {
     46            psError(PSPHOT_ERR_DATA, false, "Unable to measure image quality");
     47            psFree (psf);
     48            psFree (sources);
     49            psFree (detections);
     50            return false;
     51        }
     52    }
     53
    4054
    4155    // write NSTARS to the image header
    … …  
    4862    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER",  PS_DATA_METADATA, "header stats", header);
    4963    if (sources) {
    50         psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
     64        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
    5165    }
    5266    if (psf) {
    53         // save the psf for possible output.  if there was already an entry, it was loaded from external sources
    54         // the new one may have been updated or modified, so replace the existing entry.  We
    55         // are required to save it on the chip, but this will cause problems if we ever want to
    56         // run psphot on an unmosaiced image
    57         pmCell *cell = readout->parent;
    58         pmChip *chip = cell->parent;
    59         psMetadataAdd (chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,  "psphot psf", psf);
     67        // save the psf for possible output.  if there was already an entry, it was loaded from external sources
     68        // the new one may have been updated or modified, so replace the existing entry.  We
     69        // are required to save it on the chip, but this will cause problems if we ever want to
     70        // run psphot on an unmosaiced image
     71        pmCell *cell = readout->parent;
     72        pmChip *chip = cell->parent;
     73        psMetadataAdd (chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,  "psphot psf", psf);
    6074    }
    6175
    6276    if (psErrorCodeLast() != PS_ERR_NONE) {
    6377        psErrorStackPrint(stderr, "unexpected remaining errors");
    64         abort();
     78        abort();
    6579    }
    6680

  • branches/pap/psphot/src/psphotReadoutMinimal.c

    r23442 r25027  
    4646
    4747    // find the detections (by peak and/or footprint) in the image.
    48     pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
     48    pmDetections *detections = pmDetectionsAlloc(); // New detections; allocated to ensure pass=2
     49    detections = psphotFindDetections(detections, readout, recipe);
    4950    if (!detections) {
    5051        psLogMsg ("psphot", 3, "unable to find detections in this image");
    … …  
    5859    // find blended neighbors of very saturated stars
    5960    // XXX merge this with Basic Deblend?
    60     psphotDeblendSatstars (sources, recipe);
     61    psphotDeblendSatstars (readout, sources, recipe);
    6162
    6263    // mark blended peaks PS_SOURCE_BLEND

  • branches/pap/psphot/src/psphotRoughClass.c

    r21183 r25027  
    9191        return false;
    9292    }
    93     if (!psfClump.X || !psfClump.Y) {
     93    if (!psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) {
    9494        psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1);
    9595        return false;

  • branches/pap/psphot/src/psphotSourceFits.c

    r21519 r25027  
    224224
    225225    // recalculate the source moments using the larger extended-source moments radius
    226     if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false;
     226    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
     227    if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false;
    227228
    228229    psTrace ("psphot", 5, "trying blob...\n");
    … …  
    234235    pmSource *tmpSrc = pmSourceAlloc ();
    235236
     237    // XXX need to handle failures better here
    236238    pmModel *EXT = psphotFitEXT (readout, source, modelTypeEXT, maskVal, markVal);
    237239    okEXT = psphotEvalEXT (tmpSrc, EXT);
    238     chiEXT = EXT->chisq / EXT->nDOF;
    239 
     240    chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
     241
     242    // DBL will always be defined, but DBL->data[n] might not
    240243    psArray *DBL = psphotFitDBL (readout, source, maskVal, markVal);
    241244    okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
    … …  
    244247
    245248    // correct first model chisqs for flux trend
     249    chiDBL = NAN;
    246250    ONE = DBL->data[0];
    247     chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    248     ONE->chisqNorm = ONE->chisq / chiTrend;
    249 
    250     // save chisq for double-star/galaxy comparison
    251     chiDBL = ONE->chisq / ONE->nDOF;
     251    if (ONE) {
     252      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
     253      ONE->chisqNorm = ONE->chisq / chiTrend;
     254      chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
     255    }
    252256
    253257    // correct second model chisqs for flux trend
    254258    ONE = DBL->data[1];
    255     chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    256     ONE->chisqNorm = ONE->chisq / chiTrend;
     259    if (ONE) {
     260      chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
     261      ONE->chisqNorm = ONE->chisq / chiTrend;
     262    }
    257263
    258264    psFree (tmpSrc);
    259265
    260     // psTraceSetLevel("psModules.objects.pmSourceFitSet", 0);
     266    // (void) psTraceSetLevel("psModules.objects.pmSourceFitSet", 0);
    261267
    262268    if (okEXT && okDBL) {

  • branches/pap/psphot/src/psphotSourceStats.c

    r23219 r25027  
    8989
    9090            if (!psThreadJobAddPending(job)) {
    91                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     91                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    9292                psFree (job);
    9393                return NULL;
    9494            }
    9595            psFree(job);
    96 
    97 # if (0)
    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;
    106 # endif
    10796        }
    10897
    10998        // wait for the threads to finish and manage results
    11099        if (!psThreadPoolWait (false)) {
    111             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     100            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
    112101            return NULL;
    113102        }
    … …  
    149138    psMetadata *recipe              = job->args->data[1];
    150139
    151     float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    152     if (!status) return false;
    153     float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    154     if (!status) return false;
    155     float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     140    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     141    if (!status) return false;
     142    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     143    if (!status) return false;
     144    float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     145    if (!status) return false;
     146    float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN");
     147    if (!status) return false;
     148    float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
    156149    if (!status) return false;
    157150
    … …  
    202195
    203196        // measure basic source moments
    204         status = pmSourceMoments (source, RADIUS);
     197        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
    205198        if (status) {
    206199            Nmoments ++;
    … …  
    212205        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    213206        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    214         status = pmSourceMoments (source, BIG_RADIUS);
     207        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
    215208        if (status) {
    216209            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    … …  
    292285
    293286        // measure basic source moments
    294         status = pmSourceMoments (source, RADIUS);
     287        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
    295288        if (status) {
    296289            Nmoments ++;
    … …  
    302295        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    303296        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    304         status = pmSourceMoments (source, BIG_RADIUS);
     297        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
    305298        if (status) {
    306299            Nmoments ++;

  • branches/pap/psphot/src/psphotTest.c

    r21355 r25027  
    5050    }
    5151
    52     psTraceSetLevel ("psLib.sys.mutex", 3);
     52    (void) psTraceSetLevel ("psLib.sys.mutex", 3);
    5353
    5454    int nThreads = atoi (argv[2]);

  • branches/pap/psphot/src/psphotVisual.c

    r23259 r25027  
    471471    bool status;
    472472    Graphdata graphdata;
     473    KapaSection section;
    473474
    474475    if (!pmVisualIsVisual()) return true;
    … …  
    485486    KapaClearPlots (kapa3);
    486487    KapaInitGraph (&graphdata);
    487 
    488     // there are N regions: use the first (guaranteed to exist) to get the overall limits
    489     psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    490 
    491     float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    492     float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    493 
    494     // examine sources to set data range
    495     graphdata.xmin = -0.05;
    496     graphdata.ymin = -0.05;
    497     graphdata.xmax = 2.0*psfX;
    498     graphdata.ymax = 2.0*psfY;
    499     KapaSetLimits (kapa3, &graphdata);
    500 
    501     KapaSetFont (kapa3, "helvetica", 14);
    502     KapaBox (kapa3, &graphdata);
    503     KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    504     KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
    505 
     488    KapaSetFont (kapa3, "courier", 14);
     489
     490    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
     491
     492    // select the max psfX,Y values for the plot limits
     493    float Xmin = 0.0, Xmax = 0.0;
     494    float Ymin = 0.0, Ymax = 0.0;
     495    {
     496        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     497        for (int n = 0; n < nRegions; n++) {
     498
     499            char regionName[64];
     500            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     501            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     502
     503            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     504            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     505            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     506            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     507
     508            float X0 = psfX - 4.0*psfdX;
     509            float X1 = psfX + 4.0*psfdX;
     510            float Y0 = psfY - 4.0*psfdY;
     511            float Y1 = psfY + 4.0*psfdY;
     512
     513            if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
     514            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     515            if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
     516            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     517        }
     518    }
     519
     520    // storage vectors for data to be plotted
    506521    psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    507522    psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     523    psVector *mBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    508524    psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    509525    psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     526    psVector *mFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    510527
    511528    // construct the vectors
    … …  
    519536        xFaint->data.F32[nF] = source->moments->Mxx;
    520537        yFaint->data.F32[nF] = source->moments->Myy;
     538        mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum);
    521539        nF++;
    522540
    523541        // XXX make this a user-defined cutoff
    524         if (source->moments->SN < 50)
     542        if (source->moments->SN < SN_LIM)
    525543            continue;
    526544
    527545        xBright->data.F32[nB] = source->moments->Mxx;
    528546        yBright->data.F32[nB] = source->moments->Myy;
     547        mBright->data.F32[nB] = -2.5*log10(source->moments->Sum);
    529548        nB++;
    530549    }
    531550    xFaint->n = nF;
    532551    yFaint->n = nF;
     552    mFaint->n = nF;
    533553
    534554    xBright->n = nB;
    535555    yBright->n = nB;
     556    mBright->n = nB;
     557
     558    // three sections: MxxMyy, MagMxx, MagMyy
     559
     560    // first section: MxxMyy
     561    section.dx = 0.75;
     562    section.dy = 0.75;
     563    section.x  = 0.00;
     564    section.y  = 0.00;
     565    section.name = psStringCopy ("MxxMyy");
     566    KapaSetSection (kapa3, &section);
     567    psFree (section.name);
     568
     569    graphdata.color = KapaColorByName ("black");
     570    graphdata.xmin = Xmin;
     571    graphdata.ymin = Ymin;
     572    graphdata.xmax = Xmax;
     573    graphdata.ymax = Ymax;
     574    KapaSetLimits (kapa3, &graphdata);
     575
     576    KapaBox (kapa3, &graphdata);
     577    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
     578    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
    536579
    537580    graphdata.color = KapaColorByName ("black");
    … …  
    551594    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
    552595
    553     // XXX draw N circles to outline the clumps
     596    // second section: MagMyy
     597    section.dx = 0.75;
     598    section.dy = 0.25;
     599    section.x  = 0.00;
     600    section.y  = 0.80;
     601    section.name = psStringCopy ("MagMyy");
     602    KapaSetSection (kapa3, &section);
     603    psFree (section.name);
     604
     605    graphdata.color = KapaColorByName ("black");
     606    graphdata.xmin = -17.1;
     607    graphdata.xmax =  -7.9;
     608    graphdata.ymin = Ymin;
     609    graphdata.ymax = Ymax;
     610    KapaSetLimits (kapa3, &graphdata);
     611
     612    strcpy (graphdata.labels, "0210");
     613    KapaBox (kapa3, &graphdata);
     614    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP);
     615    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
     616
     617    graphdata.color = KapaColorByName ("black");
     618    graphdata.ptype = 0;
     619    graphdata.size = 0.3;
     620    graphdata.style = 2;
     621    KapaPrepPlot (kapa3, nF, &graphdata);
     622    KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
     623    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     624
     625    graphdata.color = KapaColorByName ("red");
     626    graphdata.ptype = 2;
     627    graphdata.size = 0.5;
     628    graphdata.style = 2;
     629    KapaPrepPlot (kapa3, nB, &graphdata);
     630    KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
     631    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     632
     633    // third section: MagMxx
     634    section.dx = 0.25;
     635    section.dy = 0.75;
     636    section.x  = 0.80;
     637    section.y  = 0.00;
     638    section.name = psStringCopy ("MagMxx");
     639    KapaSetSection (kapa3, &section);
     640    psFree (section.name);
     641
     642    graphdata.color = KapaColorByName ("black");
     643    graphdata.xmin = Xmin;
     644    graphdata.xmax = Xmax;
     645    graphdata.ymin =  -7.9;
     646    graphdata.ymax = -17.1;
     647    KapaSetLimits (kapa3, &graphdata);
     648
     649    strcpy (graphdata.labels, "2001");
     650    KapaBox (kapa3, &graphdata);
     651    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
     652    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP);
     653
     654    graphdata.color = KapaColorByName ("black");
     655    graphdata.ptype = 0;
     656    graphdata.size = 0.3;
     657    graphdata.style = 2;
     658    KapaPrepPlot (kapa3, nF, &graphdata);
     659    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
     660    KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
     661
     662    graphdata.color = KapaColorByName ("red");
     663    graphdata.ptype = 2;
     664    graphdata.size = 0.5;
     665    graphdata.style = 2;
     666    KapaPrepPlot (kapa3, nB, &graphdata);
     667    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
     668    KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
     669
     670    // draw N circles to outline the clumps
    554671    {
     672        KapaSelectSection (kapa3, "MxxMyy");
     673
    555674        // draw a circle centered on psfX,Y with size of the psf limit
    556675        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
    … …  
    562681        graphdata.color = KapaColorByName ("blue");
    563682        graphdata.style = 0;
     683
     684        graphdata.xmin = Xmin;
     685        graphdata.ymin = Ymin;
     686        graphdata.xmax = Xmax;
     687        graphdata.ymax = Ymax;
     688        KapaSetLimits (kapa3, &graphdata);
    564689
    565690        for (int n = 0; n < nRegions; n++) {
    … …  
    626751    psFree (xBright);
    627752    psFree (yBright);
     753    psFree (mBright);
    628754    psFree (xFaint);
    629755    psFree (yFaint);
     756    psFree (mFaint);
    630757
    631758    // pause and wait for user input:
Note: See TracChangeset for help on using the changeset viewer.