Index: /branches/czw_branch/20120906/ippconfig/recipes/nightly_science.config =================================================================== --- /branches/czw_branch/20120906/ippconfig/recipes/nightly_science.config (revision 34441) +++ /branches/czw_branch/20120906/ippconfig/recipes/nightly_science.config (revision 34442) @@ -223,4 +223,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG # WARP S16 60 @@ -234,4 +235,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG END @@ -244,4 +246,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG END @@ -254,4 +257,5 @@ DIFFABLE BOOL TRUE REDUCTION STR SWEETSPOT + ADDITIONAL_STACK_LABEL STR ecliptic.rp # REDUCTION STR SSTF_4SIG END Index: /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c =================================================================== --- /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c (revision 34441) +++ /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c (revision 34442) @@ -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,71 @@ +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 + // const psImage *mask = interp->mask; // Mask to interpolate + const psImage *variance = interp->variance; // Variance 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 (varianceValue && variance) { + u1 = (xh - x) * variance->data.F32[yl][xl] + (x - xl) * variance->data.F32[yl][xh]; + u2 = (xh - x) * variance->data.F32[yh][xl] + (x - xl) * variance->data.F32[yh][xh]; + + *varianceValue = (yh - y) * u1 + (y - yl) * u2; + } +#if 0 + if (maskValue && mask) { + u1 = (xh - x) * mask->data.F32[yl][xl] + (x - xl) * mask->data.F32[yl][xh]; + u2 = (xh - x) * mask->data.F32[yh][xl] + (x - xl) * mask->data.F32[yh][xh]; + + maskValue = (yh - y) * u1 + (y - yl) * u2; + } +#endif + 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 (*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] + ); + } + return PS_INTERPOLATE_STATUS_GOOD; +} + psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue, @@ -969,7 +1039,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 +1091,5 @@ switch (mode) { case PS_INTERPOLATE_FLAT: + case PS_INTERPOLATE_BILINEAR_SIMPLE: // No smearing by design return 1.0; @@ -1083,4 +1156,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/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.h =================================================================== --- /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.h (revision 34441) +++ /branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.h (revision 34442) @@ -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/czw_branch/20120906/psModules/src/detrend/pmDark.c =================================================================== --- /branches/czw_branch/20120906/psModules/src/detrend/pmDark.c (revision 34441) +++ /branches/czw_branch/20120906/psModules/src/detrend/pmDark.c (revision 34442) @@ -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/czw_branch/20120906/psModules/src/objects/Makefile.am =================================================================== --- /branches/czw_branch/20120906/psModules/src/objects/Makefile.am (revision 34441) +++ /branches/czw_branch/20120906/psModules/src/objects/Makefile.am (revision 34442) @@ -111,4 +111,5 @@ pmSourceOutputs.h \ pmSourceIO.h \ + pmSourceSatstar.h \ pmSourcePlots.h \ pmSourceVisual.h \ Index: /branches/czw_branch/20120906/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- /branches/czw_branch/20120906/psModules/src/objects/pmFootprintCullPeaks.c (revision 34441) +++ /branches/czw_branch/20120906/psModules/src/objects/pmFootprintCullPeaks.c (revision 34442) @@ -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/czw_branch/20120906/pswarp/src/pswarpArguments.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpArguments.c (revision 34441) +++ /branches/czw_branch/20120906/pswarp/src/pswarpArguments.c (revision 34442) @@ -170,15 +170,15 @@ psWarning("BACKGROUND.MODEL is not set in the %s recipe -- defaulting to FALSE.", PSWARP_RECIPE); } - int bkgXgrid = psMetadataLookupBool(&status,recipe, "BKG.XGRID"); ///< Xsize of background model + 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 = psMetadataLookupBool(&status,recipe, "BKG.YGRID"); ///< Xsize of background model + 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) Index: /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (revision 34441) +++ /branches/czw_branch/20120906/pswarp/src/pswarpLoop.c (revision 34442) @@ -159,5 +159,16 @@ } } - + pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header); + + fprintf(stderr,"WCSFI: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", + WCSF->crval1,WCSF->crval2, + WCSF->crpix1,WCSF->crpix2, + WCSF->cdelt1,WCSF->cdelt2, + WCSF->trans->x->coeff[1][0], + WCSF->trans->x->coeff[0][1], + WCSF->trans->y->coeff[1][0], + WCSF->trans->y->coeff[0][1], + -1.0,-1.0); + psFree(WCSF); psList *cells = psListAlloc(NULL); // List of cells, for concepts averaging @@ -197,4 +208,16 @@ } } + pmAstromWCS *WCSF = pmAstromWCSfromHeader(hdu->header); + + fprintf(stderr,"WCSFIB: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", + WCSF->crval1,WCSF->crval2, + WCSF->crpix1,WCSF->crpix2, + WCSF->cdelt1,WCSF->cdelt2, + WCSF->trans->x->coeff[1][0], + WCSF->trans->x->coeff[0][1], + WCSF->trans->y->coeff[1][0], + WCSF->trans->y->coeff[0][1], + -1.0,-1.0); + psFree(WCSF); @@ -574,4 +597,66 @@ bool bilevelAstrometry = false; pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa); + + pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header); + + fprintf(stderr,"WCSF: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", + WCSF->crval1,WCSF->crval2, + WCSF->crpix1,WCSF->crpix2, + WCSF->cdelt1,WCSF->cdelt2, + WCSF->trans->x->coeff[1][0], + WCSF->trans->x->coeff[0][1], + WCSF->trans->y->coeff[1][0], + WCSF->trans->y->coeff[0][1], + -1.0,-1.0); +/* #define CORRECT_INPUT_WCS 1 */ +/* #if CORRECT_INPUT_WCS */ +/* // Correct the input WCS */ +/* pmAstromWCS *WCS = pmAstromWCSfromHeader(phu->header); */ + +/* double cd1f = 1.0 * 400; */ +/* double cd2f = 1.0 * 400; */ +/* fprintf(stderr,"WCS: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", */ +/* WCS->crval1,WCS->crval2, */ +/* WCS->crpix1,WCS->crpix2, */ +/* WCS->cdelt1,WCS->cdelt2, */ +/* WCS->trans->x->coeff[1][0], */ +/* WCS->trans->x->coeff[0][1], */ +/* WCS->trans->y->coeff[1][0], */ +/* WCS->trans->y->coeff[0][1], */ +/* cd1f,cd2f); */ + +/* WCS->crpix1 = WCS->crpix1 / cd1f; */ +/* WCS->crpix2 = WCS->crpix2 / cd2f; */ + +/* WCS->cdelt1 *= cd1f; */ +/* WCS->cdelt2 *= 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] *= cd1f; */ +/* } */ +/* } */ +/* 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] *= cd2f; */ +/* } */ +/* } */ +/* fprintf(stderr,"WCS: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", */ +/* WCS->crval1,WCS->crval2, */ +/* WCS->crpix1,WCS->crpix2, */ +/* WCS->cdelt1,WCS->cdelt2, */ +/* WCS->trans->x->coeff[1][0], */ +/* WCS->trans->x->coeff[0][1], */ +/* WCS->trans->y->coeff[1][0], */ +/* WCS->trans->y->coeff[0][1], */ +/* cd1f,cd2f); */ +/* pmAstromWCStoHeader (phu->header,WCS); */ +/* // End WCS work. */ +/* #endif */ + + + + if (phu) { char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); @@ -610,8 +695,52 @@ // read WCS data from the corresponding header pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa); - psMetadataAddS32(config->arguments,PS_LIST_TAIL, - "INTERPOLATION.MODE",PS_META_REPLACE,"", PS_INTERPOLATE_BILINEAR); ///< Mode - psMetadataAddS32(config->arguments,PS_LIST_TAIL, - "INTERPOLATION.NUM",PS_META_REPLACE,"", 64); ///< Mode + + pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); + + double cd1f = 1.0 * 400; + double cd2f = 1.0 * 400; + fprintf(stderr,"WCSC: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", + WCS->crval1,WCS->crval2, + WCS->crpix1,WCS->crpix2, + WCS->cdelt1,WCS->cdelt2, + WCS->trans->x->coeff[1][0], + WCS->trans->x->coeff[0][1], + WCS->trans->y->coeff[1][0], + WCS->trans->y->coeff[0][1], + cd1f,cd2f); + + + WCS->cdelt1 *= cd1f; + WCS->cdelt2 *= cd2f; + WCS->crpix1 = WCS->crpix1 / cd1f + 46.0 / 400.0; + WCS->crpix2 = WCS->crpix2 / cd2f + 68.0 / 400.0; + + // 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++) { + fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->x->coeff[q][r]); + WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->x->coeff[q][r]); + } + } + for (int q = 0; q < WCS->trans->y->nX; q++) { + for (int r = 0; r < WCS->trans->y->nY; r++) { + fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->y->coeff[q][r]); + WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r); + fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->y->coeff[q][r]); + } + } + fprintf(stderr,"WCSO: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", + WCS->crval1,WCS->crval2, + WCS->crpix1,WCS->crpix2, + WCS->cdelt1,WCS->cdelt2, + WCS->trans->x->coeff[1][0], + WCS->trans->x->coeff[0][1], + WCS->trans->y->coeff[1][0], + WCS->trans->y->coeff[0][1], + cd1f,cd2f); + + fprintf(stderr,"Size: %d %d %d %d\n",WCS->trans->x->nX,WCS->trans->x->nY,WCS->trans->y->nX,WCS->trans->y->nY); + pmAstromWCStoHeader (hdu->header,WCS); if (bilevelAstrometry) { @@ -624,5 +753,5 @@ } else { // we use a default FPA pixel scale of 1.0 - if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) { + if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) { psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); psFree(view); @@ -631,27 +760,6 @@ } } -#define CORRECT_INPUT_WCS 1 -#if CORRECT_INPUT_WCS - // Correct the input WCS - pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header); - - double cd1f = 1.0 * 400; - double cd2f = 1.0 * 400; - - WCS->crpix1 = WCS->crpix1 / cd1f + 0.5; - WCS->crpix2 = WCS->crpix2 / cd2f - 0.5; - - 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); - // End WCS work. -#endif + // Modify structure here. + pmCell *cell; @@ -679,10 +787,10 @@ // Copy the detections from the astrometry carrier to the input, so they can be accessed by // pswarpTransformReadout - pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry - pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry - if (detections) { - psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections); - } - +/* pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry */ +/* pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry */ +/* if (detections) { */ +/* psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY | PS_META_REPLACE, "Sources from input astrometry", detections); */ +/* } */ + for (int x = 0; x < readout->image->numCols; x++) { for (int y = 0; y < readout->image->numRows; y++) { @@ -692,8 +800,8 @@ } } - psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 2); -/* psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NX",PS_META_REPLACE,"",100); ///< Number of grid points in x */ -/* psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NY",PS_META_REPLACE,"",100); ///< Number of grid points in y */ - + + psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8); + + fprintf(stderr,"Transforming Readout!\n"); pswarpTransformReadout(output, readout, config); Index: /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c (revision 34441) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c (revision 34442) @@ -172,5 +172,6 @@ } else { pswarpTransformTileArgs *args = job->args->data[0]; - // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n); + // fprintf (stderr, "finished job %d,%d, Nargs: %ld (%d %d %d %d) (%d %d %d %d)\n", args->gridX, args->gridY, job->args->n, + // xMin,xMax,yMin,yMax,args->xMin,args->xMax,args->yMin,args->yMax); goodPixels += args->goodPixels; xMin = PS_MIN(args->xMin, xMin); Index: /branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c (revision 34441) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c (revision 34442) @@ -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/czw_branch/20120906/pswarp/src/pswarpTransformTile.c =================================================================== --- /branches/czw_branch/20120906/pswarp/src/pswarpTransformTile.c (revision 34441) +++ /branches/czw_branch/20120906/pswarp/src/pswarpTransformTile.c (revision 34442) @@ -83,6 +83,19 @@ // Iterate over the output image pixels (parent frame) long goodPixels = 0; ///< Number of input pixels landing on the output image + psTrace("pswarp",3,"Size: %d %d\n",xMax,yMax); for (int y = yMin; y < yMax; y++) { for (int x = xMin; x < xMax; x++) { +/* if ((x == 931)||(x == 717)||(x == 924)) { */ +/* /\* if ((y == 1348)||(y == 1249)||(y == 1081)) { *\/ */ +/* /\* fprintf(stderr,"BEcause: %d %d [%d %d %d %d]\n",x,y,yMin,yMax,xMin,xMax); *\/ */ +/* /\* } *\/ */ +/* } */ + if (((x == 931)&&(y == 1348))|| + ((x == 717)&&(y == 1249))|| + ((x == 924)&&(y == 1081))) { + psTrace("pswarp",3,"TT: A\n"); + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + + } // Only transform those pixels requested @@ -90,4 +103,11 @@ continue; } + if (((x == 931)&&(y == 1348))|| + ((x == 717)&&(y == 1249))|| + ((x == 924)&&(y == 1081))) { + psTrace("pswarp",3,"TT: B\n"); + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + + } // pswarpMapApply converts the output coordinate (x,y) to the input coordinate. @@ -95,7 +115,22 @@ double xIn, yIn; // Input pixel coordinates pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5); + if (((x == 931)&&(y == 1348))|| + ((x == 717)&&(y == 1249))|| + ((x == 924)&&(y == 1081))) { + psTrace("pswarp",3,"TT: Ct\n"); + psTrace("pswarp",3,"TT: %d %d %d %d %g %g %d %d\n",xMin,xMax,yMin,yMax,xIn,yIn,inNumCols,inNumRows); + + } + if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) { continue; } + if (((x == 931)&&(y == 1348))|| + ((x == 717)&&(y == 1249))|| + ((x == 924)&&(y == 1081))) { + psTrace("pswarp",3,"TT: C\n"); + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + + } // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates @@ -109,4 +144,8 @@ int xOut = x - outCol0, yOut = y - outRow0; ///< Position on output image + if ((xIn >= args->interp->image->numCols - 1)|| + (yIn >= args->interp->image->numRows - 1)) { + imageValue = NAN; + } if (outImageData) { outImageData[yOut][xOut] = imageValue * jacobian; @@ -120,4 +159,43 @@ goodPixels++; + + if ((x == 931)||(x == 717)||(x == 924)) { + if ((y == 1348)||(y == 1249)||(y == 1081)) { + psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld %g %g %g %g\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, + args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp, + 0.0,0.0,0.0,0.0); +/* inImage->data.F32[(int) yIn][(int) xIn], */ +/* inImage->data.F32[(int) yIn+1][(int) xIn], */ +/* inImage->data.F32[(int) yIn+1][(int) xIn+1], */ +/* inImage->data.F32[(int) yIn][(int) xIn+1]); */ + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + } + } + + if (((xOut == 931)&&(yOut == 1348))|| + ((xOut == 717)&&(yOut == 1249))|| + ((xOut == 924)&&(yOut == 1081))) { + psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld %g %g %g %g\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, + args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp, + 0.0,0.0,0.0,0.0); +/* inImage->data.F32[(int) yIn][(int) xIn], */ +/* inImage->data.F32[(int) yIn+1][(int) xIn], */ +/* inImage->data.F32[(int) yIn+1][(int) xIn+1], */ +/* inImage->data.F32[(int) yIn][(int) xIn+1]); */ + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + + } + if ((xOut == 443)&&(yOut == 659)) { + psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, + args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp); + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + + } + if ((xOut == 524)&&(yOut == 576)) { + psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian, + args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp); + psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax); + + } } }