IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 6851


Ignore:
Timestamp:
Apr 12, 2006, 8:59:15 PM (20 years ago)
Author:
eugene
Message:

tests on non-poisson errors, updates to image background, tests of systematic biases, updates to use pmFPAfile paradigm

Location:
trunk/psphot/src
Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/Makefile.am

    r6727 r6851  
    4141        psphotModelTest.c       \
    4242        psphotFitSet.c          \
     43        psphotWeightBias.c      \
    4344        psphotCleanup.c         
    4445

  • trunk/psphot/src/psphot.c

    r6727 r6851  
    33// XXX need a better structure for handling optional sequences
    44int main (int argc, char **argv) {
    5 
    6     // headtest (argc, argv);
    75
    86    psTimerStart ("complete");
    … …  
    2725    exit (0);
    2826}
    29 
    30 /** I/O test code
    31 
    32     psFits *fits = psFitsOpen (argv[1], "r");
    33     psMetadata *header = psFitsReadHeader (NULL, fits);
    34     psFitsClose (fits);
    35 
    36     psMetadata *new = psMetadataCopy (NULL, header);
    37     psMetadataConfigWrite (new, "test.cnf");
    38 
    39     fits = psFitsOpen ("test.fits", "w");
    40     psFitsWriteHeaderNotImage (fits, new);
    41     psFitsClose (fits);
    42 
    43     exit (0);
    44 
    45 **/

  • trunk/psphot/src/psphot.h

    r6753 r6851  
    5454bool            psphotDumpMoments (psMetadata *config, psArray *sources);
    5555psMetadata     *psphotDefineHeader (psMetadata *config);
     56bool            psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf);
    5657
    5758// PSF / DBL / EXT evaluation functions
    … …  
    6364bool            psphotInitLimitsPSF (psMetadata *config);
    6465bool            psphotInitLimitsEXT (psMetadata *config);
    65 bool            psphotFitBlend (pmReadout *readout, pmSource *source);
    66 bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources);
    67 bool            psphotFitPSF (pmReadout *readout, pmSource *source);
     66bool            psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf);
     67bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf);
     68bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf);
    6869pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source);
    6970psArray        *psphotFitDBL (pmReadout *readout, pmSource *source);
    7071
     72// functions to support simultaneous multi-source fitting
     73bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode);
     74
    7175// XXX these can probably be dropped:
    72 
    73 // functions to support simultaneous multi-source fitting
    74 bool            psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
    7576
    7677# if (0)

  • trunk/psphot/src/psphotArguments.c

    r6715 r6851  
    1010
    1111    int N;
     12    bool status;
    1213
    1314    if (*argc == 1) usage ();
    … …  
    118119    }
    119120
    120     // we load all input files onto a psArray, to be parsed later
    121     psArray *input = psArrayAlloc (16);
    122     input->n = 0;
    123 
    124     // load the list of filenames the supplied file (may be a glob: "file*.fits")
    125     if ((N = psArgumentGet (*argc, argv, "-file"))) {
    126         glob_t globList;
    127         psArgumentRemove (N, argc, argv);
    128         globList.gl_offs = 0;
    129         glob (argv[N], 0, NULL, &globList);
    130         for (int i = 0; i < globList.gl_pathc; i++) {
    131             char *filename = psStringCopy (globList.gl_pathv[i]);
    132             psArrayAdd (input, 16, filename);
    133             psFree (filename);
    134         }
    135         psArgumentRemove (N, argc, argv);
    136     }
    137 
    138     // load the list from the supplied text file
    139     if ((N = psArgumentGet (*argc, argv, "-list"))) {
    140         int nItems;
    141         char line[1024]; // XXX limits the list lines to 1024 chars
    142         char word[1024];
    143         char *filename;
    144 
    145         psArgumentRemove (N, argc, argv);
    146         FILE *f = fopen (argv[N], "r");
    147         if (f == NULL) {
    148             psAbort ("psphot", "unable to open specified list file");
    149         }
    150         while (fgets (line, 1024, f) != NULL) {
    151             nItems = sscanf (line, "%s", word);
    152             switch (nItems) {
    153               case 0:
    154                 break;
    155               case 1:
    156                 filename = psStringCopy (word);
    157                 psArrayAdd (input, 16, filename);
    158                 psFree (filename);
    159                 break;
    160               default:
    161                 // rigid format, no comments allowed?
    162                 psAbort ("psphot", "error parsing input list file");
    163                 break;
    164             }
    165         }
    166         psArgumentRemove (N, argc, argv);
    167     }           
    168     if (input->n < 1) usage ();
    169 
    170     // input list gets places as an array on the config->arguements list
    171     psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "INPUT",  PS_DATA_ARRAY, "", input);
    172     psFree (input);
     121    // the input file is a required argument; if not found, we will exit
     122    status = pmConfigFileSetsMD (config->arguments, argc, argv, "INPUT", "-file", "-list");
     123    if (!status) { usage ();}
    173124
    174125    if (*argc != 2) usage ();

  • trunk/psphot/src/psphotBlendFit.c

    r6753 r6851  
    5454       
    5555        // try fitting PSFs, then try extended sources
    56         if (psphotFitBlend (readout, source)) continue;
    57         if (psphotFitBlob (readout, source, sources)) continue;
     56        if (psphotFitBlend (readout, source, psf)) continue;
     57        if (psphotFitBlob (readout, source, sources, psf)) continue;
    5858
    5959        psTrace ("psphot.blend", 5, "failed fits\n");

  • trunk/psphot/src/psphotChoosePSF.c

    r6715 r6851  
    2727    int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS");
    2828    if (!status) NSTARS = PS_MIN (sources->n, 200);
     29
     30    // use poissonian errors or local-sky errors
     31    bool POISSON_ERRORS = psMetadataLookupBool (&status, config, "POISSON_ERRORS");
     32    if (!status) POISSON_ERRORS = true;
     33    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
    2934
    3035    stars = psArrayAlloc (sources->n);
    … …  
    6368        psMetadataItem *item = psListGetAndIncrement (iter);
    6469        modelName = item->data.V;
    65         models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
     70        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS);
    6671    }
    6772    psFree (iter);

  • trunk/psphot/src/psphotEvalPSF.c

    r6715 r6851  
    109109    dSX = model->dparams->data.F32[4];
    110110    dSY = model->dparams->data.F32[5];
    111     Chi = model->chisq / model->nDOF;
     111    Chi = model->chisqNorm / model->nDOF;
    112112
    113113    // swing of sigma_x,y in sigmas

  • trunk/psphot/src/psphotFitSet.c

    r6715 r6851  
    11# include "psphot.h"
    22
    3 bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, bool PSF) {
     3bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) {
    44
    55    double x, y, Io;
    … …  
    2424    }
    2525
    26     pmSourceFitSet (source, modelSet, PSF);
     26    pmSourceFitSet (source, modelSet, mode);
    2727
    2828    // write out positive object

  • trunk/psphot/src/psphotImageMedian.c

    r6753 r6851  
    11# include "psphot.h"
    2 double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax);
    3 void fsort (float *value, int N);
    4 
    5 // random number seed to select a fraction of the image pixels
    6 static psRandom *rnd;
    7 
    8 // use no more than MAX_SAMPLE_PIXELS pixels for each median box
    9 static int MAX_SAMPLE_PIXELS;
    102
    113// generate the median in NxN boxes, clipping heavily
    … …  
    157    bool status;
    168    psRegion region;
    17     psImage *model = NULL;
     9    int MAX_SAMPLE_PIXELS;
    1810
    1911    psTimerStart ("psphot");
    2012
     13    // select the appropriate recipe information
     14    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
     15
     16    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
     17    if (!status) MAX_SAMPLE_PIXELS = 1000;
     18    psImageClippedStatsInit(MAX_SAMPLE_PIXELS);
     19
     20    // subtract this amount extra from the sky
     21    float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
     22    if (!status) SKY_BIAS = 0;
     23
    2124    // find the currently selected readout
    22     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
    23     pmFPAfile  *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    24     pmReadout  *readout = pmFPAviewThisReadout (view, input->fpa);
     25    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    2526
    2627    psImage *image = readout->image;
    2728    psImage *mask  = readout->mask;
    28 
    29     rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
    30     MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
    31     if (!status) MAX_SAMPLE_PIXELS = 1000;
    3229
    3330    // dimensions of input & output image
    … …  
    4946    int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
    5047
    51     // select model pixels, from output background model file, or create
    52     model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
     48    // select model pixels (from output background model file, or create internal file)
     49    pmReadout *model = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKMDL");
     50    if (model == NULL) {
     51        // select model pixels, from output background model file, or create
     52        model = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
     53    } else {
     54        // replace the supplied image data with an image of the desired size
     55        psImageRecycle (model->image, nx, ny, PS_TYPE_F32);
     56    }
     57    psF32 **modelData = model->image->data.F32;
    5358
    5459    // measure clipped median for subimages
    … …  
    6368            psImage *submask = psImageSubset (mask, region);
    6469
    65             model->data.F32[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.50);
     70            // XXX the value of the upper and lower cuts probably should be studied...
     71            modelData[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.75) + SKY_BIAS;
    6672
    6773            psFree (subset);
    … …  
    7682    psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
    7783
    78     // linear interpolation to full-scale
     84    // XXX temporarily until we get the pmFPAfile output finished
     85    psphotSaveImage (NULL, model->image, "model.fits");
    7986
    8087    // select background pixels, from output background file, or create
    81     psImage *background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
    82 
    83     // XXX this code skips the initial pixels
    84     for (int Iy = 0; Iy < ny-1; Iy ++) {
    85         for (int Ix = 0; Ix < nx-1; Ix ++) {
    86 
    87             float V00 = model->data.F32[Iy+0][Ix+0];
    88             float V01 = model->data.F32[Iy+0][Ix+1];
    89             float V10 = model->data.F32[Iy+1][Ix+0];
    90             float V11 = model->data.F32[Iy+1][Ix+1];
    91 
    92             // a single binned pixel quad
    93             // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
    94             // corresponding to (Ix,Iy), (Ix+1,Iy+1)
    95             int Xs = (Ix + 1)*DX - dx;
    96             int Ys = (Iy + 1)*DY - dy;
    97             int Xe = Xs + DX;
    98             int Ye = Ys + DY;
    99 
    100             for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    101                 float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
    102                 float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
    103                 float dV = (Vxe - Vxs) / DX;
    104                 float V  = Vxs;
    105                 for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    106                     background->data.F32[iy][ix] = V;
    107                     V += dV;
    108                 }
    109             }
    110         }
     88    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
     89    if (background == NULL) {
     90        background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
     91    } else {
     92        // replace the supplied image data with an image of the desired size
     93        // XXX should not have to do this: it should be allocated correctly
     94        // psImageRecycle (background->image, Nx, Ny, PS_TYPE_F32);
    11195    }
    11296
    113     // side pixels
    114     int Xs = DX - dx;
    115     int Xe = nx*DX - dx;
    116     for (int Iy = 0; Iy < ny - 1; Iy++) {
     97    psF32 **backData = background->image->data.F32;
    11798
    118         int Ys = (Iy + 1)*DY - dy;
    119         int Ye = Ys + DY;
    120 
    121         // leading edge
    122         float V0 = model->data.F32[Iy+0][0];
    123         float V1 = model->data.F32[Iy+1][0];
    124         float dV = (V1 - V0) / DY;
    125         float V = V0;
    126         for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    127             for (int ix = 0; ix < Xs; ix++) {
    128                 background->data.F32[iy][ix] = V;
    129             }
    130             V += dV;
    131         }
    132 
    133         // trailing edge
    134         V0 = model->data.F32[Iy+0][nx-1];
    135         V1 = model->data.F32[Iy+1][nx-1];
    136         dV = (V1 - V0) / DY;
    137         V = V0;
    138         for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
    139             for (int ix = Xe; ix < Nx; ix++) {
    140                 background->data.F32[iy][ix] = V;
    141             }
    142             V += dV;
    143         }
    144     }
    145 
    146     // top and bottom pixels
    147     int Ys = DY - dy;
    148     int Ye = ny*DY - dy;
    149     for (int Ix = 0; Ix < nx - 1; Ix++) {
    150 
    151         int Xs = (Ix + 1)*DX - dx;
    152         int Xe = Xs + DX;
    153 
    154         // top edge
    155         float V0 = model->data.F32[0][Ix+0];
    156         float V1 = model->data.F32[0][Ix+1];
    157         float dV = (V1 - V0) / DX;
    158         for (int iy = 0; iy < Ys; iy++) {
    159             float V = V0;
    160             for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    161                 background->data.F32[iy][ix] = V;
    162                 V += dV;
    163             }
    164         }
    165 
    166         // bottom edge
    167         V0 = model->data.F32[ny-1][Ix+0];
    168         V1 = model->data.F32[ny-1][Ix+1];
    169         dV = (V1 - V0) / DX;
    170         for (int iy = Ye; iy < Ny; iy++) {
    171             float V = V0;
    172             for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
    173                 background->data.F32[iy][ix] = V;
    174                 V += dV;
    175             }
    176         }
    177     }
    178 
    179     // the four corners
    180     {
    181         float V;
    182         // 0,0
    183         V = model->data.F32[0][0];
    184         for (int iy = 0; iy < DY - dy; iy++) {
    185             for (int ix = 0; ix < DX - dx; ix++) {
    186                 background->data.F32[iy][ix] = V;
    187             }
    188         }
    189         // Nx,0
    190         V = model->data.F32[0][nx-1];
    191         for (int iy = 0; iy < DY - dy; iy++) {
    192             for (int ix = nx*DX - dx; ix < Nx; ix++) {
    193                 background->data.F32[iy][ix] = V;
    194             }
    195         }
    196         // 0,Ny
    197         V = model->data.F32[ny-1][0];
    198         for (int iy = ny*DY - dy; iy < Ny; iy++) {
    199             for (int ix = 0; ix < DX - dx; ix++) {
    200                 background->data.F32[iy][ix] = V;
    201             }
    202         }
    203         // Nx,Ny
    204         V = model->data.F32[nx-1][ny-1];
    205         for (int iy = ny*DY - dy; iy < Ny; iy++) {
    206             for (int ix = nx*DX - dx; ix < Nx; ix++) {
    207                 background->data.F32[iy][ix] = V;
    208             }
    209         }
    210     }   
     99    // linear interpolation to full-scale
     100    psImageUnbin (background->image, model->image, DX, DY, dx, dy);
    211101    psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
    212102
    213     psphotSaveImage (NULL, background, "back.fits");
     103    // XXX temporarily until we get the pmFPAfile output finished...
     104    psphotSaveImage (NULL, background->image, "back.fits");
    214105
    215106    // back-sub image pixels, from output background file (don't create if not requested)
    216     psImage *backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32);
     107    pmReadout *backSub = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKSUB");
    217108
    218     // subtract the background model
     109    // subtract the background model (save in backSub, if requested)
    219110    for (int j = 0; j < image->numRows; j++) {
    220111        for (int i = 0; i < image->numCols; i++) {
    221112            if (!mask->data.U8[j][i]) {
    222                 image->data.F32[j][i] -= background->data.F32[j][i];
     113                image->data.F32[j][i] -= backData[j][i];
    223114                if (backSub) {
    224                     backSub->data.F32[j][i] = image->data.F32[j][i];
     115                    backSub->image->data.F32[j][i] = image->data.F32[j][i];
    225116                }
    226117            }
    … …  
    229120
    230121    psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
    231     psFree (rnd);
     122    psImageClippedStatsCleanup();
    232123
    233     // it is safe to free all of these: either they are saved on config->files or are local
    234     psFree (model);
    235     psFree (background);
    236     psFree (backSub);
     124    // the pmReadout selected in this function are all view on entries in config->files
    237125    return true;
    238126}
    239 
    240 double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax) {
    241 
    242     double value;
    243     int nx = image->numCols;
    244     int ny = image->numRows;
    245    
    246     if (nx*ny <= 0) return 0.0;
    247 
    248     int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny);
    249     int Npixels = nx*ny;
    250 
    251     psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32);
    252 
    253     float min = values->data.F32[0];
    254     float max = values->data.F32[0];
    255 
    256     int n = 0;
    257     for (int i = 0; i < Nsubset; i++) {
    258         double frnd = psRandomUniform (rnd);
    259         int pixel = Npixels * frnd;
    260         int ix = pixel % nx;
    261         int iy = pixel / nx;
    262 
    263         if (mask->data.U8[iy][ix] & maskValue) continue;
    264 
    265         value = image->data.F32[iy][ix];
    266         min = PS_MIN (value, min);
    267         max = PS_MIN (value, max);
    268         values->data.F32[n] = value;
    269         n++;
    270     }
    271     values->n = n;
    272 
    273     int imin = fmin * n;
    274     int imax = fmax * n;
    275     int npts = imax - imin + 1;
    276 
    277     fsort (values->data.F32, n);
    278 
    279     value = 0;
    280     for (int i = imin; (i <= imax) && (i < n); i++) {
    281         value += values->data.F32[i];
    282     }
    283     value = value / npts;
    284    
    285     // XXX correct for selection bias??
    286 
    287     psFree (values);
    288     return value;
    289 }
    290 
    291 void fsort (float *value, int N) {
    292 
    293     int l,j,ir,i;
    294     float temp;
    295  
    296     if (N < 2) return;
    297     l = N >> 1;
    298     ir = N - 1;
    299     for (;;) {
    300         if (l > 0) {
    301             temp = value[--l];
    302         }
    303         else {
    304             temp = value[ir];
    305             value[ir] = value[0];
    306             if (--ir == 0) {
    307                 value[0] = temp;
    308                 return;
    309             }
    310         }
    311         i = l;
    312         j = (l << 1) + 1;
    313         while (j <= ir) {
    314             if (j < ir && value[j] < value[j+1]) ++j;
    315             if (temp < value[j]) {
    316                 value[i]=value[j];
    317                 j += (i=j) + 1;
    318             }
    319             else j = ir + 1;
    320         }
    321         value[i] = temp;
    322     }
    323 }
    324 

  • trunk/psphot/src/psphotModelTest.c

    r6727 r6851  
    1010    char name[64];
    1111    pmPSF *psf = NULL;
     12    pmSourceFitMode fitMode;
    1213
    1314    psMetadataItem *item  = NULL;
     15
     16    // use poissonian errors or local-sky errors
     17    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
     18    if (!status) POISSON_ERRORS = true;
     19    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
    1420
    1521    // run model fitting tests on a single source
    … …  
    1723
    1824    // what fitting mode to use?
    19     char *psfModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
     25    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
    2026    if (!status) {
    21         psfModeWord = DEFAULT_MODE;
     27        fitModeWord = DEFAULT_MODE;
    2228    }
    23     bool psfMode = !strcasecmp (psfModeWord, "PSF");
     29    fitMode = PM_SOURCE_FIT_EXT;
     30    if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
    2431
    25     // in psfMode, psf sets the model type
    26     if (psfMode) {
     32    // in fitMode, psf sets the model type
     33    if (fitMode == PM_SOURCE_FIT_PSF) {
    2734        char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
    2835        if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
    … …  
    128135    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
    129136
    130     if (psfMode) {
     137    if (fitMode == PM_SOURCE_FIT_PSF) {
    131138        pmModel *modelPSF = pmModelFromPSF (model, psf);
    132139        psFree (model);
    … …  
    155162    char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
    156163    if (status) {
    157         status = psphotFitSet (source, model, fitset, psfMode);
     164        status = psphotFitSet (source, model, fitset, fitMode);
    158165        exit (0);
    159166    }
    160167
    161     status = pmSourceFitModel (source, model, psfMode);
     168    status = pmSourceFitModel (source, model, fitMode);
    162169
    163170    // measure the source mags

  • trunk/psphot/src/psphotParseCamera.c

    r6727 r6851  
    44
    55    bool status;
    6     psFits *fits = NULL;
    7     psMetadata *phu = NULL;
    8     psMetadata *format = NULL;
    96
    10     // psphot is supplied with a list of input images (may be only one image)
    11     psArray *infiles = psMetadataLookupPtr(&status, config->arguments, "INPUT");
    12     if (!status) psAbort (__func__, "missing INPUT entry");
    13     if (infiles->n < 1) psAbort (__func__, "empty input list");
    14 
    15     // if no camera has been specified, use the first image as a template for the rest.
    16     if (config->camera == NULL) {
    17         fits = psFitsOpen (infiles->data[0], "r");
    18         phu = psFitsReadHeader (NULL, fits);
    19         format = pmConfigCameraFormatFromHeader (config, phu);
    20         psFitsClose (fits);
     7    status = false;
     8    pmFPAfile *input = pmFPAfileFromArgs (&status, config, "PSPHOT.INPUT", "INPUT");
     9    if (!status) {
     10        psAbort (__func__, "missing INPUT entry");
    2111    }
    22     // There's no point in continuing if we can't recognize what we've got
    23     if (! config->camera) {
    24         psErrorStackPrint(stderr, "Can't find camera configuration!\n");
    25         exit(EXIT_FAILURE);
    26     }
    27 
    28     // files
    29     config->files = psMetadataAlloc ();
    30 
    31     // build the template fpa, set up the basic view
    32     // XXX : 216 leaks in pmFPAConstruct / psFree (input)
    33     pmFPA *input = pmFPAConstruct (config->camera);
    34 
    35     // assign the I/O files (potentially) needed by psphot
    36     // the output file is just a view to the file on config->files
    37     pmFPAfile *file = pmFPAfileDefine (config->files, format, input, "PSPHOT.INPUT");
    38     if (!file) {
    39         exit(EXIT_FAILURE);
    40     }
    41 
    42     for (int i = 0; i < infiles->n; i++) {
    43         if (phu == NULL) {
    44             fits = psFitsOpen (infiles->data[i], "r");
    45             phu = psFitsReadHeader (NULL, fits);
    46             pmConfigValidateCameraFormat (format, phu);
    47             psFitsClose (fits);
    48         }
    49 
    50         // set the view to the corresponding entry for this phu
    51         pmFPAview *view = pmFPAAddSource (input, phu, format);
    52 
    53         // XXX is this the correct psMD to save the filename?
    54         char *name = pmFPAfileNameFromRule (file->filextra, file, view);
    55         psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]);
    56 
    57         psFree (view);
    58         psFree (name);
    59         psFree (phu);
    60         phu = NULL;
    61     }
    62 
    63     pmFPAfileDefine (config->files, format, input, "PSPHOT.OUTPUT");
    64     pmFPAfileDefine (config->files, format, input, "PSPHOT.RESID");
    65 
    66     // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_INPUT");
    67     // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_OUTPUT");
    68 
    69     // build the template fpa, set up the basic view
    70     // supply the backgnd with a different camera?
    71     // allow alternate format?
    72     // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.BACKSUB");
    73     // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.BACKGND");
    74     // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.BACKMDL");
    75     // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.PSF_SAMPLE");
    76 
    77     // psphot is supplied with the output name
    78     // this needs to be better: supply to all output (WRITE) files?
    79     char *output = psMetadataLookupPtr(&status, config->arguments, "OUTPUT");
    80 
    81     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.OUTPUT");
    82     if (!status) psAbort (__func__, "missing OUTPUT entry");
    83     psMetadataAddStr (file->names, PS_LIST_TAIL, "OUTPUT", 0, "", output);
    84 
    85     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.RESID");
    86     if (!status) psAbort (__func__, "missing OUTPUT entry");
    87     psMetadataAddStr (file->names, PS_LIST_TAIL, "OUTPUT", 0, "", output);
    8812
    8913    // recipe override values (command-line options):
    … …  
    9620    }
    9721    psFree (iter);
     22
     23    // these calls bind the I/O handle to the specified fpa
     24    pmFPAfileDefine (config->files, config->camera, input->fpa, "PSPHOT.OUTPUT");
     25    pmFPAfileDefine (config->files, config->camera, input->fpa, "PSPHOT.RESID");
     26
     27    // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_INPUT");
     28    // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_OUTPUT");
     29
     30    // these calls construct a new fpa for the I/O handle
     31    int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND_XBIN");
     32    if (!status) {DX = 64;}
     33    int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND_YBIN");
     34    if (!status) {DY = 64;}
     35    pmFPAfileFromFPA (config, input->fpa, DX, DY, "PSPHOT.BACKMDL");
     36    pmFPAfileFromFPA (config, input->fpa,  1,  1, "PSPHOT.BACKGND");
     37    pmFPAfileFromFPA (config, input->fpa,  1,  1, "PSPHOT.BACKSUB");
     38
     39    // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.PSF_SAMPLE");
     40
     41    // psphot is supplied with the output name
     42    // this needs to be better: supply to all output (WRITE) files?
     43    char *output = psMetadataLookupPtr(&status, config->arguments, "OUTPUT");
     44    pmFPAfileAddFileNames (config->files, "OUTPUT", output, PM_FPA_MODE_WRITE);
    9845
    9946    // set default recipe values here
    … …  
    12976    # endif
    13077
    131     psFree (input);
    132     psFree (format);
    133 
    13478    psTrace(__func__, 1, "Done with psphotParseCamera...\n");
    13579    return true;

  • trunk/psphot/src/psphotReadout.c

    r6753 r6851  
    99    bool         status;
    1010
     11    // select the current recipe
     12    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
     13
    1114    // find the currently selected readout
    12     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
    13     pmFPAfile  *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    14     pmReadout  *readout = pmFPAviewThisReadout (view, input->fpa);
     15    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    1516
    1617    // XXX does this need to invoke I/O?
     18    // XXX move this input the psphotImageLoop level?
    1719    pmReadoutSetWeights (readout);
     20    psphotSaveImage (NULL, readout->weight, "weight.fits");
    1821
    1922    // I have a valid mask, now mask in the analysis region of interest
    … …  
    4952    psphotBlendFit (readout, sources, recipe, psf);
    5053
     54    // psphotWeightBias (readout, sources, recipe, psf);
     55
    5156    // replace fitted sources
    5257    psphotReplaceUnfit (sources);
    … …  
    7580    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
    7681
     82    // remove internal pmFPAfiles, if created
     83    pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
     84    pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
     85
    7786    // free up the local copies of the data
    7887    psFree (psf);

  • trunk/psphot/src/psphotSkyReplace.c

    r6727 r6851  
    55bool psphotSkyReplace (pmConfig *config, pmFPAview *view) {
    66
    7     bool status;
     7    // find the currently selected readout
     8    pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
     9    if (readout == NULL) psAbort ("psphot", "input not defined");
    810
    9     // find the currently selected readout
    10     //  psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
    11     pmFPAfile *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
    12     pmReadout *readout = pmFPAviewThisReadout (view, input->fpa);
     11    // select background pixels, from output background file, or create
     12    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
     13    if (background == NULL) psAbort ("psphot", "background not defined");
    1314
    1415    // select the corresponding images
    15     psImage *image = readout->image;
    16     psImage *mask  = readout->mask;
    17 
    18     // select background pixels, from output background file, or create
    19     psImage *background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", 0, 0, PS_TYPE_F32);
    20     if (background == NULL) {
    21         return false;
    22     }
     16    psF32 **image = readout->image->data.F32;
     17    psU8  **mask  = readout->mask->data.U8;
     18    psF32 **back  = background->image->data.F32;
    2319
    2420    // replace the background model
    25     for (int j = 0; j < image->numRows; j++) {
    26         for (int i = 0; i < image->numCols; i++) {
    27             if (!mask->data.U8[j][i]) {
    28                 image->data.F32[j][i] += background->data.F32[j][i];
     21    for (int j = 0; j < readout->image->numRows; j++) {
     22        for (int i = 0; i < readout->image->numCols; i++) {
     23            if (!mask[j][i]) {
     24                image[j][i] += back[j][i];
    2925            }
    3026        }
    3127    }
    32    
    33     psFree (background);
    3428    return true;
    3529}

  • trunk/psphot/src/psphotSourceFits.c

    r6753 r6851  
    33// given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
    44
    5 bool psphotFitPSF (pmReadout *readout, pmSource *source) {
     5bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) {
    66
    77    float x, y;
     8    double chiTrend;
    89
    910    // save the PSF model from the Ensemble fit
    … …  
    1819    // fit PSF model (set/unset the pixel mask)
    1920    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
    20     pmSourceFitModel (source, PSF, true);
     21    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF);
    2122    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
     23   
     24    // correct model chisq for flux trend
     25    chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
     26    PSF->chisqNorm = PSF->chisq / chiTrend;
     27    // PSF->chisqNorm = PSF->chisq;
    2228
    2329    // does the PSF model succeed?
    … …  
    6066}
    6167
    62 bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources) {
     68bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) {
    6369
    6470    bool okEXT, okDBL;
    6571    float chiEXT, chiDBL;
     72    double chiTrend;
     73    pmModel *ONE = NULL;
    6674
    6775    // skip the source if we don't think it is extended
    … …  
    8694    okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
    8795    okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
    88     pmModel *ONE = DBL->data[0];
    89     chiDBL = ONE->chisq / ONE->nDOF;
     96
     97    // correct first model chisqs for flux trend
     98    ONE = DBL->data[0];
     99    chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
     100    ONE->chisqNorm = ONE->chisq / chiTrend;
     101    // ONE->chisqNorm = ONE->chisq;
     102
     103    // save chisq for double-star/galaxy comparison
     104    chiDBL = ONE->chisqNorm;
     105
     106    // correct second model chisqs for flux trend
     107    ONE = DBL->data[1];
     108    chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
     109    ONE->chisqNorm = ONE->chisq / chiTrend;
     110    // ONE->chisqNorm = ONE->chisq;
    90111
    91112    psFree (tmpSrc); // XXX should I keep / save the flags set in the eval functions?
    … …  
    93114    if (okEXT && okDBL) {
    94115        psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
    95         // XXX EAM : a bogus bias: need to examine this better
     116        // XXX EAM : a bogus bias: need to examine this better 
    96117        if (3*chiEXT > chiDBL) goto keepDBL;
    97118        goto keepEXT;
    … …  
    188209    y = source->moments->y;
    189210
    190     // fit EXT (not PSF) model (set/unset the pixel mask)
     211    // fit PSF model (set/unset the pixel mask)
    191212    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
    192     pmSourceFitSet (source, modelSet, true);
     213    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
    193214    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
    194215
    … …  
    210231    // fit EXT (not PSF) model (set/unset the pixel mask)
    211232    psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
    212     pmSourceFitModel (source, EXT, false);
     233    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT);
    213234    psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
    214235
    … …  
    216237}
    217238
    218 bool psphotFitBlend (pmReadout *readout, pmSource *source) {
     239bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) {
    219240
    220241    float x, y, dR;
    … …  
    222243    // if this source is not a possible blend, just fit as PSF
    223244    if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
    224         bool status = psphotFitPSF (readout, source);
     245        bool status = psphotFitPSF (readout, source, psf);
    225246        return status;
    226247    }
    … …  
    273294    // fit PSF model (set/unset the pixel mask)
    274295    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
    275     pmSourceFitSet (source, modelSet, true);
     296    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
    276297    psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
     298
     299    // correct model chisq for flux trend
     300    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
     301    PSF->chisqNorm = PSF->chisq / chiTrend;
     302    // PSF->chisqNorm = PSF->chisq;
    277303
    278304    // evaluate the blend objects, subtract if good, free otherwise
    … …  
    280306        pmSource *src = sourceSet->data[i];
    281307        pmModel *model = modelSet->data[i];
     308
     309        // correct model chisq for flux trend
     310        chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
     311        model->chisqNorm = model->chisq / chiTrend;
     312        // model->chisqNorm = model->chisq;
    282313
    283314        // if we skip this one, free the corresponding blend entry model
Note: See TracChangeset for help on using the changeset viewer.