Index: /trunk/ppBackground/src/ppBackgroundStack.h =================================================================== --- /trunk/ppBackground/src/ppBackgroundStack.h (revision 36634) +++ /trunk/ppBackground/src/ppBackgroundStack.h (revision 36635) @@ -23,8 +23,15 @@ psImageMap *modelMap; + // These are the full extent of the input data psF32 ra_min; psF32 ra_max; psF32 dec_min; psF32 dec_max; + + // Because the binning code can't handle subsections sanely, these hold the working subsection. + psF32 x_min; + psF32 x_max; + psF32 y_min; + psF32 y_max; psArray *stack_data; Index: /trunk/ppBackground/src/ppBackgroundStackCamera.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36634) +++ /trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36635) @@ -33,5 +33,5 @@ pmConfig *config = data->config; // Because I'm reusing code. int u,v; - size_t A,P; + // size_t A,P; double RR = -9999.0 ,DD = -99999.0,rr = 9999.0,dd = 99999.0; // FIX Figure out what stacks we have to deal with. @@ -115,6 +115,6 @@ } - psMemStats(0,&A,&P); - fprintf(stderr,"fpa %ld %ld\n",A,P); +/* psMemStats(0,&A,&P); */ +/* fprintf(stderr,"fpa %ld %ld\n",A,P); */ pmChip *chip; // Chip from FPA @@ -132,6 +132,6 @@ } - psMemStats(0,&A,&P); - fprintf(stderr,"chip %ld %ld\n",A,P); +/* psMemStats(0,&A,&P); */ +/* fprintf(stderr,"chip %ld %ld\n",A,P); */ // read WCS data from the corresponding header @@ -279,11 +279,6 @@ - psMemStats(0,&A,&P); - fprintf(stderr,"chip2 %ld %ld\n",A,P); - -/* if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { */ -/* psError(psErrorCodeLast(), false, "Error saving data to files."); */ -/* return false; */ -/* } */ +/* psMemStats(0,&A,&P); */ +/* fprintf(stderr,"chip2 %ld %ld\n",A,P); */ } // end chip loop @@ -306,5 +301,5 @@ pmFPAfileActivate(config->files, false, NULL); psFree(modelContent); - // psFree(smfFile); + // psFree(smfFile); // This gets freed when we trash the config->files structure. I think. psFree(smfFileName); psFree(view); @@ -313,5 +308,5 @@ config->files = psMetadataAlloc(); // - + // } // end smf loop @@ -321,18 +316,36 @@ config->files = config_files; - // Allocate the modelMap for the region we're covering. - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = 33; // Number of samples - binning->nYruff = 33; - binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at - binning->nYfine = ceil(data->dec_max - data->dec_min) + 1; - binning->nXskip = floor(data->ra_min); - binning->nYskip = floor(data->dec_min); - psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER); - printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n", - rr,dd,RR,DD,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine); - data->modelMap = psImageMapNoImageAlloc( binning,stats); + // remove the {ra|dec}_min values from everything so we don't confuse the binning code. +/* psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL); */ +/* psMetadataItem *expItem; */ +/* while ((expItem = psMetadataGetAndIncrement(expIter))) { */ +/* if (expItem->type != PS_DATA_METADATA) { */ +/* continue; // This is the N counter */ +/* } */ +/* psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL); */ +/* psMetadataItem *chipItem; */ +/* while ((chipItem = psMetadataGetAndIncrement(chipIter))) { */ +/* psImage *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); */ +/* psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); */ +/* for (v = 0; v < ra->numRows; v++) { */ +/* for (u = 0; u < ra->numCols; u++) { */ +/* ra->data.F32[v][u] -= data->ra_min; */ +/* dec->data.F32[v][u] -= data->dec_min; */ +/* } */ +/* } */ +/* } // End chip */ +/* psFree(chipIter); */ +/* } // End smfs */ +/* psFree(expIter); */ +/* // And from the structure objects. */ +/* data->x_min -= data->ra_min; */ +/* data->y_min -= data->dec_min; */ +/* data->x_max -= data->ra_min; */ +/* data->y_max -= data->dec_min; */ +/* /\* data->ra_max -= data->ra_min; *\/ */ +/* /\* data->dec_max -= data->dec_min; *\/ */ +/* /\* data->ra_min = 0.0; *\/ */ +/* /\* data->dec_min = 0.0; *\/ */ + return true; Index: /trunk/ppBackground/src/ppBackgroundStackLoop.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36634) +++ /trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36635) @@ -68,6 +68,90 @@ } - - + } // End initialization block. + + + // Loop over the input images, and apply the models to construct the restored versions. + for (i = 0; i < data->stack_data->n; i++) { + pmFPAfile *stack = data->stack_data->data[i]; + pmFPAview *view = pmFPAviewAlloc(0); + + // PART 1: + // Determine the extent of the model map for this stack + data->x_min = 99e99; data->x_max = -99e99; + data->y_min = 99e99; data->y_max = -99e99; + // Allocate the modelMap for the region we're covering. + psPlane *pix = psPlaneAlloc(); // Pixel coordinates on chip + psPlane *tp = psPlaneAlloc(); // Focal plane coordinates + + pix->x = 0; pix->y = 0; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + pix->x = 6240; pix->y = 0; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + pix->x = 6240; pix->y = 6243; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + pix->x = 0; pix->y = 6243; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + +/* data->x_min -= data->ra_min; */ +/* data->x_max -= data->ra_min; */ +/* data->y_min -= data->dec_min; */ +/* data->y_max -= data->dec_min; */ + + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = 15; // Number of samples + binning->nYruff = 15; + // binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at + // binning->nYfine = ceil(data->dec_max - data->dec_min) + 1; + binning->nXfine = ceil(data->x_max - data->x_min) + 1; + binning->nYfine = ceil(data->y_max - data->y_min) + 1; + binning->nXskip = floor(data->x_min - data->ra_min) + 1; + binning->nYskip = floor(data->y_min - data->dec_min) + 1; + psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER); + printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n", + 0.0,0.0,0.0,0.0,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine); + printf("Sizes: corners: %f %f -> %f %f map: %d %d\n", + data->x_min,data->y_min,data->x_max,data->y_max, + binning->nXfine,binning->nYfine); + psImage *sizeImage = psImageAlloc(binning->nXfine,binning->nYfine,PS_TYPE_F32); + P_PSIMAGE_SET_COL0(sizeImage, data->x_min); + P_PSIMAGE_SET_ROW0(sizeImage, data->y_min); + data->modelMap = psImageMapAlloc(sizeImage,binning,stats); +/* psFree(sizeImage); */ +/* data->modelMap = psImageMapNoImageAlloc( binning,stats); */ +/* P_PSIMAGE_SET_COL0(data->modelMap->map, (data->x_min - binning->nXskip) / binning->nXbin); */ +/* P_PSIMAGE_SET_ROW0(data->modelMap->map, (data->y_min - binning->nYskip) / binning->nYbin); */ + + // force col0/row0 +/* data->modelMap->map->col0 = data->modelMap->binning->nXskip; */ +/* data->modelMap->map->row0 = data->modelMap->binning->nYskip; */ + + // PART 2: + // Solve the data into a model for this region of the sky. + // This seems wrong, but I need a blank modelMap object, so we fit the zero-data we've stored in the calib objects printf("Determining blank modelMap!\n"); @@ -83,41 +167,36 @@ return(false); } - } // End initialization block. - - // This is where a loop would likely start. - { - // Construct the offset information - printf("Model fit!\n"); - if (!ppBackgroundStackDataModelFit(data)) { - psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling."); - return(false); - } - - // Apply full correction - printf("Calib apply!\n"); - if (!ppBackgroundStackCalibApply(data)) { - psError(psErrorCodeLast(), false, "Error applying the calibration models."); - return(false); - } - - // Fit the new model - printf("Determining model!\n"); - if (!ppBackgroundStackModelFit(data)) { - psError(psErrorCodeLast(), false, "Error determining the modelMap object."); - return(false); - } - } // End loop - - // Loop over the input images, and apply the models to construct the restored versions. - - for (i = 0; i < data->stack_data->n; i++) { - pmFPAfile *stack = data->stack_data->data[i]; - pmFPAview *view = pmFPAviewAlloc(0); - + + // This is where an iterative solution loop would likely start. + { + // Construct the offset information + printf("Model fit!\n"); + if (!ppBackgroundStackDataModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling."); + return(false); + } + + // Apply full correction + printf("Calib apply!\n"); + if (!ppBackgroundStackCalibApply(data)) { + psError(psErrorCodeLast(), false, "Error applying the calibration models."); + return(false); + } + + // Fit the new model + printf("Determining model!\n"); + if (!ppBackgroundStackModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the modelMap object."); + return(false); + } + } // End loop + + // PART 3: + // Define output products // Define output image. Why is this always so hard to do? pmFPA *tmp_fpa1,*tmp_fpa2; tmp_fpa1 = pmFPAConstruct(config->camera,config->cameraName); tmp_fpa2 = pmFPAConstruct(config->camera,config->cameraName); - + pmFPAfile *stack_model = pmFPAfileDefineOutput(config,tmp_fpa1,"PPBACKGROUND.STACK.MODEL"); @@ -134,5 +213,5 @@ stack_model->save = true; stack_corr->save = true; - + printf("I'm about to loop over the parts of this stack: %d\n",i); // Iterate over the images. @@ -144,5 +223,4 @@ return(false); } - pmChip *chip; @@ -151,5 +229,5 @@ continue; } - + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { psError(psErrorCodeLast(), false, "load failure for Chip"); @@ -256,5 +334,6 @@ } psFree(view); - } + psFree(data->modelMap); + } Index: /trunk/ppBackground/src/ppBackgroundStackMath.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackMath.c (revision 36634) +++ /trunk/ppBackground/src/ppBackgroundStackMath.c (revision 36635) @@ -94,11 +94,15 @@ psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); - psVector *obs = psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32); - psVector *model= psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32); + psVector *obs = psVectorAllocEmpty(image->numRows, PS_TYPE_F32); + psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32); for (v = 0; v < image->numRows; v++) { for (u = 0; u < image->numCols; u++) { - obs->data.F32[v * image->numCols + u] = image->data.F32[v][u]; - model->data.F32[v * image->numCols + u] = psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]); + if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)|| + (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; } + psVectorAppend(obs,image->data.F32[v][u]); + psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u])); + // obs->data.F32[v * image->numCols + u] = image->data.F32[v][u]; + // model->data.F32[v * image->numCols + u] = psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]); } } @@ -140,4 +144,6 @@ psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg calibrated data"); psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName); + psImage *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); + psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset"); @@ -146,4 +152,7 @@ for (v = 0; v < image->numRows; v++) { for (u = 0; u < image->numCols; u++) { + if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)|| + (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; } + model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u]; } @@ -164,10 +173,11 @@ int u,v; + long used; psS16 N = psMetadataLookupS16(NULL, data->models, "N"); - psVector *X = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); - psVector *Y = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); - psVector *Z = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); - psVector *E = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); - psVector *mask = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_VECTOR_MASK); + psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *Y = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *Z = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *E = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *mask = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_VECTOR_MASK); j = 0; @@ -187,12 +197,31 @@ psImage *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); + // psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image"); for (v = 0; v < calib->numRows; v++) { for (u = 0; u < calib->numCols; u++) { - X->data.F32[j] = ra->data.F32[v][u]; - Y->data.F32[j] = dec->data.F32[v][u]; - Z->data.F32[j] = calib->data.F32[v][u]; - E->data.F32[j] = 1.0; - mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0; + if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)|| + (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { + j++; + continue; } + psVectorAppend(X,ra->data.F32[v][u]); + psVectorAppend(Y,dec->data.F32[v][u]); + psVectorAppend(Z,calib->data.F32[v][u]); + psVectorAppend(E,1.0); + psVectorAppend(mask,0); +/* X->data.F32[j] = ra->data.F32[v][u]; */ +/* Y->data.F32[j] = dec->data.F32[v][u]; */ +/* Z->data.F32[j] = calib->data.F32[v][u]; */ +/* E->data.F32[j] = 1.0; */ +/* mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0; */ + used++; j++; + +/* printf("DATA %ld %ld %f %f %f %f\n", */ +/* j,used, */ +/* ra->data.F32[v][u], */ +/* dec->data.F32[v][u], */ +/* calib->data.F32[v][u], */ +/* model->data.F32[v][u]); */ + } } @@ -200,7 +229,7 @@ psFree(chipIter); } // End smfs - + printf("%ld %ld\n", j,used); bool fitStatus; - bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 0, X, Y, Z, E); + bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E); psFree(expIter);