IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 8928


Ignore:
Timestamp:
Sep 24, 2006, 3:07:04 PM (20 years ago)
Author:
eugene
Message:

creating ppShutter tests, fixed ppShutter errors

Location:
trunk/ppShutter
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppShutter/src

    • Property svn:ignore

      •  

        old new  
        22Makefile
        33Makefile.in
        4 ppMerge
        54config.h
        65config.h.in
        76stamp-h1
         7ppShutter

  • trunk/ppShutter/src/.cvsignore

    r8901 r8928  
    22Makefile
    33Makefile.in
    4 ppMerge
    54config.h
    65config.h.in
    76stamp-h1
     7ppShutter

  • trunk/ppShutter/src/ppShutter.c

    r8901 r8928  
    55void usage () {
    66
    7     fprintf (stderr, "USAGE: ppShutter (output) (input) (input) (input) ...\n");
     7    fprintf (stderr, "USAGE: ppShutter (output) (outscale) (input) (input) (input) ...\n");
    88    fprintf (stderr, " minimum number of input images: 5\n");
    99    exit (1);
    … …  
    1919
    2020    // minimum number of input images: 5
    21     if (argc < 7) usage ();
     21    if (argc < 8) usage ();
    2222
    23     int Ninput = argc - 2;
     23    int Ninput = argc - 3;
    2424
    2525    // load the input images (for now, load them all into memory)
    … …  
    2828    psArray *images  = psArrayAlloc (Ninput);
    2929    for (int i = 0; i < Ninput; i++) {
    30         fits = psFitsOpen (argv[i-2], "r");
     30        fits = psFitsOpen (argv[i+3], "r");
    3131        files->data[i] = fits;
    3232
    … …  
    3737        images->data[i] = image;
    3838    }
     39    images->n = files->n = headers->n = Ninput;
    3940
    4041    // find the exposures times for each input images
    … …  
    4344    psVector *fluxRef = psVectorAlloc (Ninput, PS_TYPE_F32);
    4445    psVector *normCts = psVectorAlloc (Ninput, PS_TYPE_F32);
     46    psVector *normErr = psVectorAlloc (Ninput, PS_TYPE_F32);
    4547    for (int i = 0; i < Ninput; i++) {
    4648        exptime->data.F32[i] = psMetadataLookupF32 (NULL, headers->data[i], "EXPTIME");
    4749    }
     50    exptime->n = fluxRef->n = normCts->n = normErr->n = Ninput;
     51
    4852    psVector *index = psVectorSortIndex (NULL, exptime);
    4953
    … …  
    5256
    5357    // statistics used to measure initial parameters
    54     psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
     58    psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    5559
    5660    // define the reference region
    … …  
    6266        myStats = psImageStats(myStats, image, NULL, 0);
    6367        fluxRef->data.F32[i] = myStats->sampleMean;
     68        fprintf (stderr, "ref flux for image %d: %f\n", i, fluxRef->data.F32[i]);
    6469    }   
    6570
    … …  
    7984        psPlane *coord = samples->data[j];
    8085
     86        fprintf (stderr, "norm cnts for region %d: ", j);
    8187        for (int i = 0; i < Ninput; i++) {
    8288            psImage *image = images->data[index->data.S32[i]];
    … …  
    8490            region = psRegionForSquare(coord->x, coord->y, 10);
    8591            region = psRegionForImage(image, region);
     92            // fprintf (stderr, "%f,%f:%f,%f : ", region.x0, region.y0, region.x1, region.y1);
    8693
    87             myStats = psImageStats(myStats, image, NULL, 0);
     94            psImage *raster = psImageSubset (image, region);
     95            myStats = psImageStats(myStats, raster, NULL, 0);
    8896            normCts->data.F32[i] = myStats->sampleMean / fluxRef->data.F32[i];
     97            normErr->data.F32[i] = myStats->sampleStdev / fluxRef->data.F32[i];
     98            fprintf (stderr, "%f ", normCts->data.F32[i]);
     99            // fprintf (stderr, "%f\n", myStats->sampleMean);
     100           
     101            psFree (raster);
    89102        }       
     103        fprintf (stderr, "\n");
    90104       
    91105        pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, normCts);
    92         pmShutterCorrPars *params = pmShutterCorrectionFullFit (exptime, normCts, guess);
     106        pmShutterCorrPars *params = pmShutterCorrectionFullFit (exptime, normCts, normErr, guess);
    93107
    94108        refOffsetSample->data.F32[j] = params->offref;
    95109    }
     110    refOffsetSample->n = Nsamples;
    96111    myStats = psVectorStats (myStats, refOffsetSample, NULL, NULL, 0);
    97112    float refOffset = myStats->sampleMean;
     113    fprintf (stderr, "refOffset : %f +/- %f\n", myStats->sampleMean, myStats->sampleStdev);
    98114
    99115    // now run through all pixels in the images and measure the shutter parameters for each pixel
    … …  
    101117    psImage *offset = psImageAlloc (refimage->numCols, refimage->numRows, PS_TYPE_F32);
    102118
     119    int dotInterval = refimage->numRows / 32;
    103120    for (int ny = 0; ny < refimage->numRows; ny++) {
     121        if (ny % dotInterval == 0) fprintf (stderr, ".");       
    104122        for (int nx = 0; nx < refimage->numCols; nx++) {
    105123            for (int i = 0; i < Ninput; i++) {
    106124                psImage *image = images->data[index->data.S32[i]];
    107125                normCts->data.F32[i] = image->data.F32[ny][nx] / fluxRef->data.F32[i];
     126                normErr->data.F32[i] = sqrt(image->data.F32[ny][nx]) / fluxRef->data.F32[i];
    108127            }
    109             pmShutterCorrPars *params = pmShutterCorrectionLinFit (exptime, normCts, refOffset);
     128            pmShutterCorrPars *params = pmShutterCorrectionLinFit (exptime, normCts, normErr, refOffset);
    110129            pattern->data.F32[ny][nx] = params->scale;
    111130            offset->data.F32[ny][nx] = params->offset;
Note: See TracChangeset for help on using the changeset viewer.