Index: /branches/eam_branches/ipp-20121130/ppBackground/src/ppBackground.h =================================================================== --- /branches/eam_branches/ipp-20121130/ppBackground/src/ppBackground.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppBackground/src/ppBackground.h (revision 34838) @@ -37,4 +37,13 @@ ); +/// Determine the binning from the recipe if available. +psImageBinning *ppBackgroundBinningByRecipe(const psImage *image, // Image for which to generate a bg model + const pmConfig *config, // Configuration + psString recipe_name, + psString Xbin_name, + psString Ybin_name + ); + + /// Restore the background to an image bool ppBackgroundRestore( Index: /branches/eam_branches/ipp-20121130/ppBackground/src/ppBackgroundRestore.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppBackground/src/ppBackgroundRestore.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppBackground/src/ppBackgroundRestore.c (revision 34838) @@ -5,4 +5,32 @@ #include "ppBackground.h" + +psImageBinning *ppBackgroundBinningByRecipe(const psImage *image, // Image for which to generate a bg model + const pmConfig *config, // Configuration + psString recipe_name, + psString Xbin_name, + psString Ybin_name + ) +{ + bool status = true; + + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, recipe_name); + assert (recipe); + + // I have the fine image size, I know the binning factor, determine the ruff image size + psImageBinning *binning = psImageBinningAlloc(); + binning->nXfine = image->numCols; + binning->nYfine = image->numRows; + binning->nXbin = psMetadataLookupS32(&status, recipe, Xbin_name); + binning->nYbin = psMetadataLookupS32(&status, recipe, Ybin_name); + + psImageBinningSetRuffSize(binning, PS_IMAGE_BINNING_CENTER); + psImageBinningSetSkip(binning, image); + + return binning; +} + + + bool ppBackgroundRestore(pmChip *chip, const pmChip *background, const pmChip *pattern, @@ -27,5 +55,18 @@ if (background) { pmReadout *bgRO = pmFPAviewThisReadout(view, background->parent); // Readout with background - psImageBinning *binning = psphotBackgroundBinning(image, config); + + psImageBinning *binning; + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPBACKGROUND_RECIPE); // Recipe + if (!recipe) { + binning = psphotBackgroundBinning(image, config); + } + else { + binning = ppBackgroundBinningByRecipe(image,config, + psMetadataLookupStr(NULL,recipe,"BINNING_RECIPE"), + psMetadataLookupStr(NULL,recipe,"BINNING_XNAME"), + psMetadataLookupStr(NULL,recipe,"BINNING_YNAME")); + } + + fprintf(stderr,"%d %d %d %d\n",binning->nXfine,binning->nYfine,binning->nXbin,binning->nYbin); if (!binning) { psError(psErrorCodeLast(), false, "Unable to find background binning"); @@ -49,5 +90,9 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { + // if ((y % 250 == 0)&&(x % 250 == 0)) { + // printf("%d %d %g\n",x,y,bgImage->data.F32[y][x]); + // } image->data.F32[y][x] += bgImage->data.F32[y][x]; + } } Index: /branches/eam_branches/ipp-20121130/ppSkycell/src/ppSkycellCamera.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppSkycell/src/ppSkycellCamera.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppSkycell/src/ppSkycellCamera.c (revision 34838) @@ -103,11 +103,27 @@ return false; } - if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN1")) { + pmFPAfile *bin1 = pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN1"); + if (!bin1) { psError(psErrorCodeLast(), false, "Unable to define output."); return false; } - if (!pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN2")) { + pmFPAfile *bin2 = pmFPAfileDefineOutput(data->config, NULL, "PPSKYCELL.BIN2"); + if (!bin2) { psError(psErrorCodeLast(), false, "Unable to define output."); return false; + } + if (data->masksName) { + pmFPAfile *mask1 = pmFPAfileDefineOutput(data->config, bin1->fpa , "PPSKYCELL.BIN1.MASK"); + if (!mask1) { + psError(psErrorCodeLast(), false, "Unable to define output."); + return false; + } + mask1->save = true; + pmFPAfile *mask2 = pmFPAfileDefineOutput(data->config, bin2->fpa , "PPSKYCELL.BIN2.MASK"); + if (!mask2) { + psError(psErrorCodeLast(), false, "Unable to define output."); + return false; + } + mask2->save = true; } Index: /branches/eam_branches/ipp-20121130/ppSkycell/src/ppSkycellLoop.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppSkycell/src/ppSkycellLoop.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppSkycell/src/ppSkycellLoop.c (revision 34838) @@ -15,6 +15,6 @@ ) { - base->x0 = PS_MIN(base->x0, compare->x0); - base->x1 = PS_MAX(base->x1, compare->x1); + base->x0 = PS_MAX(base->x0, compare->x0); + base->x1 = PS_MIN(base->x1, compare->x1); base->y0 = PS_MIN(base->y0, compare->y0); base->y1 = PS_MAX(base->y1, compare->y1); @@ -27,5 +27,5 @@ { psPlane *fromCoords = psPlaneAlloc(), *toCoords = psPlaneAlloc(); // Coordinates for transforms - psRegion *region = psRegionAlloc(INFINITY, -INFINITY, INFINITY, -INFINITY); // Region for skycell + psRegion *region = psRegionAlloc(-INFINITY, INFINITY, INFINITY, -INFINITY); // Region for skycell // Limit the region using the nominated point (X,Y) @@ -38,6 +38,6 @@ toCoords->x += wcs->crpix1; \ toCoords->y += wcs->crpix2; \ - region->x0 = PS_MIN(region->x0, toCoords->x); \ - region->x1 = PS_MAX(region->x1, toCoords->x); \ + region->x0 = PS_MAX(region->x0, toCoords->x); \ + region->x1 = PS_MIN(region->x1, toCoords->x); \ region->y0 = PS_MIN(region->y0, toCoords->y); \ region->y1 = PS_MAX(region->y1, toCoords->y); @@ -132,5 +132,6 @@ psVector *target = psVectorAlloc(data->numInputs, PS_TYPE_S32); // Target for each input psArray *imageRegions = psArrayAlloc(data->numInputs); // Region for image - + psArray *regionHDUs = psArrayAlloc(data->numInputs); + // Determine which projection cells we have to deal with. for (int i = 0; i < data->numInputs; i++) { pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", i); // File to examine @@ -178,5 +179,5 @@ target->data.S32[i] = numProj; psTrace("ppSkycell", 3, "Image %d uses new projection\n", i); - + psArrayAdd(regionHDUs,1,hdu); numProj++; } @@ -184,10 +185,10 @@ pmFPAfileActivate(data->config->files, false, NULL); - + // Loop over projections for (int i = 0; i < numProj; i++) { psRegion *projRegion = projRegions->data[i]; // Region for skycell projection psTrace("ppSkycell", 2, "Projection %d: [%.0f:%.0f,%.0f:%.0f]\n", i, projRegion->x0, projRegion->x1, projRegion->y0, projRegion->y1); - int xSize = projRegion->x1 - projRegion->x0 + 1; // Size of unbinned image + int xSize = -projRegion->x1 + projRegion->x0 + 1; // Size of unbinned image int ySize = projRegion->y1 - projRegion->y0 + 1; // Size of unbinned image // Size of binned image 1 @@ -198,7 +199,7 @@ psImage *image1 = psImageAlloc(numCols1, numRows1, PS_TYPE_F32); // Binned image psImage *image2 = psImageAlloc(numCols2, numRows2, PS_TYPE_F32); // Binned image - psImageInit(image1, 0); - psImageInit(image2, 0); - + psImageInit(image1,NAN); + psImageInit(image2,NAN); + psImage *mask1 = NULL, *mask2 = NULL; // Binned masks if (data->masksName) { @@ -208,10 +209,12 @@ psImageInit(mask2, 0xFF); } - + pmHDU *projhdu = NULL; + int modify_wcs1 = 1; + int modify_wcs2 = 1; + // Loop over inputs to this projection. for (int j = 0; j < data->numInputs; j++) { if (target->data.S32[j] != i) { continue; } - pmFPAfileActivateSingle(data->config->files, true, "PPSKYCELL.IMAGE", j); if (data->masksName) { @@ -227,19 +230,24 @@ pmFPAfile *file = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.IMAGE", j); + if (!projhdu) { + projhdu = file->fpa->hdu; + // psMetadataCopy(projhdu->header,file->fpa->hdu->header); + } +#if 0 + else { + if ((psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX1") >= + psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))&& + (psMetadataLookupF32(NULL,file->fpa->hdu->header,"CRPIX2") >= + psMetadataLookupF32(NULL,projhdu->header,"CRPIX1"))) { + projhdu = file->fpa->hdu; + psMetadataCopy(projhdu->header,file->fpa->hdu->header); + } + } +#endif + pmReadout *inRO = pmFPAviewThisReadout(view, file->fpa); // Readout with input psFree(view); - // Flip images; no idea why this has to be done, but apparently it does - { - psImage *rot = psImageRotate(NULL, inRO->image, M_PI, NAN, PS_INTERPOLATE_BILINEAR); - psFree(inRO->image); - inRO->image = rot; - } - if (inRO->mask) { - psImage *rot = psImageRotate(NULL, inRO->mask, M_PI, 0, PS_INTERPOLATE_FLAT); - psFree(inRO->mask); - inRO->mask = rot; - } - + // data->maskVal = 0xffff; pmReadout *bin1RO = pmReadoutAlloc(NULL), *bin2RO = pmReadoutAlloc(NULL); // Binned readouts if (!pmReadoutRebin(bin1RO, inRO, data->maskVal, data->bin1, data->bin1)) { @@ -256,15 +264,14 @@ psRegion *imageRegion = imageRegions->data[j]; // Region for image // Offsets for image on skycell - int xOffset1 = (imageRegion->x0 - projRegion->x0) / (float)data->bin1; - int yOffset1 = (imageRegion->y0 - projRegion->y0) / (float)data->bin1; + int xOffset1 = (-imageRegion->x0 + projRegion->x0) / (float)data->bin1; + int yOffset1 = (-imageRegion->y1 + projRegion->y1) / (float)data->bin1; int xOffset2 = xOffset1 / (float)data->bin2, yOffset2 = yOffset1 / (float)data->bin2; - // XXX Completely neglecting rotations // The skycells are divided up neatly with them all having the same orientation - psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "="); - psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "="); + psImageOverlaySection(image1, bin1RO->image, xOffset1, yOffset1, "E"); + psImageOverlaySection(image2, bin2RO->image, xOffset2, yOffset2, "E"); if (data->masksName) { - psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "="); - psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "="); + psImageOverlaySection(mask1, bin1RO->mask, xOffset1, yOffset1, "M"); + psImageOverlaySection(mask2, bin2RO->mask, xOffset2, yOffset2, "M"); } @@ -287,4 +294,8 @@ pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1"); pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2"); + if (data->masksName) { + pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN1.MASK"); + pmFPAfileActivate(data->config->files, true, "PPSKYCELL.BIN2.MASK"); + } } @@ -313,36 +324,161 @@ if (data->doFits) { + + pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0); + // Copy header from projection root hdu + fits1->fpa->hdu = pmHDUAlloc(NULL); + fits1->fpa->hdu->header = psMetadataAlloc(); + psMetadataCopy(fits1->fpa->hdu->header,projhdu->header); + + // Change wcs here +#define WCS_DEBUG 0 + if (modify_wcs1) { + pmAstromWCS *WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); + double cd1f = 1.0 * data->bin1; + double cd2f = 1.0 * data->bin1; +#if WCS_DEBUG + fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2, + cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, + WCS->crpix1,WCS->crpix2 + ); +#endif + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; +#if WCS_DEBUG + fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g) %d %d %d %d\n",data->bin1,data->bin2, + cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, + WCS->crpix1,WCS->crpix2, + WCS->trans->x->nX, WCS->trans->x->nY, + WCS->trans->y->nX, WCS->trans->y->nY + ); +#endif + for (int q = 0; q <= WCS->trans->x->nX; q++) { + for (int r = 0; r <= WCS->trans->x->nY; r++) { + WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); +#if WCS_DEBUG + fprintf(stderr,"PC: %d %d %g %g\n", + q,r,pow(cd1f,q) * pow(cd2f,r), + WCS->trans->x->coeff[q][r]); +#endif + } + } + for (int q = 0; q <= WCS->trans->y->nX; q++) { + for (int r = 0; r <= WCS->trans->y->nY; r++) { + WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); +#if WCS_DEBUG + fprintf(stderr,"PC: %d %d %g %g\n", + q,r,pow(cd1f,q) * pow(cd2f,r), + WCS->trans->y->coeff[q][r]); +#endif + } + } + pmAstromWCStoHeader (fits1->fpa->hdu->header,WCS); + WCS = pmAstromWCSfromHeader(fits1->fpa->hdu->header); +#if WCS_DEBUG + fprintf(stderr,">>> %d %d (%g %g) (%g %g) (%g %g)\n",data->bin1,data->bin2, + cd1f,cd2f,WCS->cdelt1,WCS->cdelt2, + WCS->crpix1,WCS->crpix2 + ); +#endif + modify_wcs1 = 0; + } + + + pmChip *Fchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1"); + psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); + psMetadataAddS32(Fchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); + pmCell *Fcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1"); // Rebinned cell 1 - pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2 - - // This is a hack to get a functioning header created so the fits images can be written out. +/* // This is a hack to get a functioning header created so the fits images can be written out. */ psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); + + pmReadout *Fro1 = pmReadoutAlloc(Fcell1); + Fro1->image = image1; + Fro1->mask = mask1; + Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; + + fits1->save = true; + fits1->fileIndex = i; + + // Do the mask if we need to + if (data->masksName) { +/* pmFPAfile *Mask1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1.MASK", 0); */ +/* // Mask1->fpa->hdu = pmHDUAlloc(NULL); */ +/* // Mask1->fpa->hdu->header = psMetadataAlloc(); */ +/* // psMetadataCopy(Mask1->fpa->hdu->header,fits1->fpa->hdu->header); */ +/* pmChip *Mchip1 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN1.MASK"); */ +/* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); */ +/* psMetadataAddS32(Mchip1->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); */ + +/* pmCell *Mcell1 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN1.MASK"); // Rebinned cell 1 */ +/* /\* // This is a hack to get a functioning header created so the fits images can be written out. *\/ */ +/* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); */ +/* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); */ +/* psMetadataAddS32(Mcell1->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); */ + +/* pmReadout *Mro1 = pmReadoutAlloc(Mcell1); */ +/* Mro1->image = image1; */ +/* Mro1->mask = mask1; */ +/* Mro1->data_exists = Mcell1->data_exists = Mcell1->parent->data_exists = true; */ + +/* Mask1->save = true; */ +/* Mask1->fileIndex = i; */ + } + + + + // Repeat with second binned image + pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0); + fits2->fpa->hdu = pmHDUAlloc(NULL); + fits2->fpa->hdu->header = psMetadataAlloc(); + psMetadataCopy(fits2->fpa->hdu->header,projhdu->header); + + if (modify_wcs2) { + pmAstromWCS *WCS = pmAstromWCSfromHeader(fits2->fpa->hdu->header); + double cd1f = 1.0 * data->bin2 * data->bin1; + double cd2f = 1.0 * data->bin2 * data->bin1; + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + + for (int q = 0; q < WCS->trans->x->nX; q++) { + for (int r = 0; r < WCS->trans->x->nY; r++) { + WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + } + } + for (int q = 0; q < WCS->trans->y->nX; q++) { + for (int r = 0; r < WCS->trans->y->nY; r++) { + WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + } + } + pmAstromWCStoHeader (fits2->fpa->hdu->header,WCS); + modify_wcs2 = 0; + } + + pmChip *Fchip2 = pmFPAfileThisChip(data->config->files, view, "PPSKYCELL.BIN2"); + psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); + psMetadataAddS32(Fchip2->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); + + + pmCell *Fcell2 = pmFPAfileThisCell(data->config->files, view, "PPSKYCELL.BIN2"); // Rebinned cell 2 psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.XPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.YPARITY", PS_META_REPLACE,"",1); psMetadataAddS32(Fcell2->concepts,PS_LIST_TAIL,"CELL.READDIR", PS_META_REPLACE,"",1); - - psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); - psMetadataAddS32(Fcell1->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); - psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.XPARITY", PS_META_REPLACE,"",1); - psMetadataAddS32(Fcell2->parent->concepts,PS_LIST_TAIL,"CHIP.YPARITY", PS_META_REPLACE,"",1); - - pmReadout *Fro1 = pmReadoutAlloc(Fcell1), *Fro2 = pmReadoutAlloc(Fcell2); // Binned readouts - - Fro1->image = image1; + + pmReadout *Fro2 = pmReadoutAlloc(Fcell2); Fro2->image = image2; - Fro1->mask = mask1; Fro2->mask = mask2; - - Fro1->data_exists = Fcell1->data_exists = Fcell1->parent->data_exists = true; Fro2->data_exists = Fcell2->data_exists = Fcell2->parent->data_exists = true; - - pmFPAfile *fits1 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN1", 0); - fits1->save = true; - fits1->fileIndex = i; - pmFPAfile *fits2 = pmFPAfileSelectSingle(data->config->files, "PPSKYCELL.BIN2", 0); + fits2->save = true; fits2->fileIndex = i; + + } Index: /branches/eam_branches/ipp-20121130/ppStack/src/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/Makefile.am (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/Makefile.am (revision 34838) @@ -1,3 +1,3 @@ -bin_PROGRAMS = ppStack +bin_PROGRAMS = ppStack ppStackMedian if HAVE_SVNVERSION @@ -25,4 +25,7 @@ ppStack_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_LIBS) $(PPSTATS_LIBS) $(PPSTACK_LIBS) +ppStackMedian_CFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS) +ppStackMedian_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_LIBS) $(PPSTATS_LIBS) $(PPSTACK_LIBS) + ppStack_SOURCES = \ ppStack.c \ @@ -43,4 +46,35 @@ ppStackCombinePrepare.c \ ppStackCombineInitial.c \ + ppStackCombineAlternate.c \ + ppStackReject.c \ + ppStackCombineFinal.c \ + ppStackCleanup.c \ + ppStackPhotometry.c \ + ppStackFinish.c \ + ppStackTarget.c \ + ppStackUpdateHeader.c \ + ppStackJPEGs.c \ + ppStackStats.c \ + ppStackErrorCodes.c + +ppStackMedian_SOURCES = \ + ppStackMedian.c \ + ppStackArguments.c \ + ppStackCamera.c \ + ppStackFiles.c \ + ppStackMedianLoop.c \ + ppStackPSF.c \ + ppStackReadout.c \ + ppStackVersion.c \ + ppStackMatch.c \ + ppStackSources.c \ + ppStackThread.c \ + ppStackOptions.c \ + ppStackSetup.c \ + ppStackPrepare.c \ + ppStackConvolve.c \ + ppStackCombinePrepare.c \ + ppStackCombineInitial.c \ + ppStackCombineAlternate.c \ ppStackReject.c \ ppStackCombineFinal.c \ Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStack.h =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStack.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStack.h (revision 34838) @@ -38,5 +38,8 @@ PPSTACK_FILES_STACK, // Stack files PPSTACK_FILES_UNCONV, // Unconvolved stack files - PPSTACK_FILES_PHOT // Files for photometry + PPSTACK_FILES_PHOT, // Files for photometry + PPSTACK_FILES_BKG, // Files for bkg + PPSTACK_FILES_MEDIAN_IN, // Files for median only stacks. + PPSTACK_FILES_MEDIAN_OUT // Files for median only stacks. } ppStackFileList; @@ -111,4 +114,7 @@ bool ppStackReadoutFinalThread(psThreadJob *job // Job to process ); +// Perform median stacking for background +bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); +bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); // Return software version Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackArguments.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackArguments.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackArguments.c (revision 34838) @@ -324,4 +324,5 @@ } + psTrace("ppStack", 1, "Done parsing arguments\n"); Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCamera.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCamera.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCamera.c (revision 34838) @@ -156,5 +156,6 @@ psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources - + psString bkgmodel = psMetadataLookupStr(&mdok, input, "BKGMODEL"); // Name of warped background model + pmFPAfile *imageFile = defineFile(config, NULL, "PPSTACK.INPUT", image, PM_FPA_FILE_IMAGE); // File for image @@ -215,4 +216,21 @@ } + // Grab bkgmodel information here + if (!bkgmodel || strlen(bkgmodel) == 0) { + // We have no background models. + } + else { + pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel, + PM_FPA_FILE_IMAGE); + if (!inputBKG) { + psError(psErrorCodeLast(), false, + "Unable to define file from bkgmodel %d (%s)",i,bkgmodel); + return(false); + } + }// End bkgmodel + + + + i++; } @@ -460,4 +478,25 @@ jpeg2->save = true; + // Output background + pmFPAfile *outBkg = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGMODEL"); + if (!outBkg) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGMODEL")); + return(false); + } + if (outBkg->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGMODEL is not of type IMAGE"); + return(false); + } + outBkg->save = true; + pmFPAfile *outBkgRest = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGREST"); + if (!outBkgRest) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGREST")); + return(false); + } + if (outBkgRest->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGREST is not of type IMAGE"); + return(false); + } + // For photometry, we operate on the chip-mosaicked image // we create a copy of the mosaicked image for psphot so we can write out a clean image @@ -496,4 +535,6 @@ } + // Define output file here. + return true; } Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCleanup.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCleanup.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCleanup.c (revision 34838) @@ -22,10 +22,24 @@ options->outRO = NULL; - options->expRO->data_exists = false; - options->expRO->parent->data_exists = false; - options->expRO->parent->parent->data_exists = false; - psFree(options->expRO); - options->expRO = NULL; + if (options->expRO) { + options->expRO->data_exists = false; + if (options->expRO->parent) { + options->expRO->parent->data_exists = false; + options->expRO->parent->parent->data_exists = false; + } + psFree(options->expRO); + options->expRO = NULL; + } + if (options->bkgRO) { + options->bkgRO->data_exists = false; + if (options->bkgRO->parent) { + options->bkgRO->parent->data_exists = false; + options->bkgRO->parent->parent->data_exists = false; + } + psFree(options->bkgRO); + options->bkgRO = NULL; + } + for (int i = 0; i < options->num; i++) { pmCellFreeData(options->cells->data[i]); @@ -62,7 +76,6 @@ bool ppStackCleanup (pmConfig *config, ppStackOptions *options) { - psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code - + // Ensure everything closes if (config) { @@ -72,4 +85,6 @@ ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true); ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); + ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_IN, true); + ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_OUT, true); if (!ppStackFilesIterateUp(config)) { psError(psErrorCodeLast(), false, "Unable to close files."); Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCombineAlternate.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCombineAlternate.c (revision 34838) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCombineAlternate.c (revision 34838) @@ -0,0 +1,115 @@ +#include "ppStack.h" + +// This is the doomsday switch. +// #define TESTING // Enable test output +bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_BKGMED"); + + + pmReadout *outRO = options->outRO; + + psArray *inputs = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + inputs->data[i] = ro; + pmFPAfileClose(file,view); + } + if (!pmStackSimpleMedianCombine(outRO,inputs)) { + psFree(inputs); + return(false); + } +#if 0 + if (!ppStackWriteImage("/tmp/test_forced.median.fits", + outRO->parent->parent->parent->hdu->header, + outRO->image, + config)) { + fprintf(stderr,"Failed to write image because fail.\n"); + } +#endif + for (int i = 0; i < options->num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + pmFPAview *view = ppStackFilesIterateDown(config); + bool success = pmFPAfileClose(file,view); + if (!success) { + psTrace("ppStack",5,"I failed at closing a file.\n"); + } + + psFitsClose(file->fits); + file->fits = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + file->wrote_phu = false; + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, false, i); + } + psFree(inputs); + outRO->data_exists = true; + outRO->parent->data_exists = true; + outRO->parent->parent->data_exists = true; + + return(true); +} + + +bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_BKGMED"); + + + pmReadout *bkgRO = options->bkgRO; + + psArray *inputs = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + inputs->data[i] = ro; + pmFPAfileClose(file,view); + } + if (!pmStackSimpleMedianCombine(bkgRO,inputs)) { + psFree(inputs); + return(false); + } +#if 0 + if (!ppStackWriteImage("/tmp/test_forced.bkgmdl.fits", + bkgRO->parent->parent->parent->hdu->header, + bkgRO->image, + config)) { + fprintf(stderr,"Failed to write image because fail.\n"); + } +#endif + for (int i = 0; i < options->num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); + pmFPAview *view = ppStackFilesIterateDown(config); + bool success = pmFPAfileClose(file,view); + if (!success) { + psTrace("ppStack",5,"I failed at closing a file.\n"); + } + + psFitsClose(file->fits); + file->fits = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + file->wrote_phu = false; + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, false, i); + } + psFree(inputs); + bkgRO->data_exists = true; + bkgRO->parent->data_exists = true; + bkgRO->parent->parent->data_exists = true; + + return(true); +} + Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCombinePrepare.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCombinePrepare.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackCombinePrepare.c (revision 34838) @@ -1,5 +1,5 @@ #include "ppStack.h" -bool ppStackCombinePrepare(const char *outName, const char *expName, +bool ppStackCombinePrepare(const char *outName, const char *expName, const char *bkgName, ppStackFileList files, ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) @@ -27,6 +27,46 @@ options->outRO = pmReadoutAlloc(cell); // Output readout - pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell - options->expRO = pmReadoutAlloc(expCell); // Output readout + if (expName) { + pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell + options->expRO = pmReadoutAlloc(expCell); //Output readout + } +/* else { */ +/* options->expRO = NULL; */ +/* } */ + pmCell *bkgCell; + int bkg_r0,bkg_c0; + int bkg_nC,bkg_nR; + if (bkgName) { + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, 0); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", 0); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + + bkg_r0 = ro->image->row0; + bkg_c0 = ro->image->col0; + bkg_nC = ro->image->numCols; + bkg_nR = ro->image->numRows; + bkgCell = pmFPAfileThisCell(config->files, view, bkgName); // Bkg cell + + options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout + // if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) { + options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL); + if (!options->bkgRO->parent->parent->parent->hdu) { + fprintf(stderr,"failed to generate a HDU for this thing.\n"); + } + options->bkgRO->parent->parent->parent->hdu->header = psMetadataCopy(options->bkgRO->parent->parent->parent->hdu->header, + ro->parent->parent->parent->hdu->header); + + options->bkgRO->parent->concepts = psMetadataCopy(options->bkgRO->parent->concepts, + ro->parent->concepts); + options->bkgRO->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->concepts, + ro->parent->parent->concepts); + options->bkgRO->parent->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->parent->concepts, + ro->parent->parent->parent->concepts); + + } + else { + options->bkgRO = NULL; + } psFree(view); @@ -42,7 +82,19 @@ } - if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) { + if (expName) { + if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) { psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); return false; + } + } + + if (bkgName) { + if (!pmReadoutStackDefineOutput(options->bkgRO, bkg_c0, bkg_r0, bkg_nC, bkg_nR, false, false, 0)) { + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); + return false; + } + + + } Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackFiles.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackFiles.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackFiles.c (revision 34838) @@ -14,4 +14,14 @@ /// Files required for the convolution static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL }; + +/// Files required for the background +static char *filesBkg[] = { "PPSTACK.INPUT.BKGMODEL", NULL }; + +/// Files required for median only stacking +static char *filesMedianIn[] = { "PPSTACK.INPUT", + NULL }; +/// Files required for median only stacking +static char *filesMedianOut[] = { "PPSTACK.OUTPUT", + NULL }; /// Regular (convolved) stack files @@ -19,4 +29,5 @@ "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT", "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", + "PPSTACK.OUTPUT.BKGMODEL", NULL }; /// Unconvolved stack files @@ -41,4 +52,7 @@ case PPSTACK_FILES_UNCONV: return filesUnconv; case PPSTACK_FILES_PHOT: return filesPhot; + case PPSTACK_FILES_BKG: return filesBkg; + case PPSTACK_FILES_MEDIAN_IN: return filesMedianIn; + case PPSTACK_FILES_MEDIAN_OUT: return filesMedianOut; default: psAbort("Unrecognised file list: %x", list); Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackLoop.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackLoop.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackLoop.c (revision 34838) @@ -63,5 +63,5 @@ // Prepare for combination - if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK, stack, options, config)) { + if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.BKGMODEL", PPSTACK_FILES_STACK, stack, options, config)) { psError(psErrorCodeLast(), false, "Unable to prepare for combination."); psFree(stack); @@ -84,5 +84,5 @@ pmCellFreeData(options->cells->data[i]); } - psFree(stack); + // psFree(stack); // Pixel rejection @@ -142,4 +142,11 @@ } + // Generate median background stack here. + if (!ppStackCombineBackground(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to generate median of background images."); + psFree(stack); + return false; + } + ppStackFileActivation(config,PPSTACK_FILES_BKG ,false); // Photometry psTrace("ppStack", 1, "Photometering stacked image....\n"); @@ -171,5 +178,5 @@ return false; } - psFree(stack); + // psFree(stack); psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("cleanup"); @@ -186,5 +193,5 @@ // Prepare for combination - if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV, + if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV, stack, options, config)) { psError(psErrorCodeLast(), false, "Unable to prepare for combination."); Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackLoop.h =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackLoop.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackLoop.h (revision 34838) @@ -4,4 +4,9 @@ // Loop over the inputs, doing the combination bool ppStackLoop( + pmConfig *config, // Configuration + ppStackOptions *options // Options for stacking + ); +// Median only loop. +bool ppStackMedianLoop( pmConfig *config, // Configuration ppStackOptions *options // Options for stacking @@ -32,4 +37,5 @@ const char *outName, // Name of output file const char *expName, // Name of exposure file + const char *bkgName, // Name of background file ppStackFileList files, // Files of interest ppStackThreadData *stack, // Stack Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackMedian.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackMedian.c (revision 34838) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackMedian.c (revision 34838) @@ -0,0 +1,51 @@ +#include "ppStack.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + psTimerStart("PPSTACK"); + psTimerStart("PPSTACK_STEPS"); + + pmErrorRegister(); + ppStackErrorRegister(); + psphotErrorRegister(); + + ppStackOptions *options = NULL; // Options for stacking + + pmConfig *config = pmConfigRead(&argc, argv, PPSTACK_RECIPE); // Configuration + if (!config) { + ppStackCleanup(config, options); + } + + ppStackVersionPrint(); + + if (!pmModelClassInit()) { + psError(PPSTACK_ERR_PROG, false, "Unable to initialise model classes."); + ppStackCleanup(config, options); + } + + if (!psphotInit()) { + psError(PPSTACK_ERR_PROG, false, "Error initialising psphot."); + ppStackCleanup(config, options); + } + + if (!ppStackArgumentsSetup(argc, argv, config)) { + ppStackCleanup(config, options); + } + + if (!ppStackCamera(config)) { + ppStackCleanup(config, options); + } + + if (!ppStackArgumentsParse(config)) { + ppStackCleanup(config, options); + } + + options = ppStackOptionsAlloc(); + if (!ppStackMedianLoop(config, options)) { + ppStackCleanup(config, options); + } + + ppStackCleanup(config, options); +} + Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackMedianLoop.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackMedianLoop.c (revision 34838) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackMedianLoop.c (revision 34838) @@ -0,0 +1,204 @@ +#include "ppStack.h" + +// static functions are defined below +/* static int stackSummary(const ppStackOptions *options, const char *place); */ + +bool ppStackMedianLoop(pmConfig *config, ppStackOptions *options) +{ + assert(config); + + psTimerStart("PPSTACK_TOTAL"); + psTimerStart("PPSTACK_STEPS"); + + // Setup + psTrace("ppStack", 1, "Setup....\n"); + if (!ppStackSetup(options, config)) { + psError(psErrorCodeLast(), false, "Unable to setup."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("setup"); + + // Duplicate code that is in ppStackConvolve + options->cells = psArrayAlloc(options->num); // Cells for convolved images --- a handle for reading again + pmFPAview *view = ppStackFilesIterateDown(config); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest + + // pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + if (!view) { + return false; + } + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout + pmCell *inCell = readout->parent; // Input cell */ + options->cells->data[i] = psMemIncrRefCounter(inCell); + + } + psFree(view); + + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, options->convolve); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Preparation for stacking + if (!ppStackCombinePrepare("PPSTACK.OUTPUT",NULL,NULL,PPSTACK_FILES_MEDIAN_OUT,stack,options,config)) { + psError(psErrorCodeLast(), false, "Unabel to combine data."); + psFree(stack); + return false; + } + + if (!ppStackCombineMedian(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unabel to combine data."); + psFree(stack); + return false; + } +/* for (int x = 0; x < options->outRO->image->numCols; x++) { */ +/* for (int y = 0; y < options->outRO->image->numRows; y++) { */ +/* options->outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; */ +/* } */ +/* } */ + psFree(options->outRO->mask); + // Do this before performing photometry so that the cmf header gets all of the information. + options->zpInput = NULL; + options->expTimeInput = NULL; + options->airmassInput = NULL; + if (!ppStackUpdateHeader(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to update header."); + psFree(stack); + return false; + } +/* // Generate median background stack here. */ +/* if (!ppStackCombineBackground(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to generate median of background images."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +// ppStackFileActivation(config,PPSTACK_FILES_BKG ,false); +/* // Photometry */ +/* psTrace("ppStack", 1, "Photometering stacked image....\n"); */ +/* if (!ppStackPhotometry(options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to perform photometry."); */ +/* return false; */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); */ +/* ppStackMemDump("photometry"); */ + + // Create JPEGS +/* if (!ppStackJPEGs(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to make jpegs."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* // Assemble Stats */ +/* if (!ppStackStats(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to assemble statistics."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ + + // Clean up + + psTrace("ppStack", 2, "Cleaning up after combination....\n"); + if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_MEDIAN_OUT, PPSTACK_FILES_MEDIAN_OUT, false)) { + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + // psFree(stack); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("cleanup"); + +/* // Unconvolved stack --- it's cheap to calculate, compared to everything else! */ +/* // XXX unconvolved stack is currently using the convolved mask! oops! */ +/* if (options->convolve) { */ +/* // Start threading */ +/* ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false); */ +/* if (!stack) { */ +/* psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); */ +/* return false; */ +/* } */ + +/* // Prepare for combination */ +/* if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV, */ +/* stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to prepare for combination."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ + +/* // generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been */ +/* psTrace("ppStack", 2, "Stack of unconvolved images....\n"); */ +/* if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false)) { */ +/* psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); */ +/* ppStackMemDump("unconv"); */ + +/* // Update Header */ +/* if (!ppStackUpdateHeader(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to update header."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* // Clean up unconvolved stack */ +/* psTrace("ppStack", 2, "Cleaning up after unconvolved stack....\n"); */ +/* if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_UNCONV, PPSTACK_FILES_NONE, false)) { */ +/* psError(psErrorCodeLast(), false, "Unable to clean up."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* psFree(stack); */ +/* } */ + psFree(options->cells); options->cells = NULL; + + // Finish up + psTrace("ppStack", 1, "Finishing up....\n"); + if (!ppStackFinish(options, config)) { + psError(psErrorCodeLast(), false, "Unable to finish up."); + return false; + } + ppStackMemDump("finish"); + + return true; +} + +/* /// Print a summary of the inputs, and return the number of good inputs */ +/* static int stackSummary(const ppStackOptions *options, const char *place) */ +/* { */ +/* int numGood = 0; // Number of good inputs */ +/* psString summary = NULL; // Summary of images */ +/* for (int i = 0; i < options->num; i++) { */ +/* psString reason = NULL; // Reason for rejecting */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) { */ +/* psStringAppend(&reason, " Good."); */ +/* numGood++; */ +/* } else { */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CAL) { */ +/* psStringAppend(&reason, " Calibration failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_PSF) { */ +/* psStringAppend(&reason, " PSF measurement failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_MATCH) { */ +/* psStringAppend(&reason, " PSF matching failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CHI2) { */ +/* psStringAppend(&reason, " PSF matching chi^2 deviant."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_REJECT) { */ +/* psStringAppend(&reason, " Rejection exceeded threshold."); */ +/* } */ +/* } */ +/* psStringAppend(&summary, "Image %d: %s\n", i, reason); */ +/* psFree(reason); */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Summary of images for %s:\n%s", place, summary); */ +/* psFree(summary); */ + +/* return numGood; */ +/* } */ Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackOptions.h =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackOptions.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackOptions.h (revision 34838) @@ -48,4 +48,7 @@ // Rejection psArray *rejected; // Rejected pixels + // Background + pmReadout *bkgRO; // Output background readout + psArray *bkgImages; // Input background images } ppStackOptions; Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackPrepare.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackPrepare.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackPrepare.c (revision 34838) @@ -132,4 +132,5 @@ psVectorInit(options->inputMask, 0); options->exposures = psVectorAlloc(options->num, PS_TYPE_F32); + psVectorInit(options->exposures, NAN); @@ -150,4 +151,8 @@ options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); + if ((options->exposures->data.F32[i] == 0)|| + (!(isfinite(options->exposures->data.F32[i])))){ + options->exposures->data.F32[i] = psMetadataLookupF32(NULL,recipe,"DEFAULT.EXPTIME"); + } options->sumExposure += options->exposures->data.F32[i]; Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackReadout.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackReadout.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackReadout.c (revision 34838) @@ -194,5 +194,4 @@ - bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackSetup.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackSetup.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackSetup.c (revision 34838) @@ -89,5 +89,8 @@ options->origMasks = psArrayAlloc(num); options->origVariances = psArrayAlloc(num); + options->bkgImages = psArrayAlloc(num); pmFPAview *view = pmFPAviewAlloc(0); + int nullMasks = 0; + int nullVariances = 0; for (int i = 0; i < num; i++) { { @@ -97,11 +100,34 @@ { // We want the convolved mask, since that defines the area that has been tested for outliers + if (options->convolve) { options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]); + } + else { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.MASK", i); + options->origMasks->data[i] = pmFPAfileName(file, view, config); + } + if (!(options->origMasks->data[i])) { + nullMasks++; + } } { pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i); options->origVariances->data[i] = pmFPAfileName(file, view, config); + if (!(options->origVariances->data[i])) { + nullVariances++; + } } +/* { */ +/* pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); */ +/* options->bkgImages->data[i] = pmFPAfileName(file, view, config); */ +/* } */ } + if (nullMasks == num) { + psFree(options->origMasks); + } + if (nullVariances == num) { + psFree(options->origVariances); + } + psFree(view); Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackThread.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackThread.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackThread.c (revision 34838) @@ -40,4 +40,5 @@ psFree(stack->maskFits); psFree(stack->varianceFits); + psFree(stack->bkgFits); return; } @@ -67,5 +68,5 @@ PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL); } - + ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree); @@ -77,4 +78,5 @@ stack->maskFits = psArrayAlloc(numInputs); stack->varianceFits = psArrayAlloc(numInputs); + stack->bkgFits = psArrayAlloc(numInputs); for (int i = 0; i < numInputs; i++) { if (!cells->data[i]) { @@ -95,8 +97,12 @@ psFree(resolved); \ } - + IMAGE_OPEN(imageNames, stack->imageFits, i); - IMAGE_OPEN(maskNames, stack->maskFits, i); - IMAGE_OPEN(varianceNames, stack->varianceFits, i); + if (maskNames) { + IMAGE_OPEN(maskNames, stack->maskFits, i); + } + if (varianceNames) { + IMAGE_OPEN(varianceNames, stack->varianceFits, i); + } } Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackThread.h =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackThread.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackThread.h (revision 34838) @@ -25,4 +25,5 @@ psArray *maskFits; // FITS file pointers for masks psArray *varianceFits; // FITS file pointers for variances + psArray *bkgFits; // FITS file pointers for background models } ppStackThreadData; Index: /branches/eam_branches/ipp-20121130/ppStack/src/ppStackUpdateHeader.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppStack/src/ppStackUpdateHeader.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppStack/src/ppStackUpdateHeader.c (revision 34838) @@ -6,9 +6,10 @@ pmReadout *outRO = options->outRO; // Output readout + pmReadout *expRO = options->expRO; // Propagate WCS bool wcsDone = false; // Have we done the WCS? for (int i = 0; i < options->num && !wcsDone; i++) { - if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + if ((options->inputMask)&&(options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i])) { continue; } @@ -123,5 +124,23 @@ snprintf (field, 64, "AIR_%04d", i); psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value); - } + } + + // Copy information into expRO, because it should be there too. + if ((expRO)&&(expRO->parent)) { + pmHDU *expROhdu = pmHDUFromCell(expRO->parent); + if (!expROhdu->header) { + expROhdu->header = psMetadataAlloc(); + } + expRO->parent->parent->parent->hdu->header = psMetadataCopy(expRO->parent->parent->parent->hdu->header, + outRO->parent->parent->parent->hdu->header); + + expRO->parent->concepts = psMetadataCopy(expRO->parent->concepts, + outRO->parent->concepts); + expRO->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->concepts, + outRO->parent->concepts); + expRO->parent->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->parent->concepts, + outRO->parent->parent->parent->concepts); + } + return true; } Index: /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMops.h =================================================================== --- /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMops.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMops.h (revision 34838) @@ -62,5 +62,5 @@ } ppMopsDetections; -ppMopsDetections *ppMopsDetectionsAlloc(); +ppMopsDetections *ppMopsDetectionsAlloc(void); /// Copy a detection Index: /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsDetections.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsDetections.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsDetections.c (revision 34838) @@ -23,5 +23,5 @@ } -ppMopsDetections *ppMopsDetectionsAlloc() +ppMopsDetections *ppMopsDetectionsAlloc(void) { ppMopsDetections *det = psAlloc(sizeof(ppMopsDetections)); // Detections, to return Index: /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsRead.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsRead.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsRead.c (revision 34838) @@ -28,5 +28,5 @@ for (int i = 0; i < num; i++) { const char *name = inNames->data[i]; - printf("%s\n", name); + printf("Input filename:\n%s\n", name); psFits *fits = psFitsOpen(name, "r"); // FITS file @@ -83,5 +83,6 @@ continue; } - ppMopsDetections *det = ppMopsDetectionsAlloc(size); + ppMopsDetections *det = ppMopsDetectionsAlloc(); + det->platescale = NAN; detections->data[i] = det; det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension @@ -280,4 +281,5 @@ if (isfinite(scale)) { det->platescale = scale; + //printf("platescale = %g\n", scale); } numGood++; @@ -306,4 +308,5 @@ } + printf("Detection platescale = %g\n", det->platescale); psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name); } Index: /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsWrite.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsWrite.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppTranslate/src/ppMopsWrite.c (revision 34838) @@ -51,4 +51,5 @@ float seeing = 0.; int totalGood = 0; + int platescale_has_been_set = 0; for (int d = 0; d < detections->n; d++) { if ( (detections->data[d] != NULL) && (isfinite(((ppMopsDetections*) detections->data[d])->seeing)) ) { @@ -56,9 +57,19 @@ totalGood += 1; if (isfinite(((ppMopsDetections*) detections->data[d])->platescale)) { + if (platescale_has_been_set == 0) { det->platescale = ((ppMopsDetections*) detections->data[d])->platescale; - } - } - } - printf("plate-scale = [%g]\n", det->platescale); + platescale_has_been_set = 1; + } else { + if (det->platescale != ((ppMopsDetections*) detections->data[d])->platescale) { + printf("Different values for platescale: %g - %g in detection %d\n", + det->platescale, + ((ppMopsDetections*) detections->data[d])->platescale, + d); + } + } + } + } + } + //printf("plate-scale = [%g]\n", det->platescale); seeing /= (float) totalGood; if (det != NULL) { Index: /branches/eam_branches/ipp-20121130/ppViz/src/ppVizPSF/ppVizPSFCamera.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppViz/src/ppVizPSF/ppVizPSFCamera.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppViz/src/ppVizPSF/ppVizPSFCamera.c (revision 34838) @@ -19,4 +19,5 @@ files->data[0] = psStringCopy(name); if (psMetadataLookup(config->arguments, file)) { + psMetadataRemoveKey(config->arguments, file); } @@ -42,5 +43,6 @@ fileArguments("SOURCES", data->sourcesName, "Input sources", data->config); pmFPAfile *srcs = pmFPAfileBindFromArgs(&status, psf, data->config, - "PSPHOT.INPUT.CMF", "SOURCES"); // File + "PSPHOT.INPUT.CMF", "SOURCES"); // File + fprintf(stderr,"%ld %d\n",(long) srcs, status); if (!status || !srcs) { psError(PS_ERR_IO, false, "Failed to build file from PSPHOT.INPUT.CMF"); Index: /branches/eam_branches/ipp-20121130/ppViz/src/ppVizPSF/ppVizPSFLoop.c =================================================================== --- /branches/eam_branches/ipp-20121130/ppViz/src/ppVizPSF/ppVizPSFLoop.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/ppViz/src/ppVizPSF/ppVizPSFLoop.c (revision 34838) @@ -20,5 +20,5 @@ return NULL; } - + fprintf(stderr,"Woo!\n"); pmChip *chip; // Chip from FPA while ((chip = pmFPAviewNextChip(view, psfFile->fpa, 1))) { @@ -56,15 +56,13 @@ psWarning("More than one readout present for chip %d, cell %d", view->chip, view->cell); } - + fprintf(stderr,"Woo!\n"); pmReadout *readout; // Readout from cell while ((readout = pmFPAviewNextReadout(view, psfFile->fpa, 1))) { + fprintf(stderr,"Woo?\n"); if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { psError(PS_ERR_UNKNOWN, false, "Error loading data from files."); return false; } - if (!readout->data_exists) { - continue; - } - + fprintf(stderr,"Woo2?\n"); pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF assert(psf); @@ -72,4 +70,12 @@ bool mdok; // Status of MD lookup psArray *sources = psMetadataLookupPtr(&mdok, readout->analysis, "PSPHOT.SOURCES"); // Sources + + if (sources && !readout->data_exists) { // This fails if -sources not specified + continue; + } + readout->data_exists = true; + fprintf(stderr,"Woo! %ld\n", sources ? sources->n : -1); + + fprintf(stderr,"%d\n",mdok); int numCols = 0, numRows = 0; // Size of image psVector *xOffset = NULL, *yOffset = NULL; // Offset from source to true position @@ -80,5 +86,5 @@ psLogMsg("ppVizPSF", PS_LOG_INFO, "Generating %dx%d image", numCols, numRows); } - if (sources) { + if (sources && !sources->n) { psMemIncrRefCounter(sources); psLogMsg("ppVizPSF", PS_LOG_INFO, "Using %ld input sources from CMF file", sources->n); @@ -89,4 +95,5 @@ if (data->fakeNum > 0 && isfinite(data->fakeMag)) { + fprintf(stderr,"Here! fakes\n"); long numOld = 0; // Old number of sources long numNew = -1; // New number of sources @@ -118,4 +125,5 @@ } if (!sources && !data->input) { + fprintf(stderr,"Here! default\n"); // Generate fake image with only a single realisation of the PSF sources = psArrayAlloc(1); @@ -151,6 +159,7 @@ psFree(psf->residuals); psf->residuals = NULL; - + fprintf(stderr,"still here"); if (sources) { + fprintf(stderr,"Here! sources?\n"); if (!pmReadoutFakeFromSources(readout, numCols, numRows, sources, 0, xOffset, yOffset, psf, data->minFlux, 0, false, true)) { @@ -159,7 +168,9 @@ } } else if (data->input) { + psVector *x = data->input->data[0]; // x coordinates psVector *y = data->input->data[1]; // y coordinates psVector *mag = data->input->data[2]; // Magnitudes + fprintf(stderr,"Here! input? %d %d %g %d %ld\n",numCols,numRows,data->minFlux,1,x->n); if (!pmReadoutFakeFromVectors(readout, numCols, numRows, x, y, mag, xOffset, yOffset, psf, data->minFlux, 0, false, true)) { @@ -168,5 +179,5 @@ } } - + psFree(sources); psFree(xOffset); Index: /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImageInterpolate.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImageInterpolate.c (revision 34838) @@ -31,5 +31,5 @@ #include "psImageConvolve.h" -# define IS_BILIN_SEPARABLE 0 +# define IS_BILIN_SEPARABLE 1 #include "psImageInterpolate.h" @@ -189,4 +189,5 @@ switch (mode) { case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR_SIMPLE: // Nothing to pre-compute break; @@ -287,6 +288,8 @@ # endif case PS_INTERPOLATE_BIQUADRATIC: + case PS_INTERPOLATE_BILINEAR_SIMPLE: // 2D kernel, would cost too much memory to pre-calculate break; + default: psAbort("Unsupported interpolation mode: %x", mode); @@ -937,4 +940,59 @@ +psImageInterpolateStatus interpolateJustWork(double *imageValue, double *varianceValue, + psImageMaskType *maskValue, float x, float y, + const psImageInterpolation *interp) { + float u1,u2; + int xl,xh; + int yl,yh; + const psImage *image = interp->image; // Image to interpolate + + xl = floor(x); + if (xl < 0) { xl = 0; } + if (xl >= image->numCols) {xl = image->numCols - 1; } + xh = xl + 1; + if (xh >= image->numCols) {xh = image->numCols - 1; } + + yl = floor(y); + if (yl < 0) { yl = 0; } + if (yl >= image->numRows) {yl = image->numRows - 1; } + yh = yl + 1; + if (yh >= image->numRows) {yh = image->numRows - 1; } + + if (imageValue && image) { + if (image->data.F32[yl][xl] == image->data.F32[yl][xh]) { + u1 = image->data.F32[yl][xh]; + } + else { + u1 = (xh - x) * image->data.F32[yl][xl] + (x - xl) * image->data.F32[yl][xh]; + } + if (image->data.F32[yh][xl] == image->data.F32[yh][xh]) { + u2 = image->data.F32[yh][xh]; + } + else { + u2 = (xh - x) * image->data.F32[yh][xl] + (x - xl) * image->data.F32[yh][xh]; + } + if (u1 == u2) { + *imageValue = u1; + } + else { + *imageValue = (yh - y) * u1 + (y - yl) * u2; + } + if ((floor(x) >= image->numCols)|| + (floor(y) >= image->numRows)) { + *imageValue = NAN; + } + } +#if (0) + if (*imageValue == 0.0) { + fprintf(stderr,"IJK: Zero!: %g %g [%d %d %d %d] %g %g (%g %g %g %g)\n", + x,y,xl,xh,yl,yh,u1,u2, + image->data.F32[yl][xl],image->data.F32[yl][xh],image->data.F32[yh][xl],image->data.F32[yh][xh] + ); + } +#endif + return PS_INTERPOLATE_STATUS_GOOD; +} + psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, @@ -969,7 +1027,9 @@ case PS_INTERPOLATE_BIQUADRATIC: return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); + case PS_INTERPOLATE_BILINEAR_SIMPLE: + return interpolateJustWork(imageValue, varianceValue, maskValue, x, y, interp); case PS_INTERPOLATE_BILINEAR: # if (!IS_BILIN_SEPARABLE) - return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); + return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp); # endif case PS_INTERPOLATE_GAUSS: @@ -1019,4 +1079,5 @@ switch (mode) { case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR_SIMPLE: // No smearing by design return 1.0; @@ -1083,4 +1144,5 @@ if (!strcasecmp(name, "FLAT")) return PS_INTERPOLATE_FLAT; if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR; + if (!strcasecmp(name, "SIMPLEBILINEAR")) return PS_INTERPOLATE_BILINEAR_SIMPLE; if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC; if (!strcasecmp(name, "GAUSS")) return PS_INTERPOLATE_GAUSS; Index: /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImageInterpolate.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImageInterpolate.h (revision 34838) @@ -33,4 +33,5 @@ PS_INTERPOLATE_LANCZOS3, ///< Sinc interpolation with 6x6 pixel kernel PS_INTERPOLATE_LANCZOS4, ///< Sinc interpolation with 8x8 pixel kernel + PS_INTERPOLATE_BILINEAR_SIMPLE, ///< Simple manual bilinear interpolation } psImageInterpolateMode; Index: /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImagePixelManip.c =================================================================== --- /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImagePixelManip.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psLib/src/imageops/psImagePixelManip.c (revision 34838) @@ -216,15 +216,54 @@ } - +#define psImageOverlayLoopClean(DATATYPE,OP) { \ + for (int row=y0;rowdata.DATATYPE[row]; \ + ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \ + for (int col=x0;coldata.DATATYPE[row]; \ + ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \ + for (int col=x0;coldata.DATATYPE[row]; \ ps##DATATYPE* overlayRow = overlay->data.DATATYPE[row-y0]; \ - for (int col=x0;coldata.PS_TYPE_IMAGE_MASK_DATA[y][x] != 0)) { + for (int j = 0; j < Nbits; j++) { + psImageMaskType M = (psImageMaskType) pow(2,j); + if (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & M) { + bitcounter[j]++; + } + } + } + if (inMask && (inMask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { continue; @@ -65,7 +84,10 @@ sum += inImage->data.F32[y][x]; numPix++; + + } } - + + // Values to set float imageValue; @@ -79,5 +101,11 @@ } outImage->data.F32[yOut][xOut] = imageValue; - outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue; + // outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = maskValue; + outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] = 0; + for (int j = 0; j < Nbits; j++) { + if (bitcounter[j] > 0.5 * pxlcount) { + outMask->data.PS_TYPE_IMAGE_MASK_DATA[yOut][xOut] |= (int) pow(2,j); + } + } xStart = xStop; } Index: /branches/eam_branches/ipp-20121130/psModules/src/camera/pmFPAfileIO.c =================================================================== --- /branches/eam_branches/ipp-20121130/psModules/src/camera/pmFPAfileIO.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psModules/src/camera/pmFPAfileIO.c (revision 34838) @@ -76,5 +76,4 @@ while ((item = psMetadataGetAndIncrement (iter)) != NULL) { pmFPAfile *file = item->data.V; - switch (place) { case PM_FPA_BEFORE: Index: /branches/eam_branches/ipp-20121130/psModules/src/detrend/pmDark.c =================================================================== --- /branches/eam_branches/ipp-20121130/psModules/src/detrend/pmDark.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psModules/src/detrend/pmDark.c (revision 34838) @@ -353,12 +353,35 @@ // retrieve the required parameter vectors - psArray *values = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES"); - psAssert(values, "values not supplied"); + psArray *in_values = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES"); + psAssert(in_values, "values not supplied"); psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert(roMask, "roMask not supplied"); - psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS"); - psAssert(orders, "orders not supplied"); - - psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting + psVector *max_orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS"); + psAssert(max_orders, "orders not supplied"); + + psArray *values_set = psArrayAlloc(max_orders->n); + psArray *poly_set = psArrayAlloc(max_orders->n); + psVector *logL = psVectorAlloc(max_orders->n,PS_TYPE_F64); + + for (int i = 0; i < max_orders->n; i++) { + psVector *orders = psVectorAlloc(i+1,PS_TYPE_U8); + for (int j = 0; j < orders->n; j++) { + orders->data.U8[j] = max_orders->data.U8[j]; + } + poly_set->data[i] = psPolynomialMDAlloc(orders); // Polynomial for fitting + + psArray *values = psArrayAlloc(in_values->n); + + for (int j = 0; j < values->n; j++) { + psVector *these_values = psVectorAlloc(i+1,PS_TYPE_F32); + psVector *input_values = in_values->data[j]; + + for (int k = 0; k < orders->n; k++) { + these_values->data.F32[k] = input_values->data.F32[k]; + } + values->data[j] = these_values; + } + values_set->data[i] = values; + } // retrieve the norm vector, if supplied @@ -383,5 +406,5 @@ } - pmDarkVisualInit(values); + pmDarkVisualInit(values_set->data[max_orders->n - 1]); pmReadout *outReadout = output->readouts->data[0]; @@ -423,15 +446,53 @@ } - if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) { + int k_best = 0; + for (int k = 0; k < max_orders->n; k++) { + psPolynomialMD *poly = poly_set->data[k]; + psArray *values = values_set->data[k]; + + if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) { psErrorClear(); // Nothing we can do about it psVectorInit(poly->coeff, NAN); - } - - pmDarkVisualPixelFit(pixels, mask); - pmDarkVisualPixelModel(poly, values); - - for (int k = 0; k < poly->coeff->n; k++) { + } + + pmDarkVisualPixelFit(pixels, mask); + pmDarkVisualPixelModel(poly, values); + + // Insert math here to choose optimum model. + logL->data.F64[k] = 0.0; + psPolynomialMD *polySig = poly_set->data[0]; + for (int m = 0; m < poly->deviations->n; m++) { + logL->data.F64[k] += pow(poly->deviations->data.F32[m] / polySig->stdevFit,2); +/* if ((xOut == 20) && (yOut == 256)) { */ +/* psTrace("psModules.detrend",3,"pmDarkCombine DEV: %d %d: input %d models: Norders: %d logL: %g value: %g\n", */ +/* xOut,yOut,m,k,logL->data.F64[k],poly->deviations->data.F32[m]); */ +/* } */ + } + if (k > 0) { + if ( ( logL->data.F64[k - 1] - logL->data.F64[k] ) > 1) { // Hard coded criterion for a ~5% limit with one degree of freedom + k_best = k; + } + } + if ((xOut <= 600) && (yOut <= 600)) { + psTrace("psModules.detrend",3,"pmDarkCombine: %d %d: models: Norders: %d logL: %g BestOrders: %d\n", + xOut,yOut,k,logL->data.F64[k],k_best); + } + } + if (k_best > 1) { + k_best = 1; + } +/* k_best = 1; */ + // Select the polynomial that seems best. + psPolynomialMD *poly = poly_set->data[k_best]; + + // for (int k = 0; k < poly->coeff->n; k++) { + for (int k = 0; k < max_orders->n + 1; k++) { // There is one more coefficient than is stored here. pmReadout *ro = output->readouts->data[k]; // Readout of interest - ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k]; + if (k < poly->coeff->n) { + ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k]; + } + else { + ro->image->data.F32[yOut][xOut] = 0.0; + } } counts->data.U16[yOut][xOut] = poly->numFit; Index: /branches/eam_branches/ipp-20121130/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/eam_branches/ipp-20121130/psModules/src/imcombine/pmStack.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psModules/src/imcombine/pmStack.c (revision 34838) @@ -1358,4 +1358,57 @@ } +bool pmStackSimpleMedianCombine( + pmReadout *combined, + psArray *input) { + int num = input->n; + // int numCols, numRows; + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine + int xSize, ySize; // Size of the output image + + psArray *stack = psArrayAlloc(num); // Stack of readouts + for (int i = 0; i < num; i++) { + // pmStackData *data = input->data[i]; // Stack data for this input + pmReadout *ro = input->data[i]; // data->readout; // Readout of interest + if (!ro) { + continue; + } + stack->data[i] = psMemIncrRefCounter(ro); + } + + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, + stack)) { + psError(psErrorCodeLast(), false, "Input stack is not valid."); + psFree(stack); + return false; + } + + psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32); + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { + for (int i = 0; i < input->n; i++) { + pmReadout *ro = stack->data[i]; + psImage *image = ro->image; + pixelData->data.F32[i] = image->data.F32[y][x]; + } + if (!psVectorStats(stats,pixelData,NULL,NULL,0)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate median"); + psFree(stats); + psFree(pixelData); + psFree(stack); + return(false); + } + combined->image->data.F32[y][x] = stats->robustMedian; + combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + } + } + + psFree(stats); + psFree(pixelData); + psFree(stack); + return (true); +} + /// Stack input images bool pmStackCombine( Index: /branches/eam_branches/ipp-20121130/psModules/src/imcombine/pmStack.h =================================================================== --- /branches/eam_branches/ipp-20121130/psModules/src/imcombine/pmStack.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/psModules/src/imcombine/pmStack.h (revision 34838) @@ -41,4 +41,8 @@ float addVariance ///< Additional variance when rejecting ); +/// Stack input images simply +bool pmStackSimpleMedianCombine(pmReadout *combined, ///< Combined readout (output) + psArray *input ///< Input array of pmStackData + ); /// Stack input images Index: /branches/eam_branches/ipp-20121130/psModules/src/objects/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20121130/psModules/src/objects/Makefile.am (revision 34837) +++ /branches/eam_branches/ipp-20121130/psModules/src/objects/Makefile.am (revision 34838) @@ -102,5 +102,4 @@ pmSourceDiffStats.h \ pmSourceExtendedPars.h \ - pmSourceSatstar.h \ pmSourceUtils.h \ pmSourceSky.h \ @@ -111,4 +110,5 @@ pmSourceOutputs.h \ pmSourceIO.h \ + pmSourceSatstar.h \ pmSourcePlots.h \ pmSourceVisual.h \ Index: /branches/eam_branches/ipp-20121130/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- /branches/eam_branches/ipp-20121130/psModules/src/objects/pmFootprintCullPeaks.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psModules/src/objects/pmFootprintCullPeaks.c (revision 34838) @@ -91,5 +91,16 @@ // max flux is above threshold for brightest peak - pmPeak *maxPeak = fp->peaks->data[0]; + pmPeak *maxPeak = NULL; + for (int i = 0; i < fp->peaks->n; i++) { + pmPeak *testPeak = fp->peaks->data[i]; + float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux; + + if (isfinite(this_peak)) { + maxPeak = fp->peaks->data[i]; + break; + } + } + psAssert(maxPeak,"maxPeak was not set in these peaks"); + // = fp->peaks->data[0]; float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux; Index: /branches/eam_branches/ipp-20121130/psphot/src/psphotOutput.c =================================================================== --- /branches/eam_branches/ipp-20121130/psphot/src/psphotOutput.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psphot/src/psphotOutput.c (revision 34838) @@ -201,8 +201,13 @@ int nCR = 0; int nEXT = 0; + int nForced = 0; + int nDetections = sources != NULL ? sources->n : 0; // count the number of sources which will be written and other sub-types for (int i = 0; (sources != NULL) && (i < sources->n); i++) { pmSource *source = (pmSource *) sources->data[i]; + if (source->mode2 & PM_SOURCE_MODE2_MATCHED) { + nForced ++; + } pmModel *model = pmSourceGetModel (NULL, source); if (model == NULL) @@ -217,4 +222,6 @@ } + // XXX: This loop cauess nSrc to be twice as large as it should be. This is kept + // for compatability for (int i = 0; (sources != NULL) && (i < sources->n); i++) { pmSource *source = (pmSource *) sources->data[i]; @@ -226,6 +233,8 @@ psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", PS_META_REPLACE, "Number of sources", nSrc); + psMetadataAddS32 (header, PS_LIST_TAIL, "NDET", PS_META_REPLACE, "Number of detections", nDetections); psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT); psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_CR", PS_META_REPLACE, "Number of cosmic rays", nCR); + psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_FRC", PS_META_REPLACE, "Number of Forced detections", nForced); return true; } @@ -306,6 +315,8 @@ psMetadataItemSupplement (&status, header, analysis, "NSTARS"); + psMetadataItemSupplement (&status, header, analysis, "NDET"); psMetadataItemSupplement (&status, header, analysis, "NDET_EXT"); psMetadataItemSupplement (&status, header, analysis, "NDET_CR"); + psMetadataItemSupplement (&status, header, analysis, "NDET_FRC"); psMetadataItemSupplement (&status, header, analysis, "PSPHOT.PSF.APRESID"); Index: /branches/eam_branches/ipp-20121130/psphot/src/psphotSourceMatch.c =================================================================== --- /branches/eam_branches/ipp-20121130/psphot/src/psphotSourceMatch.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psphot/src/psphotSourceMatch.c (revision 34838) @@ -7,4 +7,6 @@ int order; bool matchAll; + long nSources; + float yRatioMax; char filterID[MATCH_INFO_FILTER_STR_LEN]; } psphotStackMatchInfo; @@ -218,7 +220,49 @@ pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); + psAssert (detections->allSources, "all sources not defined?"); detArrays->data[i] = psMemIncrRefCounter(detections); readouts->data[i] = psMemIncrRefCounter(readout); + } + + // We are having a problem with large number of false y band detections causing large numbers of + // matched detections to be created. + // The matchInfo recipe each entry contains a value for Y.RATIO.MAX. For each filter with a non-zero value + // we compare the ratio of the number of y detections to the number of detections in that filter to Y.RATIO.MAX. + // If the ratio in all filters with a non-zero is above its limit we consider the y band stack to + // be suspect and so disable the y band matchAll recipe (if set) + + int iy = -1; + // note now we are looping over matchInfo + for (int i = 0; i < nImages; i++) { + // Remember which input is y + if (!strcmp(matchInfo[i].filterID, "y")) { + iy = i; + } + pmDetections *detections = detArrays->data[matchInfo[i].inputNum]; + matchInfo[i].nSources = detections->allSources->n; + } + + if (iy >= 0) { + int nCut = 0; + int nTry = 0; + for (int i = 0; i < nImages; i++) { + if (i == iy) continue; + if (matchInfo[i].yRatioMax != 0.0 && matchInfo[i].nSources > 0) { + ++nTry; + float ratio = (float) matchInfo[iy].nSources / (float) matchInfo[i].nSources; + psLogMsg ("psphot", PS_LOG_DETAIL, "nSrc_y / nSrc_%s: %6.3f max: %6.3f", matchInfo[i].filterID, ratio, + matchInfo[i].yRatioMax); + if (ratio > matchInfo[i].yRatioMax) { + // This one is above it's threshold + ++nCut; + } + } + } + if (nTry > 0 && nCut == nTry) { + // All filters with a cut were above the threshold. Disable matchAll for y band + psLogMsg ("psphot", PS_LOG_INFO, "Clearing y band matchAll because y ratio was too high in %d filters", nTry); + matchInfo[iy].matchAll = false; + } } @@ -516,4 +560,7 @@ bool matchAll = (strcasecmp("true", matchAllStr) == 0); + psString yRatioStr = psMetadataLookupStr (&status, item->data.md, "Y.RATIO.MAX"); + psF32 yRatioMax = atof(yRatioStr); + // look for this filter in the inputs for (int i = 0; i < num; i++) { @@ -524,4 +571,6 @@ matchInfo[i].order = order; matchInfo[i].matchAll = matchAll; + matchInfo[i].yRatioMax = yRatioMax; + matchInfo[i].nSources = 0; psLogMsg ("psphot", PS_LOG_DETAIL, "input: %d %s match order: %d match all: %d\n", i, thisFilter, order, matchAll); Index: /branches/eam_branches/ipp-20121130/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/ipp-20121130/psphot/src/psphotSourceSize.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/psphot/src/psphotSourceSize.c (revision 34838) @@ -11,4 +11,6 @@ float nSigmaMoments; float nSigmaCR; + bool altDiffExt; + float altDiffExtThresh; float soft; int grow; @@ -120,4 +122,11 @@ assert (status); + // Optional extended source measurement algorithm to improve diff image trails + options.altDiffExt = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.DIFF.ALTERNATE"); + assert (status); + // Threshold for this alternate method + options.altDiffExtThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.DIFF.ALTERNATE.THRESH"); + assert (status); + // location of a single test source options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST"); @@ -601,5 +610,24 @@ continue; } - + // Alternate extended source limit calculation + if (options->altDiffExt) { + // ratio of major to minor axes + // MRV = Major / Minor + // MRV = (0.5 * (Mxx + Myy) + 0.5 * sqrt( (Mxx + Myy)^2 + 4 Mxy^2)) / + // (0.5 * (Mxx + Myy) - 0.5 * sqrt( (Mxx + Myy)^2 + 4 Mxy^2)) + // MRV = (2 * 0.5 * (Mxx + Myy) - Minor) / Minor + float momentRatioVeres = (Mxx + Myy - Mminor) / Mminor; + bool isAltEXT = (momentRatioVeres > options->altDiffExtThresh); + if (isAltEXT) { + psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g ALTEXT\t%g %g\n", + source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid, + momentRatioVeres,options->altDiffExtThresh); + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + Next ++; + continue; + } + } + // Everything else should just be treated as a PSF psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g PSF\n", Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarp.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarp.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarp.c (revision 34838) @@ -46,4 +46,9 @@ goto DIE; } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + if (!pswarpDefineBackground(config)) { + goto DIE; + } + } // Open the statistics file @@ -66,4 +71,10 @@ if (!pswarpLoop(config, stats)) { goto DIE; + } + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND.MODEL")) { + if (!pswarpLoopBackground(config, stats)) { + fprintf(stderr,"Dying!\n"); + goto DIE; + } } Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarp.h =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarp.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarp.h (revision 34838) @@ -62,4 +62,9 @@ psImage *region; + /** values which are needed to control the background model warping. */ + bool background_warping; + double offset_x; + double offset_y; + /** input values for this tile */ int gridX; @@ -81,5 +86,7 @@ bool pswarpParseCamera (pmConfig *config); bool pswarpDefine (pmConfig *config); +bool pswarpDefineBackground (pmConfig *config); bool pswarpLoop (pmConfig *config, psMetadata *stats); +bool pswarpLoopBackground (pmConfig *config, psMetadata *stats); psExit pswarpExitCode(psExit exitValue); bool pswarpTransformReadout (pmReadout *output, pmReadout *input, pmConfig *config); Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpArguments.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpArguments.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpArguments.c (revision 34838) @@ -86,5 +86,6 @@ pmConfigFileSetsMD (config->arguments, &argc, argv, "MASK", "-mask", "-masklist"); pmConfigFileSetsMD (config->arguments, &argc, argv, "VARIANCE", "-variance", "-variancelist"); - + pmConfigFileSetsMD (config->arguments, &argc, argv, "BACKGROUND", "-background", "-bkglist"); + if (argc != 3) { usage(); @@ -164,4 +165,21 @@ } + bool doBKG = psMetadataLookupBool(&status,recipe, "BACKGROUND.MODEL"); ///< Generate the warped background model? + if (!status) { + doBKG = false; + psWarning("BACKGROUND.MODEL is not set in the %s recipe -- defaulting to FALSE.", PSWARP_RECIPE); + } + int bkgXgrid = psMetadataLookupS32(&status,recipe, "BKG.XGRID"); ///< Xsize of background model + if (!status) { + bkgXgrid = 10; + psWarning("BKG.XGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgXgrid); + } + int bkgYgrid = psMetadataLookupS32(&status,recipe, "BKG.YGRID"); ///< Xsize of background model + if (!status) { + bkgYgrid = 10; + psWarning("BKG.YGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgYgrid); + } + + // Set recipe values in the recipe (since we've possibly altered some) psMetadataAddS32(recipe, PS_LIST_TAIL, "GRID.NX", PS_META_REPLACE, @@ -176,5 +194,9 @@ "Fraction of bad flux for a pixel to be marked as poor", poorFrac); psMetadataAddBool(recipe, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF); - + psMetadataAddBool(recipe, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG); + psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid); + psMetadataAddS32(recipe, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid); + + // Set recipe values in the arguments psMetadataAddS32(config->arguments, PS_LIST_TAIL, "GRID.NX", 0, @@ -189,4 +211,7 @@ "Fraction of bad flux for a pixel to be marked as poor", poorFrac); psMetadataAddBool(config->arguments, PS_LIST_TAIL, "PSF", PS_META_REPLACE, "Generate a PSF Model?", PSF); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "BACKGROUND.MODEL", PS_META_REPLACE, "Generate the warped background model?", doBKG); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.XGRID", PS_META_REPLACE, "Xsize of background model", bkgXgrid); + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BKG.YGRID", PS_META_REPLACE, "Ysize of background model", bkgYgrid); psTrace("pswarp", 1, "Done with pswarpArguments...\n"); Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpDefine.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpDefine.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpDefine.c (revision 34838) @@ -41,5 +41,5 @@ return false; } - + // open the full skycell file; no need to defer different depths. only load the header data pmFPAview *view = pmFPAviewAlloc (0); @@ -86,4 +86,6 @@ psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0"); psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0"); + + } @@ -123,9 +125,153 @@ } } - + view->chip = view->cell = view->readout = -1; pmFPAAddSourceFromView(output->fpa, view, output->format); + psFree (view); return true; } + +bool pswarpDefineBackground (pmConfig *config) { + + // load the PSWARP recipe + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, true, "Can't find PSWARP recipe!\n"); + return false; + } + + // select the input data sources + pmFPAfile *skycell = psMetadataLookupPtr (NULL, config->files, "PSWARP.SKYCELL"); + if (!skycell) { + psError(PSWARP_ERR_CONFIG, false, "Can't find skycell data!\n"); + return false; + } + pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PSWARP.OUTPUT.BKGMODEL"); + if (!output) { + psError(PSWARP_ERR_CONFIG, false, "Can't find output data!\n"); + return false; + } + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.INPUT"); + if (!input) { + psError(PSWARP_ERR_CONFIG, false, "Can't find input data!\n"); + return false; + } + + // open the full skycell file; no need to defer different depths. only load the header data + pmFPAview *view = pmFPAviewAlloc (0); + pmFPAfileOpen (skycell, view, config); + + // Read header and create target + { + if (!pmFPAReadHeaderSet(skycell->fpa, skycell->fits, config)) { + psError(psErrorCodeLast(), false, "Unable to read headers for skycell."); + psFree(view); + return false; + } + view->chip = 0; + view->cell = 0; + view->readout = 0; + pmCell *source = pmFPAfileThisCell(config->files, view, "PSWARP.SKYCELL"); ///< Source cell + pmHDU *hdu = pmHDUFromCell(source); ///< HDU for source + if (!hdu || !hdu->header) { + psError(PM_ERR_PROG, false, "Unable to find header for sky cell."); + psFree(view); + return false; + } + int numCols = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); ///< Number of columns + int numRows = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); ///< Number of rows + if ((numCols == 0) || (numRows == 0)) { + psError(PSWARP_ERR_DATA, false, "skycell has invalid dimensions %d x %d", numCols, numRows); + psFree(view); + return false; + } + + pmCell *target = pmFPAviewThisCell(view, output->fpa); ///< Target cell + pmReadout *readout = pmReadoutAlloc(target); ///< Target readout + readout->image = psImageAlloc(numCols / output->xBin, numRows / output->yBin, PS_TYPE_F32); + psImageInit(readout->image, NAN); + psFree(readout); // Drop reference + + bool status = false; + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XBIN"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YBIN"); +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XBIN",PS_META_REPLACE,"",output->xBin); */ +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YBIN",PS_META_REPLACE,"",output->yBin); */ + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XSIZE"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YSIZE"); +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.XSIZE",PS_META_REPLACE,"",numCols / output->xBin); */ +/* psMetadataAddS32(target->concepts,PS_LIST_TAIL,"CELL.YSIZE",PS_META_REPLACE,"",numRows / output->yBin); */ + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.XPARITY"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.YPARITY"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.X0"); + psMetadataItemSupplement(&status, target->concepts, source->concepts, "CELL.Y0"); + + + } + + // XXX this is not a sufficient test + view->chip = 0; + view->cell = 0; + view->readout = -1; + pmHDU *phu = pmFPAviewThisPHU(view, skycell->fpa); ///< Skycell PHU + pmHDU *hdu = pmFPAviewThisHDU(view, skycell->fpa); ///< Skycell header + + pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); + + double cd1f = 1.0 * output->xBin; + double cd2f = 1.0 * output->yBin; + + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + + WCS->trans->x->coeff[1][0] *= cd1f; + WCS->trans->x->coeff[0][1] *= cd2f; + WCS->trans->y->coeff[1][0] *= cd1f; + WCS->trans->y->coeff[0][1] *= cd2f; + + + pmAstromWCStoHeader (hdu->header,WCS); + + bool bilevelAstrometry = false; + if (phu) { + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp(&ctype[4], "-DIS"); + } + } + + // We read from the skycell into the output. i.e., the output receives the desired astrometry. + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output + if (bilevelAstrometry) { + if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for skycell."); + psFree(view); + return false; + } + if (!pmAstromReadBilevelChip(outputChip, hdu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for skycell."); + psFree(view); + return false; + } + } else { + // we use a default FPA pixel scale of 1.0 + if (!pmAstromReadWCS(output->fpa, outputChip, hdu->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for skycell."); + psFree(view); + return false; + } + } + + + view->chip = view->cell = view->readout = -1; + pmFPAAddSourceFromView(output->fpa, view, output->format); + + + psFree (view); + return true; +} Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpFileNames.h =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpFileNames.h (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpFileNames.h (revision 34838) @@ -12,4 +12,5 @@ "PSWARP.MASK", "PSWARP.VARIANCE", + "PSWARP.BKGMODEL", NULL }; @@ -20,4 +21,5 @@ "PSWARP.OUTPUT.MASK", "PSWARP.OUTPUT.VARIANCE", + "PSWARP.OUTPUT.BKGMODEL", NULL }; Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpLoop.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpLoop.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpLoop.c (revision 34838) @@ -179,4 +179,6 @@ // read WCS data from the corresponding header pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); + + if (bilevelAstrometry) { if (!pmAstromReadBilevelChip (chip, hdu->header)) { @@ -195,5 +197,5 @@ } } - + pmCell *cell; while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { @@ -227,5 +229,5 @@ pswarpTransformReadout(output, readout, config); - + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { psError(psErrorCodeLast(), false, "Unable to write files."); @@ -255,5 +257,5 @@ goto DONE; } - + pmCell *outCell = output->parent; ///< Output cell pmChip *outChip = outCell->parent; ///< Output chip @@ -450,2 +452,354 @@ return true; } + +// Loop over the inputs, warp them to the output skycell and then write out the output. +bool pswarpLoopBackground(pmConfig *config, psMetadata *stats) +{ + bool status; + bool mdok; // Status of MD lookup + const char *skyCamera = psMetadataLookupStr(NULL, config->arguments, + "SKYCELL.CAMERA"); // Name of camera for skycell + pmConfigCamerasCull(config, skyCamera); + pmConfigRecipesCull(config, "PSWARP,PPSTATS,PSPHOT,PSASTRO,MASKS,JPEG"); + + // load the recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSWARP_RECIPE); + if (!recipe) { + psError(PSWARP_ERR_CONFIG, false, "missing recipe %s", PSWARP_RECIPE); + return false; + } + + if (!pswarpSetMaskBits(config)) { + psError(psErrorCodeLast(), false, "failed to set mask bits"); + return NULL; + } + + // select the input data sources + pmFPAfile *input = psMetadataLookupPtr(NULL, config->files, "PSWARP.BKGMODEL"); + if (!input) { + psError(PSWARP_ERR_CONFIG, true, "Can't find input data!\n"); + return false; + } + + // use the external astrometry source if supplied + pmFPAfile *astrom = psMetadataLookupPtr(NULL, config->files, "PSWARP.ASTROM"); + if (!astrom) { + astrom = input; + } + + if (astrom->camera != input->camera) { + psError(PSWARP_ERR_DATA, true, "Input camera and astrometry camera do not match."); + return false; + } + + // select the output readout + pmFPAview *view = pmFPAviewAlloc(0); + view->chip = 0; + view->cell = 0; + view->readout = 0; + pmReadout *output = pmFPAfileThisReadout(config->files, view, "PSWARP.OUTPUT.BKGMODEL"); + if (!output) { + psError(PSWARP_ERR_CONFIG, true, "Can't find output background data!\n"); + return false; + } + psFree (view); + // Turn all skycell files on to generate them, and then turn them off for the loop over the input images + // the input, which is in a different format. + { + pswarpFileActivation(config, detectorFiles, false); + pswarpFileActivation(config, photFiles, false); + pswarpFileActivation(config, independentFiles, false); + pswarpFileActivation(config, skycellFiles, true); + if (!pswarpIOChecksBefore(config)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + pswarpFileActivation(config, skycellFiles, false); + } + // Read the input astrometry + // XXX rather than use the activations here, this should just explicitly loop over the desired filerule + { + + pmFPAfileActivate(config->files, true, "PSWARP.ASTROM"); + + pmChip *chip; + pmFPAview *view = pmFPAviewAlloc(0); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { +#if 0 + // This needs to be removed because otherwise it throws an error of duplicate PSPHOT.DETECTIONS. + if (!chip->process || !chip->file_exists) { continue; } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } +#endif + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "ACell %d: %x %x %d\n", view->cell, cell->file_exists, cell->process,psErrorCodeLast()); + if (!cell->process || !cell->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE) || + !pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto DONE; + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto DONE; + } + psFree(view); + pswarpFileActivation(config, detectorFiles, true); + pmFPAfileActivate(config->files, false, "PSWARP.ASTROM"); + } + + // Don't care about the skycell anymore --- we've read it, and that's all we need to do. + pmFPAfileActivate(config->files, false, "PSWARP.SKYCELL"); + view = pmFPAviewAlloc(0); + + // find the FPA phu + bool bilevelAstrometry = false; + pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa); + + // pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header); + + if (phu) { + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + if (bilevelAstrometry) { + if (!pmAstromReadBilevelMosaic(input->fpa, phu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA."); + psFree(view); + psFree(stats); + goto DONE; + } + } + + psList *cells = psListAlloc(NULL); // List of cells, for concepts averaging + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + pmChip *chip; + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "DChip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + + // Pull information from the header of the background files so we can use it to set values. + pmHDU *hdu = pmFPAviewThisHDU(view,input->fpa); + psMetadata *header = hdu->header; + + int IMAXIS1 = psMetadataLookupS32(NULL,header,"IMNAXIS1"); + int IMAXIS2 = psMetadataLookupS32(NULL,header,"IMNAXIS2"); + int NAXIS1 = psMetadataLookupS32(NULL,header,"NAXIS1"); + int NAXIS2 = psMetadataLookupS32(NULL,header,"NAXIS2"); + char *CCDSUM = psMetadataLookupStr(NULL,header,"CCDSUM"); + int CCDSUM1 = atoi(strtok(CCDSUM," ")); + int CCDSUM2 = atoi(strtok(NULL," ")); + + psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_XOFFSET", PS_META_REPLACE, + "xoffset for background model data", (NAXIS1 * CCDSUM1 - IMAXIS1) / (2.0 * CCDSUM1)); + psMetadataAddF32(config->arguments,PS_LIST_TAIL,"BKG_WARP_YOFFSET", PS_META_REPLACE, + "yoffset for background model data", (NAXIS2 * CCDSUM2 - IMAXIS2) / (2.0 * CCDSUM2)); + psTrace("pswarp",5,"%d %d %d %d %d %d %g %g %d %d", + psMetadataLookupS32(NULL,header,"IMNAXIS1"), + psMetadataLookupS32(NULL,header,"IMNAXIS2"), + psMetadataLookupS32(NULL,header,"NAXIS1"), + psMetadataLookupS32(NULL,header,"NAXIS2"), + CCDSUM1,CCDSUM2, + psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"), + psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"), + psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID"), + psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID")); + + + // read WCS data from the corresponding header + hdu = pmFPAviewThisHDU (view, astrom->fpa); + + pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); + + double cd1f = (1.0 * CCDSUM1);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.XGRID")); + double cd2f = (1.0 * CCDSUM2);// * (1.0 * psMetadataLookupS32(NULL,config->arguments,"BKG.YGRID")); + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f; + WCS->crpix2 = WCS->crpix2 / cd2f; + + // WCS->trans->x->nX/nY + for (int q = 0; q <= WCS->trans->x->nX; q++) { + for (int r = 0; r <= WCS->trans->x->nY; r++) { + WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + } + } + for (int q = 0; q <= WCS->trans->y->nX; q++) { + for (int r = 0; r <= WCS->trans->y->nY; r++) { + WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + } + } + pmAstromWCStoHeader (hdu->header,WCS); + + if (bilevelAstrometry) { + if (!pmAstromReadBilevelChip (chip, hdu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel chip astrometry for input FPA."); + psFree(view); + psFree(stats); + goto DONE; + } + } else { + // we use a default FPA pixel scale of 1.0 + if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); + psFree(view); + psFree(stats); + goto DONE; + } + } + // Modify structure here. + + + pmCell *cell; + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psTrace ("pswarp", 4, "DCell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + + psListAdd(cells, PS_LIST_TAIL, cell); + + // process each of the readouts + pmReadout *readout; + while ((readout = pmFPAviewNextReadout(view, input->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Unable to read files."); + goto DONE; + } + if (!readout->data_exists) { + continue; + } + + for (int x = 0; x < readout->image->numCols; x++) { + for (int y = 0; y < readout->image->numRows; y++) { + readout->image->data.F32[y][x] = readout->image->data.F32[y][x] * (cd1f * cd2f) / + (psMetadataLookupS32(&mdok,config->arguments,"BKG.XGRID") * + psMetadataLookupS32(&mdok,config->arguments,"BKG.YGRID")); + } + } + + psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8); + + psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", true); + pswarpTransformReadout(output, readout, config); + psMetadataAddBool(config->arguments,PS_LIST_TAIL, "BACKGROUND_WARPING", PS_META_REPLACE, "", false); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto DONE; + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto DONE; + } + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto DONE; + } + } + if (!output->data_exists) { + psWarning("No overlap between input and skycell."); + psphotFilesActivate(config, false); + psFree(cells); + psFree(view); + goto DONE; + } + pmCell *outCell = output->parent; ///< Output cell + pmChip *outChip = outCell->parent; ///< Output chip + pmFPA *outFPA = outChip->parent; ///< Output FP + +/* if (!pswarpPixelsLit(output, stats, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to calculate pixel regions."); */ +/* psFree(cells); */ +/* psFree(view); */ +/* goto DONE; */ +/* } */ + psRegion *trimsec = psMetadataLookupPtr(NULL, outCell->concepts, "CELL.TRIMSEC"); ///< Trim section + trimsec->x0 = trimsec->x1 = trimsec->y0 = trimsec->y1 = 0; ///< All pixels + + if (!psMetadataCopy(outFPA->concepts, input->fpa->concepts)) { + psError(psErrorCodeLast(), false, "Unable to copy FPA concepts from input to output."); + psFree(stats); + psFree(view); + goto DONE; + } + + pmHDU *hdu = outFPA->hdu; ///< HDU for the output warped image + + // Copy header from target + { + pmFPAview *skyView = pmFPAviewAlloc(0); ///< View into skycell + skyView->chip = skyView->cell = 0; + pmCell *cell = pmFPAfileThisCell(config->files, skyView, "PSWARP.SKYCELL"); // Skycell cell + psFree(skyView); + pmHDU *skyHDU = pmHDUFromCell(cell); ///< HDU + if (!skyHDU) { + psError(PSWARP_ERR_DATA, false, "Unable to find skycell HDU."); + psFree(view); + goto DONE; + } + hdu->header = psMetadataCopy(hdu->header, skyHDU->header); + } + pswarpVersionHeader(hdu->header); + + if (!pmAstromWriteWCS(hdu->header, outFPA, outChip, WCS_NONLIN_TOL)) { + psError(psErrorCodeLast(), false, "Unable to generate WCS header."); + psFree(stats); + goto DONE; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Unable to write files."); + goto DONE; + } + // Done with the detector side of things + pswarpFileActivation(config, detectorFiles, false); + pswarpFileActivation(config, independentFiles, false); + + + // Add MD5 information for readout + const char *chipName = psMetadataLookupStr(NULL, output->parent->parent->concepts, "CHIP.NAME"); + const char *cellName = psMetadataLookupStr(NULL, output->parent->concepts, "CELL.NAME"); + psString headerName = NULL; ///< Header name for MD5 + psStringAppend(&headerName, "MD5_%s_%s_%d", chipName, cellName, view->readout); + psVector *md5 = psImageMD5(output->image); ///< md5 hash + psString md5string = psMD5toString(md5); ///< String + psFree(md5); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, headerName, PS_META_REPLACE, + "Image MD5", md5string); + psFree(md5string); + psFree(headerName); + psFree(view); + DONE: + + return true; +} Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpParseCamera.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpParseCamera.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpParseCamera.c (revision 34838) @@ -84,4 +84,10 @@ } + pmFPAfile *inBackground = defineInputFile(config, NULL, "PSWARP.BKGMODEL", "BACKGROUND", + PM_FPA_FILE_IMAGE); + if (!inBackground) { + psLogMsg("pswarp", PS_LOG_INFO, "No background models supplied"); + } + // The input skycell is a required argument: it defines the output image // XXX we may need a different skycell structure here @@ -134,4 +140,20 @@ } + if (inBackground && psMetadataLookupBool(&mdok, recipe, "BACKGROUND.MODEL")) { + pmFPAfile *outBackground = pmFPAfileDefineSkycell(config, NULL, "PSWARP.OUTPUT.BKGMODEL"); +/* pmFPAfile *outBackground = pmFPAfileDefineFromFPA(config,output->fpa, */ +/* psMetadataLookupS32(&mdok,recipe,"BKG.XGRID"), */ +/* psMetadataLookupS32(&mdok,recipe,"BKG.YGRID"), */ +/* "PSWARP.OUTPUT.BKGMODEL"); */ + outBackground->xBin = psMetadataLookupS32(&mdok, recipe, "BKG.XGRID"); + outBackground->yBin = psMetadataLookupS32(&mdok, recipe, "BKG.YGRID"); + + if (!outBackground) { + psError(psErrorCodeLast(), false, "Failed to build FPA from PSWARP.OUTPUT.BKGMODEL"); + return false; + } + outBackground->save = true; + } + if (psMetadataLookupBool(&mdok, recipe, "PSF")) { // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformReadout.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformReadout.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformReadout.c (revision 34838) @@ -29,4 +29,5 @@ psImageInterpolateMode interpolationMode = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.MODE"); ///< Mode for interp + int numKernels = psMetadataLookupS32(NULL, config->arguments, "INTERPOLATION.NUM"); ///< Number of kernels @@ -60,5 +61,6 @@ // output coordinates to input coordinates pswarpMapGrid *grid = pswarpMapGridFromImage(input, output, nGridX, nGridY); - + // if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { + // XXX optionally modify the grid based on this result and force the maxError < XXX double maxError = pswarpMapGridMaxError(grid); // Maximum (positional) error from using grid @@ -130,4 +132,11 @@ args->goodPixels = 0; + if (psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { + args->background_warping = true; + args->offset_x = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_XOFFSET"); + args->offset_y = psMetadataLookupF32(NULL,config->arguments,"BKG_WARP_YOFFSET"); + } + + // allocate a job psThreadJob *job = psThreadJobAlloc ("PSWARP_TRANSFORM_TILE"); @@ -171,5 +180,4 @@ } else { pswarpTransformTileArgs *args = job->args->data[0]; - // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n); goodPixels += args->goodPixels; xMin = PS_MIN(args->xMin, xMin); @@ -183,4 +191,6 @@ jacobian += args->jacobian * args->goodPixels; } + + } psFree(job); @@ -201,8 +211,10 @@ if (goodPixels > 0 && psMetadataLookupBool(&mdok, recipe, "SOURCES")) { + if (!psMetadataLookupBool(NULL,config->arguments,"BACKGROUND_WARPING")) { if (!pswarpTransformSources(output, input, config)) { - psError(psErrorCodeLast(), false, "Unable to interpolate image."); - return false; + psError(psErrorCodeLast(), false, "Unable to interpolate image."); + return false; } + } } Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformSources.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformSources.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformSources.c (revision 34838) @@ -41,5 +41,5 @@ if (!outDetections) { outDetections = pmDetectionsAlloc(); - psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Warped sources", outDetections); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY , "Warped sources", outDetections); } psArray *outSources = outDetections->allSources; Index: /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformTile.c =================================================================== --- /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformTile.c (revision 34837) +++ /branches/eam_branches/ipp-20121130/pswarp/src/pswarpTransformTile.c (revision 34838) @@ -46,4 +46,8 @@ args->jacobian = NAN; + args->background_warping = false; + args->offset_x = 0.0; + args->offset_y = 0.0; + return args; } @@ -83,16 +87,32 @@ // Iterate over the output image pixels (parent frame) long goodPixels = 0; ///< Number of input pixels landing on the output image + for (int y = yMin; y < yMax; y++) { for (int x = xMin; x < xMax; x++) { - // Only transform those pixels requested if (region && region->data.PS_TYPE_IMAGE_MASK_DATA[y][x]) { continue; } - // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. // both are in the parent frames of the input and output images. double xIn, yIn; // Input pixel coordinates pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5); + + // This needs to use a more reliable method to do this offset and limiting + // if (args->interp->mode == 8) { + if (args->background_warping) { + // double xOffset = 177.0 / 400.0; // (modelsize * modelbinning - xsize) / 2.0 + // double yOffset = 166.0 / 400.0; // (modelsize * modelbinning - ysize) / 2.0 + xIn += args->offset_x; + yIn += args->offset_y; + + if ((xIn > inNumCols - args->offset_x)|| + (yIn > inNumRows - args->offset_y)|| + (xIn < args->offset_x)|| + (yIn < args->offset_y)) { + continue; + } + } + if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) { continue;