Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotArguments.c =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphotArguments.c (revision 42506) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotArguments.c (revision 42507) @@ -97,4 +97,6 @@ } +void psphotSaveConfig (pmConfig *config); + pmConfig *psphotArguments(int argc, char **argv) { @@ -244,4 +246,7 @@ psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + // save a copy for deep retrieval + psphotSaveConfig (config); + psTrace("psphot", 1, "Done with psphotArguments...\n"); return (config); Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceArguments.c =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceArguments.c (revision 42506) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceArguments.c (revision 42507) @@ -3,4 +3,5 @@ static void usage(const char *program, psMetadata *arg, pmConfig *config, int exitCode); static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile); +void psphotSaveConfig (pmConfig *config); pmConfig *psphotFullForceArguments(int argc, char **argv) { @@ -111,4 +112,6 @@ // output position is fixed psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + psphotSaveConfig (config); psTrace("psphot", 1, "Done with psphotFullForceArguments...\n"); Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceReadout.c =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceReadout.c (revision 42506) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotFullForceReadout.c (revision 42507) @@ -106,4 +106,23 @@ psphotKronFlux (config, view, filerule); + // Option to do the non-linear fitting for the brighter sources + if (1) { + // identify CRs and extended sources (only unmeasured sources are measured) + psphotSourceSize (config, view, filerule, true); // pass 1 (detections->allSources) + + // non-linear PSF and EXT fit to brighter sources + // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) + // XXX: can leave faulted job in done queue + psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources) + + // replace all sources + psphotReplaceAllSources (config, view, filerule, false); // pass 1 (detections->allSources) + + // linear fit to include all sources (subtract again) + // NOTE : apply to ALL sources (extended + psf) + // do NOT skip negative-flux sources (last arg = false) + psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources) + } + psphotChipParams (config, view, filerule); Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c (revision 42506) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotMakeResiduals.c (revision 42507) @@ -1,3 +1,5 @@ # include "psphotInternal.h" +pmConfig *psphotGetConfig (); +bool _psphotSaveCube (char *basename, char *extname, psArray *cube); # define RESIDUAL_SOFTENING 0.005 @@ -93,4 +95,12 @@ psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32); +# define TESTRESID 1 +# if (TESTRESID) + psArray *sourceRawCube = psArrayAllocEmpty (1); + psArray *sourceSigCube = psArrayAllocEmpty (1); + psArray *sourceVarCube = psArrayAllocEmpty (1); + psArray *sourceMskCube = psArrayAllocEmpty (1); +# endif + // build (DATA - MODEL) [an image] for each psf star psArray *input = psArrayAllocEmpty (100); @@ -108,4 +118,11 @@ psImage *mask = psImageCopy (NULL, source->maskView ? source->maskView : source->maskObj, PS_TYPE_IMAGE_MASK); psImage *variance = psImageCopy (NULL, source->variance ? source->variance : source->pixels, PS_TYPE_F32); + +# if (TESTRESID) + psImage *tmpImage = psImageCopy(NULL, image, PS_TYPE_F32); + psArrayAdd (sourceRawCube, 1, tmpImage); + psFree (tmpImage); +# endif + pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal); @@ -114,4 +131,10 @@ psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32)); psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32)); + +# if (TESTRESID) + psArrayAdd (sourceSigCube, 1, image); + psArrayAdd (sourceVarCube, 1, variance); + psArrayAdd (sourceMskCube, 1, mask); +# endif // we interpolate the image and variance - include the mask or not? @@ -141,4 +164,32 @@ pmResiduals *resid = pmResidualsAlloc (xSize, ySize, xBin, yBin); psImageInit (resid->mask, 0); + + +# if (TESTRESID) + psArray *inputSigCube = psArrayAllocEmpty (1); + psArray *inputVarCube = psArrayAllocEmpty (1); + psArray *inputMskCube = psArrayAllocEmpty (1); + + // below, we are generating the residual images with interpolation + // to save these in an output cube, I need to generate place-holder images here + for (int i = 0; i < input->n; i++) { + psImage *inputSig = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32); + psImage *inputVar = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32); + psImage *inputMsk = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_IMAGE_MASK); + + psImageInit (inputSig, NAN); + psImageInit (inputVar, NAN); + psImageInit (inputMsk, 0); + + psArrayAdd (inputSigCube, 1, inputSig); + psArrayAdd (inputVarCube, 1, inputVar); + psArrayAdd (inputMskCube, 1, inputMsk); + + // XXX free these here, right? + psFree (inputSig); + psFree (inputVar); + psFree (inputMsk); + } +# endif // x(resid) = (x(image) - Xo)*xBin + xCenter @@ -194,5 +245,15 @@ nGoodPixel ++; } + +# if (TESTRESID) + psImage *inputSig = inputSigCube->data[i]; inputSig->data.F32[oy][ox] = flux; + psImage *inputVar = inputVarCube->data[i]; inputVar->data.F32[oy][ox] = flux; + psImage *inputMsk = inputMskCube->data[i]; inputMsk->data.PS_TYPE_IMAGE_MASK_DATA[oy][ox] = fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]; +# endif } + + if (nGoodPixel < 0.05*input->n) { + fprintf (stderr, "warning: %d, %d : residual pixel with few input pixels: %d vs %d\n", ox, oy, nGoodPixel, (int) input->n); + } // skip pixels with insufficient data @@ -220,4 +281,8 @@ } + if (fabs(fluxClip->robustMedian) > 1.5) { + fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->robustMedian); + } + // mark input pixels which are more than N sigma from the median int nKeep = 0; @@ -233,4 +298,8 @@ if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++; } + + if (nKeep < 5) { + fprintf (stderr, "warning: %d, %d : residual pixel with few good pixels: %d vs %d\n", ox, oy, nKeep, (int) input->n); + } if (SPATIAL_ORDER == 0) { @@ -242,4 +311,5 @@ } + // XXX this test should go at the top of the loop float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); if (radius > radiusMax) { @@ -304,4 +374,8 @@ float dRo = sqrt(A->data.F32[0][0]); + if (fabs(B->data.F64[0]) > 1.5) { + fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, B->data.F64[0]); + } + if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) { resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; @@ -313,4 +387,43 @@ } } + +# if (TESTRESID) + + pmConfig *config = psphotGetConfig(); + char *output = psMetadataLookupStr (&status, config->arguments, "OUTPUT"); + + // XXX // save the star coordinate vectors + // XXX + // XXX psString filename = NULL; + // XXX psStringAppend (&filename, "%s.instar.dat", output); + + char filename[1024]; + sprintf (filename, "%s.instar.dat", output); + + FILE *fout = fopen (filename, "w"); + for (int i = 0; i < xC->n; i++) { + fprintf (fout, "%d %f %f\n", i, xC->data.F32[i], yC->data.F32[i]); + } + fclose (fout); + + _psphotSaveCube (output, "raw.raw.fits", sourceRawCube); + _psphotSaveCube (output, "sig.raw.fits", sourceSigCube); + _psphotSaveCube (output, "var.raw.fits", sourceVarCube); + _psphotSaveCube (output, "msk.raw.fits", sourceMskCube); + + _psphotSaveCube (output, "sig.fits", inputSigCube); + _psphotSaveCube (output, "var.fits", inputVarCube); + _psphotSaveCube (output, "msk.fits", inputMskCube); + + psFree (sourceRawCube); + psFree (sourceSigCube); + psFree (sourceVarCube); + psFree (sourceMskCube); + + psFree (inputSigCube); + psFree (inputVarCube); + psFree (inputMskCube); + +# endif psFree (A); @@ -343,2 +456,18 @@ return (resid != NULL) ? true : false; } + +bool _psphotSaveCube (char *basename, char *extname, psArray *cube) { + + psString filename = NULL; + psStringAppend (&filename, "%s.%s", basename, extname); + + psFits *fits = psFitsOpen (filename, "w"); + if (!psFitsWriteImageCube (fits, NULL, cube, NULL)) { + fprintf (stderr, "failed to write the cube %s\n", filename); + } + psFitsClose (fits); + + psFree (filename); + + return true; +} Index: /branches/eam_branches/ipp-20230313/psphot/src/psphotOutput.c =================================================================== --- /branches/eam_branches/ipp-20230313/psphot/src/psphotOutput.c (revision 42506) +++ /branches/eam_branches/ipp-20230313/psphot/src/psphotOutput.c (revision 42507) @@ -1,3 +1,13 @@ # include "psphotInternal.h" + +pmConfig *staticConfig = NULL; + +void psphotSaveConfig (pmConfig *config) { + staticConfig = config; +} + +pmConfig *psphotGetConfig () { + return staticConfig; +} // convert filerule to filerule.NUM and look up in the config->arguments metadata