Index: trunk/ppShutter/src/.cvsignore =================================================================== --- trunk/ppShutter/src/.cvsignore (revision 8928) +++ (revision ) @@ -1,7 +1,0 @@ -.deps -Makefile -Makefile.in -config.h -config.h.in -stamp-h1 -ppShutter Index: trunk/ppShutter/src/Makefile.am =================================================================== --- trunk/ppShutter/src/Makefile.am (revision 8928) +++ (revision ) @@ -1,19 +1,0 @@ -bin_PROGRAMS = ppShutter - -ppShutter_CFLAGS += $(PPSTATS_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -ppShutter_LDFLAGS += $(PPSTATS_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -Wl,-Bdynamic - -ppShutter_SOURCES = \ - ppShutter.c - -#noinst_HEADERS = \ -# ppShutter.h - -CLEANFILES = *~ - -clean-local: - -rm -f TAGS - -# Tags for emacs -tags: - etags `find . -name \*.[ch] -print` Index: trunk/ppShutter/src/ppShutter.c =================================================================== --- trunk/ppShutter/src/ppShutter.c (revision 8928) +++ (revision ) @@ -1,145 +1,0 @@ -#include -#include -#include - -void usage () { - - fprintf (stderr, "USAGE: ppShutter (output) (outscale) (input) (input) (input) ...\n"); - fprintf (stderr, " minimum number of input images: 5\n"); - exit (1); -} - -// simple test: ppShutter (output) (input) (input) (input) ... -int main(int argc, char **argv) -{ - - psFits *fits; - psRegion fullimage = {0, 0, 0, 0}; - psRegion region; - - // minimum number of input images: 5 - if (argc < 8) usage (); - - int Ninput = argc - 3; - - // load the input images (for now, load them all into memory) - psArray *files = psArrayAlloc (Ninput); - psArray *headers = psArrayAlloc (Ninput); - psArray *images = psArrayAlloc (Ninput); - for (int i = 0; i < Ninput; i++) { - fits = psFitsOpen (argv[i+3], "r"); - files->data[i] = fits; - - psMetadata *header = psFitsReadHeader (NULL, fits); - headers->data[i] = header; - - psImage *image = psFitsReadImage (fits, fullimage, 0); - images->data[i] = image; - } - images->n = files->n = headers->n = Ninput; - - // find the exposures times for each input images - // need to resort the array based on the exposure times - psVector *exptime = psVectorAlloc (Ninput, PS_TYPE_F32); - psVector *fluxRef = psVectorAlloc (Ninput, PS_TYPE_F32); - psVector *normCts = psVectorAlloc (Ninput, PS_TYPE_F32); - psVector *normErr = psVectorAlloc (Ninput, PS_TYPE_F32); - for (int i = 0; i < Ninput; i++) { - exptime->data.F32[i] = psMetadataLookupF32 (NULL, headers->data[i], "EXPTIME"); - } - exptime->n = fluxRef->n = normCts->n = normErr->n = Ninput; - - psVector *index = psVectorSortIndex (NULL, exptime); - - // use the first image as a coord reference - psImage *refimage = images->data[0]; - - // statistics used to measure initial parameters - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - - // define the reference region - region = psRegionForSquare(0.5*refimage->numCols, 0.5*refimage->numRows, 10); - - // first, measure the reference region statistic from each image - for (int i = 0; i < Ninput; i++) { - psImage *image = images->data[index->data.S32[i]]; - myStats = psImageStats(myStats, image, NULL, 0); - fluxRef->data.F32[i] = myStats->sampleMean; - fprintf (stderr, "ref flux for image %d: %f\n", i, fluxRef->data.F32[i]); - } - - // define the sample regions - int Nsamples = 4; - psArray *samples = psArrayAlloc (Nsamples); - for (int i = 0; i < 4; i++) { - psPlane *coord = psPlaneAlloc (); - coord->x = (i % 2) ? 10 : refimage->numCols - 10; - coord->y = (i > 1) ? 10 : refimage->numRows - 10; - samples->data[i] = coord; - } - - // use several regions to make a measurement of dTo - psVector *refOffsetSample = psVectorAlloc (Nsamples, PS_TYPE_F32); - for (int j = 0; j < Nsamples; j++) { - psPlane *coord = samples->data[j]; - - fprintf (stderr, "norm cnts for region %d: ", j); - for (int i = 0; i < Ninput; i++) { - psImage *image = images->data[index->data.S32[i]]; - - region = psRegionForSquare(coord->x, coord->y, 10); - region = psRegionForImage(image, region); - // fprintf (stderr, "%f,%f:%f,%f : ", region.x0, region.y0, region.x1, region.y1); - - psImage *raster = psImageSubset (image, region); - myStats = psImageStats(myStats, raster, NULL, 0); - normCts->data.F32[i] = myStats->sampleMean / fluxRef->data.F32[i]; - normErr->data.F32[i] = myStats->sampleStdev / fluxRef->data.F32[i]; - fprintf (stderr, "%f ", normCts->data.F32[i]); - // fprintf (stderr, "%f\n", myStats->sampleMean); - - psFree (raster); - } - fprintf (stderr, "\n"); - - pmShutterCorrPars *guess = pmShutterCorrectionGuess (exptime, normCts); - pmShutterCorrPars *params = pmShutterCorrectionFullFit (exptime, normCts, normErr, guess); - - refOffsetSample->data.F32[j] = params->offref; - } - refOffsetSample->n = Nsamples; - myStats = psVectorStats (myStats, refOffsetSample, NULL, NULL, 0); - float refOffset = myStats->sampleMean; - fprintf (stderr, "refOffset : %f +/- %f\n", myStats->sampleMean, myStats->sampleStdev); - - // now run through all pixels in the images and measure the shutter parameters for each pixel - psImage *pattern = psImageAlloc (refimage->numCols, refimage->numRows, PS_TYPE_F32); - psImage *offset = psImageAlloc (refimage->numCols, refimage->numRows, PS_TYPE_F32); - - int dotInterval = refimage->numRows / 32; - for (int ny = 0; ny < refimage->numRows; ny++) { - if (ny % dotInterval == 0) fprintf (stderr, "."); - for (int nx = 0; nx < refimage->numCols; nx++) { - for (int i = 0; i < Ninput; i++) { - psImage *image = images->data[index->data.S32[i]]; - normCts->data.F32[i] = image->data.F32[ny][nx] / fluxRef->data.F32[i]; - normErr->data.F32[i] = sqrt(image->data.F32[ny][nx]) / fluxRef->data.F32[i]; - } - pmShutterCorrPars *params = pmShutterCorrectionLinFit (exptime, normCts, normErr, refOffset); - pattern->data.F32[ny][nx] = params->scale; - offset->data.F32[ny][nx] = params->offset; - } - } - - // save the resulting offset image - fits = psFitsOpen (argv[1], "w"); - psFitsWriteImage (fits, NULL, offset, 0, NULL); - psFitsClose (fits); - - // save the resulting illumination pattern image - fits = psFitsOpen (argv[2], "w"); - psFitsWriteImage (fits, NULL, pattern, 0, NULL); - psFitsClose (fits); - - exit (0); -}