Index: trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- trunk/psphot/src/psphotMakeResiduals.c (revision 35559) +++ trunk/psphot/src/psphotMakeResiduals.c (revision 42842) @@ -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,14 @@ psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32); +# define TESTRESID 0 +# if (TESTRESID) + psArray *sourceRawCube = psArrayAllocEmpty (1); + psArray *sourceSigCube = psArrayAllocEmpty (1); + psArray *sourceVarCube = psArrayAllocEmpty (1); + psArray *sourceMskCube = psArrayAllocEmpty (1); + + psImage *refImage = NULL; +# endif + // build (DATA - MODEL) [an image] for each psf star psArray *input = psArrayAllocEmpty (100); @@ -108,4 +120,24 @@ 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); + if (refImage == NULL) { + refImage = tmpImage; + psArrayAdd (sourceRawCube, 1, tmpImage); + } else { + if (tmpImage->numCols != refImage->numCols) { + fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + continue; + } + if (tmpImage->numRows != refImage->numRows) { + fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + continue; + } + psArrayAdd (sourceRawCube, 1, tmpImage); + psFree (tmpImage); + } +# endif + pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal); @@ -114,4 +146,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? @@ -142,4 +180,32 @@ 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 @@ -148,7 +214,4 @@ psVector *fmasks = psVectorAlloc (input->n, PS_TYPE_VECTOR_MASK); - // statistic to use to determine baseline for clipping - psStats *fluxClip = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psStats *fluxClipDef = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // statistic to use to determine output flux // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV? @@ -156,4 +219,8 @@ psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV); + // use this for the IRLS fitting below (defines niter) + psStats *statsIRLS = psStatsAlloc(PS_STAT_CLIPPED_MEAN); statsIRLS->clipIter = 10; // max number of iterations + psPolynomial2D *polyIRLS = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1); + // Use psF64 to minimize overflow problems? psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix @@ -166,4 +233,11 @@ int nGoodPixel = 0; // pixel is off the image + + // skip pixels outside of the requested radius + float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); + if (radius > radiusMax) { + resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; + continue; + } // build the vector of data values for this output pixel @@ -194,5 +268,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 @@ -207,15 +291,21 @@ // measure the robust median to determine a baseline reference value - *fluxClip = *fluxClipDef; + psStats *fluxClip = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) { psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats"); + psFree (fluxClip); return false; } - if (isnan(fluxClip->robustMedian)) { + if (isnan(fluxClip->sampleMedian)) { resid->Ro->data.F32[oy][ox] = 0.0; resid->Rx->data.F32[oy][ox] = 0.0; resid->Ry->data.F32[oy][ox] = 0.0; resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask; + psFree (fluxClip); continue; + } + + if (fabs(fluxClip->sampleMedian) > 1.5) { + fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->sampleMedian); } @@ -223,5 +313,5 @@ int nKeep = 0; for (int i = 0; i < fluxes->n; i++) { - float delta = fluxes->data.F32[i] - fluxClip->robustMedian; + float delta = fluxes->data.F32[i] - fluxClip->sampleMedian; float sigma = sqrt (dfluxes->data.F32[i]); float swing = fabs(delta) / sigma; @@ -233,4 +323,9 @@ if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++; } + psFree (fluxClip); + + 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,17 +337,9 @@ } - float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); - if (radius > radiusMax) { - resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; - continue; - } + resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption); - resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; - if (isnan(resid->Ro->data.F32[oy][ox])) { resid->Ro->data.F32[oy][ox] = 0.0; - resid->Rx->data.F32[oy][ox] = 0.0; - resid->Ry->data.F32[oy][ox] = 0.0; resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask; continue; @@ -265,52 +352,71 @@ assert (SPATIAL_ORDER == 1); - float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); - if (radius > radiusMax) { - resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; - continue; - } - - psImageInit(A, 0.0); - psVectorInit(B, 0.0); - for (int i = 0; i < fluxes->n; i++) { - if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i]; - B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i]; - B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i]; - - A->data.F64[0][0] += 1.0/dfluxes->data.F32[i]; - A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i]; - A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i]; - - A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i]; - A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i]; - A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i]; - } - - A->data.F64[0][1] = A->data.F64[1][0]; - A->data.F64[0][2] = A->data.F64[2][0]; - A->data.F64[2][1] = A->data.F64[1][2]; - - if (!psMatrixGJSolve(A, B)) { - resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; - psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox); - continue; - } - - resid->Ro->data.F32[oy][ox] = B->data.F64[0]; - resid->Rx->data.F32[oy][ox] = B->data.F64[1]; - resid->Ry->data.F32[oy][ox] = B->data.F64[2]; - - float dRo = sqrt(A->data.F32[0][0]); + // we have the following vectors to describe this pixel: + // fluxes, dfluxes, fmasks, xC, yC + + psVectorIRLSFitPolynomial2D (polyIRLS, statsIRLS, fmasks, fmaskVal, fluxes, dfluxes, xC, yC); + + resid->Ro->data.F32[oy][ox] = polyIRLS->coeff[0][0]; + resid->Rx->data.F32[oy][ox] = polyIRLS->coeff[1][0]; + resid->Ry->data.F32[oy][ox] = polyIRLS->coeff[0][1]; + + // is the error well-defined? + float dRo = polyIRLS->coeffErr[0][0]; + + if (fabs(resid->Ro->data.F32[oy][ox]) > 1.5) { + fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, resid->Ro->data.F32[oy][ox]); + } if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) { + // fprintf (stderr, "mask : %d, %d : %f +/- %f (%f for %d) : %d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; resid->Ro->data.F32[oy][ox] = 0.0; resid->Rx->data.F32[oy][ox] = 0.0; resid->Ry->data.F32[oy][ox] = 0.0; - } + } else { + // fprintf (stderr, "keep : %d, %d : %f +/- %f (%f for %d) : %d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); + } } } } + +# 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); @@ -329,6 +435,4 @@ psFree (fluxStats); psFree (fluxStatsDef); - psFree (fluxClip); - psFree (fluxClipDef); if (resid != NULL && psTraceGetLevel("psphot") > 5) { @@ -340,5 +444,24 @@ } + psFree (statsIRLS); + psFree (polyIRLS); + psf->residuals = resid; 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; +}