Index: branches/eam_branches/ipp-20111122/psModules/src/camera/pmReadoutFake.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/camera/pmReadoutFake.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/camera/pmReadoutFake.c (revision 33638) @@ -210,5 +210,5 @@ const pmPSF *psf = args->data[7]; // PSF float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux - float radius = PS_SCALAR_VALUE(args->data[9], F32); // Minimum radius + float radius = PS_SCALAR_VALUE(args->data[9], S32); // Minimum radius - typecast to float from S32 outside of PS_SCALAR_VALUE otherwise sets 0.0 bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF? bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak? @@ -314,5 +314,5 @@ } } - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); psFree(groups); Index: branches/eam_branches/ipp-20111122/psModules/src/concepts/pmConceptsStandard.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/concepts/pmConceptsStandard.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/concepts/pmConceptsStandard.c (revision 33638) @@ -751,13 +751,17 @@ bool has_video_cell = false; - if (concept->type != PS_DATA_STRING) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type for %s (%x) is not string\n", - concept->name, concept->type); - return NULL; - } - - char *Vptr = strchr(concept->data.V,'V'); - if (Vptr) { - has_video_cell = true; + if (concept->type == PS_DATA_BOOL) { + has_video_cell = concept->data.B; + } else { + if (concept->type != PS_DATA_STRING) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Type for %s (%x) is not string\n", + concept->name, concept->type); + return NULL; + } + + char *Vptr = strchr(concept->data.V,'V'); + if (Vptr) { + has_video_cell = true; + } } Index: branches/eam_branches/ipp-20111122/psModules/src/detrend/pmBias.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/detrend/pmBias.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/detrend/pmBias.c (revision 33638) @@ -154,5 +154,5 @@ if (threaded) { // wait here for the threaded jobs to finish - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction."); return false; Index: branches/eam_branches/ipp-20111122/psModules/src/detrend/pmDark.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/detrend/pmDark.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/detrend/pmDark.c (revision 33638) @@ -601,5 +601,5 @@ if (threaded) { // wait here for the threaded jobs to finish - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to apply dark."); psFree(orders); Index: branches/eam_branches/ipp-20111122/psModules/src/detrend/pmFlatField.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/detrend/pmFlatField.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/detrend/pmFlatField.c (revision 33638) @@ -161,5 +161,5 @@ if (threaded) { // wait here for the threaded jobs to finish - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to flat-field image."); return false; Index: branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.c (revision 33638) @@ -6,4 +6,6 @@ #include "pmPattern.h" + +#define PATTERN_ROW_BKG_FIX 1 @@ -89,4 +91,22 @@ psImageInit(corr, NAN); +#ifdef PATTERN_ROW_BKG_FIX + // CZW: 2011-11-30 + // Define the vectors to hold the "x" and "y" slope trends. + // Briefly, the slope trend in the y-axis is a due to variations in the 0-th order term + // of the PATTERN.ROW fit between individual rows across the cell. Similarly, the 1-st + // order term of the PATTERN.ROW fit defines the trend in the x-axis (as that's what we + // are fitting with PATTERN.ROW in the first place). However, the thing we're trying to + // fix with PATTERN.ROW is the detector level bias wiggles. These should be overlaid on + // the true sky level. Therefore, simply applying the PATTERN.ROW correction will + // introduce cell-to-cell sky variations as these two trends are removed. To avoid this, + // We store the 0th and 1st order values used for each row, and then fit a polynomial to + // these results. By re-adding these systematic trends back, we can remove the row-to-row + // variations without improperly removing the real sky trend. + psVector *yaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the constant term + psVector *yaxisMask = psVectorAlloc(numRows, PS_TYPE_VECTOR_MASK); // Mask for rows with no fit + psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term + psVectorInit(yaxisMask, 0); +#endif for (int y = 0; y < numRows; y++) { psVectorInit(clipMask, 0); @@ -105,4 +125,8 @@ // Not enough points to fit patternMaskRow(ro, y, maskBad); +#ifdef PATTERN_ROW_BKG_FIX + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; +#endif continue; } @@ -111,8 +135,22 @@ psErrorClear(); patternMaskRow(ro, y, maskBad); - continue; - } - - poly->coeff[0] -= background; +#ifdef PATTERN_ROW_BKG_FIX + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; +#endif + continue; + } +#ifndef PATTERN_ROW_BKG_FIX + poly->coeff[0] -= background; +#else + // Store the results we found for this row. + yaxisData->data.F32[y] = poly->coeff[0]; + xaxisData->data.F32[y] = poly->coeff[1]; + psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]); + + // yaxisData->data.F32[y] = 0.0; +/* xaxisData->data.F32[y] = 0.0; */ + +#endif memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector @@ -121,4 +159,7 @@ psErrorClear(); patternMaskRow(ro, y, maskBad); +#ifdef PATTERN_ROW_BKG_FIX + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; +#endif continue; } @@ -126,8 +167,94 @@ for (int x = 0; x < numCols; x++) { image->data.F32[y][x] -= solution->data.F32[x]; + psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]); } psFree(solution); } +#ifdef PATTERN_ROW_BKG_FIX + // Put the global trends back that were removed by the PATTERN.ROW correction. + // Set up the indices for the polynomial + psVector *yaxisIndices = psVectorAlloc(numRows, PS_TYPE_F32); + norm = 2.0 / (float)numRows; + for (int y = 0; y < numRows; y++) { + yaxisIndices->data.F32[y] = y * norm - 1.0; + psTrace("psModules.detrend.pattern",10,"%d %f %f\n",y,yaxisIndices->data.F32[y],yaxisData->data.F32[y]); + } + + // Fit the trend of the constant term, producing the y-axis global trend + psStatsInit(clip); + psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit. + if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) { + psWarning("Unable to fit polynomial to y-axis trend"); + psErrorClear(); + // If we've failed, we need to do something, so add back in the background level, and + // expect that the final image will have background mismatches. + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += background; + corr->data.F64[y][0] -= background; + } + } + } + else { + psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices); + if (!solution) { + psWarning("Unable to evaluate polynomial"); + psErrorClear(); + // If we've failed, we need to do something, so add back in the background level, and + // expect that the final image will have background mismatches. + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += background; + corr->data.F64[y][0] -= background; + } + } + } + else { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += solution->data.F32[y]; + corr->data.F64[y][0] -= solution->data.F32[y]; + psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]); + } + } + } + psFree(solution); + } + + // Fit the trend of the linear term, producing the x-axis global trend + // We can use the same mask vector, as the same rows failed the row-fit earlier. + psStatsInit(clip); + psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit. + if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) { + psWarning("Unable to fit polynomial to x-axis trend"); + psErrorClear(); + } + else { + psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices); + if (!solution) { + psWarning("Unable to evaluate polynomial"); + psErrorClear(); + } + else { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x]; + corr->data.F64[y][1] -= solution->data.F32[y] ; + psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],indices->data.F32[x]); + } + } + } + psFree(solution); + } + psFree(yaxisPoly); + psFree(xaxisPoly); + psFree(yaxisIndices); + psFree(yaxisMask); + psFree(yaxisData); + psFree(xaxisData); + // End PATTERN_ROW_BKG_FIX global trend replacement +#endif + psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE, "Pattern row correction", corr); @@ -382,2 +509,432 @@ + +bool pmPatternContinuity(pmChip *chip, const psVector *tweak, psStatsOptions bgStat, psStatsOptions cellStat, + psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth) +{ + PS_ASSERT_PTR_NON_NULL(chip, false); + PS_ASSERT_VECTOR_NON_NULL(tweak, false); + PS_ASSERT_VECTOR_SIZE(tweak, chip->cells->n, false); + PS_ASSERT_VECTOR_TYPE(tweak, PS_TYPE_U8, false); + + int numCells = tweak->n; // Number of cells + + psVector *meanMask = psVectorAlloc(numCells, PS_TYPE_VECTOR_MASK); // Mask for means + psVectorInit(meanMask, 0); + + // Mask bits + enum { + PM_PATTERN_IGNORE = 0x01, // Ignore this cell + PM_PATTERN_TWEAK = 0x02, // Tweak this cell + PM_PATTERN_ERROR = 0x04, // Error in calculating background + PM_PATTERN_ALL = 0xFF, // All causes + }; + + // Count number of cells to tweak + int numTweak = 0; // Number of cells to tweak + int numIgnore = 0; // Number of cells to ignore + for (int i = 0; i < numCells; i++) { + pmCell *cell = chip->cells->data[i]; // Cell of interest + if (!cell || !cell->data_exists || !cell->process || + cell->readouts->n == 0 || cell->readouts->n > 1 || !cell->readouts->data[0]) { + numIgnore++; + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_IGNORE; + continue; + } + if (tweak->data.U8[i]) { + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_TWEAK; + numTweak++; + } + } + if (numTweak == 0) { + // Nothing to do + psFree(meanMask); + return true; + } + + // Measure mean of each cell edge, and use that to determine the cell offsets. + + psStats *bgStats = psStatsAlloc(bgStat); // Statistics on background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + + psRegion region = {0,0,0,0}; + + /* These images hold the edge data for the OTA structure. */ + psImage *A = psImageAlloc(8,8,PS_TYPE_F64); // Top edge + psImage *B = psImageAlloc(8,8,PS_TYPE_F64); // Bottom edge + psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge + psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge + psImageInit(A,0.0); + psImageInit(B,0.0); + psImageInit(C,0.0); + psImageInit(D,0.0); + + for (int i = 0; i < numCells; i++) { + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) { + continue; + } + pmCell *cell = chip->cells->data[i]; // Cell of interest + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + + psStatsInit(bgStats); + + // Convert cell iterator i into an xy coordinate on the grid of cells + int y = (i % 8); + int x = (i - y) / 8; + + for (int j = 0; j < 4; j++) { + if (j == 0) { // Region B + region = psRegionSet(0,ro->image->numCols, + 0,edgeWidth); + } + else if (j == 1) { // Region A + region = psRegionSet(0,ro->image->numCols, + ro->image->numRows - edgeWidth,ro->image->numRows); + } + else if (j == 2) { // Region D + region = psRegionSet(0,edgeWidth, + 0,ro->image->numRows); + } + else if (j == 3) { // Region C + region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols, + 0,ro->image->numRows); + } + psImage *subset = psImageSubset(ro->image,region); + psImage *submask = psImageSubset(ro->mask,region); + + if (!psImageBackground(bgStats, NULL, subset, submask, maskVal, rng)) { + psWarning("Unable to measure background for cell %d on edge %d\n", i, j); + psErrorClear(); + meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR; + if (j == 0) { B->data.F64[y][x] = NAN; } + else if (j == 1) { A->data.F64[y][x] = NAN; } + else if (j == 2) { C->data.F64[y][x] = NAN; } + else if (j == 3) { D->data.F64[y][x] = NAN; } + psFree(subset); + psFree(submask); + continue; // Move on to next edge, as only part of this cell may be a problem + } + + // If the returned value is zero, assume something is wrong. Do I still need this? + if (psStatsGetValue(bgStats,bgStat) < 1e-6) { + if (j == 0) { B->data.F64[y][x] = NAN; } + else if (j == 1) { A->data.F64[y][x] = NAN; } + else if (j == 2) { C->data.F64[y][x] = NAN; } + else if (j == 3) { D->data.F64[y][x] = NAN; } + } + // If we have an error for this cell/edge, make sure we mask the value + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) { + if (j == 0) { B->data.F64[y][x] = NAN; } + else if (j == 1) { A->data.F64[y][x] = NAN; } + else if (j == 2) { C->data.F64[y][x] = NAN; } + else if (j == 3) { D->data.F64[y][x] = NAN; } + } + else { // Set the value to match what we got from the edge box. + if (j == 0) { B->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); } + else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); } + else if (j == 2) { C->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); } + else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); } + } + + for (int u = 0; u < subset->numCols; u++) { + for (int v = 0; v < subset->numRows; v++) { + psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d", + i,j,x,y,u,v,subset->data.F32[v][u],submask->data.PS_TYPE_IMAGE_MASK_DATA[v][u]); + } + } + + psFree(subset); + psFree(submask); + + } + psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) A: %f B: %f C: %f D: %f", + i,x,y, + A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]); + } + psFree(bgStats); + psFree(rng); + + // We've now allocated all the edge values, so we can now minimize the offsets. + // This involves solving the equation A x = b, where + // A is the (64x64 for GPC1) matrix containing the edges that match for each cell + // x is the solution vector + // b is the combination of offsets across each cell boundary for each cell. + // Below "XX" is used as the matrix A, and "solution" is used as both b and x + // (due to the way psMatrixLUSolve operates). + psVector *solution = psVectorAlloc(64,PS_TYPE_F64); + psImage *XX = psImageAlloc(64,64,PS_TYPE_F64); + psVectorInit(solution,0.0); + psImageInit(XX,0.0); + + for (int i = 0; i < numCells; i++) { + // Accumulate all the possible edge differences we can for this cell. + // As we do so, make a note of the correlations by incrementing the element of the matrix. + int y = (i % 8); + int x = (i - y) / 8; + int j; + double critical_value = 0.0; + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) { + continue; + } + if (x + 1 < 8) { // We have a neighbor adjacent in the +x direction + j = 8 * (x + 1) + y; // Determine that neighbor's index + if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) { + critical_value = 2.0 * fabs(D->data.F64[y][x+1]); + } + else { + critical_value = 2.0 * fabs(C->data.F64[y][x]); + } + if (critical_value < 25) { critical_value = 25; } + psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic + i,x,y,j, + C->data.F64[y][x], + D->data.F64[y][x+1], + critical_value + ); + if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&& // If there are no errors with the neighbor, + (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&& // and all edges have valid values, + (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value) // and there are no large discontinuities, + ) { + solution->data.F64[i] += C->data.F64[y][x] - D->data.F64[y][x+1]; // Take the difference + XX->data.F64[i][i] += 1; // increment our relation with ourself + XX->data.F64[i][j] += -1; // decrement our relation with the neighbor + } + } + if (x - 1 > -1) { // etc. + j = 8 * (x - 1) + y; + if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) { + critical_value = 2.0 * fabs(D->data.F64[y][x]); + } + else { + critical_value = 2.0 * fabs(C->data.F64[y][x-1]); + } + if (critical_value < 25) { critical_value = 25; } + psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g", + i,x,y,j, + D->data.F64[y][x], + C->data.F64[y][x-1], + critical_value + ); + + if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&& + (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&& + (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value) + ) { + solution->data.F64[i] += D->data.F64[y][x] - C->data.F64[y][x-1]; + XX->data.F64[i][i] += 1; + XX->data.F64[i][j] += -1; + } + } + if (y + 1 < 8) { + j = 8 * x + (y + 1); + psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g", + i,x,y,j, + A->data.F64[y][x], + B->data.F64[y+1][x] + ); + if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) { + critical_value = 2.0 * fabs(B->data.F64[y+1][x]); + } + else { + critical_value = 2.0 * fabs(A->data.F64[y][x]); + } + if (critical_value < 25) { critical_value = 25; } + if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&& + (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&& + (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value) + ) { + solution->data.F64[i] += A->data.F64[y][x] - B->data.F64[y+1][x]; + XX->data.F64[i][i] += 1; + XX->data.F64[i][j] += -1; + } + } + if (y - 1 > -1) { + j = 8 * x + (y - 1); + psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g", + i,x,y,j, + B->data.F64[y][x], + A->data.F64[y-1][x] + ); + if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) { + critical_value = 2.0 * fabs(B->data.F64[y][x]); + } + else { + critical_value = 2.0 * fabs(A->data.F64[y-1][x]); + } + if (critical_value < 25) { critical_value = 25; } + if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&& + (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&& + (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value) + ) { + solution->data.F64[i] += B->data.F64[y][x] - A->data.F64[y-1][x]; + XX->data.F64[i][i] += 1; + XX->data.F64[i][j] += -1; + } + } + } + double max_XX = 0; + double solution_V = 0; + int i_peak = -1; + for (int i = 0; i < numCells; i++) { // If any cells have no value of themself, set the matrix to 1.0. + if (XX->data.F64[i][i] == 0.0) { + XX->data.F64[i][i] = 1.0; + } + if (XX->data.F64[i][i] > max_XX) { + max_XX = XX->data.F64[i][i]; + solution_V = solution->data.F64[i]; + i_peak = i; + } + } + psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n", + i_peak,solution_V); + + for (int i = 0; i < numCells; i++) { +/* if (!((XX->data.F64[i][i] == 1.0)&& */ +/* (solution->data.F64[i] == 0.0))) { */ + solution->data.F64[i] -= solution_V; + if (i != i_peak) { + for (int j = 0; j < numCells; j++) { + XX->data.F64[i][j] -= XX->data.F64[i_peak][j]; + } + } +/* } */ + } + for (int i = 0; i < numCells; i++) { + XX->data.F64[i_peak][i] = 0.0; + } + XX->data.F64[i_peak][i_peak] = 1.0; + + +#if (1) + for (int i = 0; i < numCells; i++) { // print matrix A + psTrace("psModules.detrend.cont",5,"A: %3d % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f", + i, + XX->data.F64[i][0], XX->data.F64[i][1], XX->data.F64[i][2], XX->data.F64[i][3], + XX->data.F64[i][4], XX->data.F64[i][5], XX->data.F64[i][6], XX->data.F64[i][7], + XX->data.F64[i][8], XX->data.F64[i][9], XX->data.F64[i][10], XX->data.F64[i][11], + XX->data.F64[i][12], XX->data.F64[i][13], XX->data.F64[i][14], XX->data.F64[i][15], + XX->data.F64[i][16], XX->data.F64[i][17], XX->data.F64[i][18], XX->data.F64[i][19], + XX->data.F64[i][20], XX->data.F64[i][21], XX->data.F64[i][22], XX->data.F64[i][23], + XX->data.F64[i][24], XX->data.F64[i][25], XX->data.F64[i][26], XX->data.F64[i][27], + XX->data.F64[i][28], XX->data.F64[i][29], XX->data.F64[i][30], XX->data.F64[i][31], + XX->data.F64[i][32], XX->data.F64[i][33], XX->data.F64[i][34], XX->data.F64[i][35], + XX->data.F64[i][36], XX->data.F64[i][37], XX->data.F64[i][38], XX->data.F64[i][39], + XX->data.F64[i][40], XX->data.F64[i][41], XX->data.F64[i][42], XX->data.F64[i][43], + XX->data.F64[i][44], XX->data.F64[i][45], XX->data.F64[i][46], XX->data.F64[i][47], + XX->data.F64[i][48], XX->data.F64[i][49], XX->data.F64[i][50], XX->data.F64[i][51], + XX->data.F64[i][52], XX->data.F64[i][53], XX->data.F64[i][54], XX->data.F64[i][55], + XX->data.F64[i][56], XX->data.F64[i][57], XX->data.F64[i][58], XX->data.F64[i][59], + XX->data.F64[i][60], XX->data.F64[i][61], XX->data.F64[i][62], XX->data.F64[i][63] + ); + } + + for (int i = 0; i < numCells; i++) { // print vector b + psTrace("psModules.detrend.cont",5,"b: %d %f", + i, + solution->data.F64[i] + ); + } +#endif + + // Solve the Ax=b equation + // psMatrixLUSolve(XX,solution); + psMatrixGJSolve(XX,solution); +#if (1) + for (int i = 0; i < numCells; i++) { // print vector b + psTrace("psModules.detrend.cont",5,"x: %d %f", + i, + solution->data.F64[i] + ); + } +#endif + + /* old code to remove the minimum solution value from the set, to give a "minimal set of offsets." Mathematically unnecessary. */ +/* double min = 99e99; */ +/* for (int i = 0; i < numCells; i++) { */ +/* if (solution->data.F64[i] < min) { */ +/* min = solution->data.F64[i]; */ +/* } */ +/* psTrace("psModules.detrend.cont",5,"x: %d %f %f ", */ +/* i, */ +/* solution->data.F64[i],min */ +/* ); */ +/* } */ +/* for (int i = 0; i < numCells; i++) { */ +/* if (solution->data.F64[i] != 0.0) { */ +/* solution->data.F64[i] -= min; */ +/* } */ +/* } */ + + // Cleanup + psFree(XX); + psFree(A); + psFree(B); + psFree(C); + psFree(D); + + // Correct cells based on the offsets calculated, and store the result in the analysis metadata. + for (int i = 0; i < numCells; i++) { + if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) { + continue; + } + if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_TWEAK)) { + continue; + } + pmCell *cell = chip->cells->data[i]; // Cell of interest + pmReadout *ro = cell->readouts->data[0]; // Readout of interest + + float correction = solution->data.F64[i]; + const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell + psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of cell %s by %f", + cellName, correction); + psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32)); + psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE, + "Pattern cell correction solution", correction); + } + + psFree(solution); + psFree(meanMask); + + return true; +} + +bool pmPatternContinuityApply(pmReadout *ro, psImageMaskType maskBad) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + + bool mdok; // Status of MD lookup + float corr = psMetadataLookupF32(&mdok, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Correction to apply + if (!mdok) { + // No correction to apply + return true; + } + + psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest + int numCols = image->numCols, numRows = image->numRows; // Size of image + + if (!isfinite(corr)) { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] = NAN; + } + } + if (mask) { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad; + } + } + } + } else { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += corr; + } + } + } + + return true; +} + + Index: branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.h =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.h (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/detrend/pmPattern.h (revision 33638) @@ -54,4 +54,21 @@ ); +/// Fix the background on cells known to be troublesome +bool pmPatternContinuity( + pmChip *chip, ///< Chip to correct + const psVector *tweak, ///< U8 vector indicating whether to tweak the corresponding cell + psStatsOptions bgStat, ///< Statistic to use for background measurement + psStatsOptions cellStat, ///< Statistic to use for combination of cell background measurements + psImageMaskType maskVal, ///< Mask value to use + psImageMaskType maskBad, ///< Mask value to give bad pixels + int edgeWidth ///< Size of box to use + ); + +/// Apply previously measured cell pattern correction +bool pmPatternContinuityApply(pmReadout *ro, ///< Readout to correct + psImageMaskType maskBad ///< Mask value to give bad pixels + ); + + /// @} Index: branches/eam_branches/ipp-20111122/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/detrend/pmShutterCorrection.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/detrend/pmShutterCorrection.c (revision 33638) @@ -805,5 +805,5 @@ if (threaded) { // wait here for the threaded jobs to finish - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to apply shutter correction."); psFree(shutterImage); Index: branches/eam_branches/ipp-20111122/psModules/src/extras/psVectorBracket.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/extras/psVectorBracket.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/extras/psVectorBracket.c (revision 33638) @@ -51,6 +51,6 @@ } } - // at this point, index[Nhi] >= key > index[Nlo] - N = Nhi; + N = (Nhi >= index->n) ? Nhi - 1 : Nhi; + // at this point, index[N] >= key > index[Nlo] while ((index->data.F32[N] >= key) && (N > Nlo)) { N--; Index: branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmStackReject.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmStackReject.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmStackReject.c (revision 33638) @@ -313,5 +313,5 @@ } - if (!psThreadPoolWait(false)) { + if (!psThreadPoolWait(false, true)) { psError(psErrorCodeLast(), false, "Unable to grow bad pixels."); psFree(source); Index: branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtraction.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtraction.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtraction.c (revision 33638) @@ -796,5 +796,5 @@ pmSubtractionStamp *stamp = job->args->data[0]; // List of stamps pmSubtractionKernels *kernels = job->args->data[1]; // Kernels - int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index + int footprint = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index -- MEH - it is? return pmSubtractionConvolveStamp(stamp, kernels, footprint); @@ -832,9 +832,11 @@ #ifdef TESTING + //MEH - index conflict or changed in past? for (int j = 0; j < kernels->num; j++) { if (stamp->convolutions1) { psString convName = NULL; - psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j); - psFits *fits = psFitsOpen(convName, "w"); + //psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j); + psStringAppend(&convName, "conv1_xxx_%03d.fits", j); + psFits *fits = psFitsOpen(convName, "w"); psFree(convName); psKernel *conv = stamp->convolutions1->data[j]; @@ -845,6 +847,7 @@ if (stamp->convolutions2) { psString convName = NULL; - psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j); - psFits *fits = psFitsOpen(convName, "w"); + //psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j); + psStringAppend(&convName, "conv2_xxx_%03d.fits", j); + psFits *fits = psFitsOpen(convName, "w"); psFree(convName); psKernel *conv = stamp->convolutions2->data[j]; @@ -905,5 +908,5 @@ } } - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(psErrorCodeLast(), false, "Error waiting for threads."); return false; @@ -1427,5 +1430,5 @@ } - if (!psThreadPoolWait(false)) { + if (!psThreadPoolWait(false, true)) { psError(psErrorCodeLast(), false, "Error waiting for threads."); return false; Index: branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.c (revision 33638) @@ -958,5 +958,5 @@ } - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(psErrorCodeLast(), false, "Error waiting for threads."); return false; Index: branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.v0.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionEquation.v0.c (revision 33638) @@ -882,5 +882,5 @@ } - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(psErrorCodeLast(), false, "Error waiting for threads."); return false; Index: branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionMatch.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionMatch.c (revision 33638) @@ -53,4 +53,5 @@ fprintf(stderr, " Memory in use: %zd\n", totalSize); fprintf(stderr, " Largest block: %ld\n", largest); + //MEH -- osx may not like sbrk fprintf(stderr, " sbrk(): %zd\n", (size_t)sbrk(0)); #endif @@ -122,5 +123,5 @@ PS_ASSERT_FLOAT_LESS_THAN(sysError, 1.0, false); } - if (isfinite(sysError)) { + if (isfinite(skyError)) { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(skyError, 0.0, false); } @@ -1091,5 +1092,5 @@ } - if (!psThreadPoolWait(true)) { + if (!psThreadPoolWait(true, true)) { psError(psErrorCodeLast(), false, "Error waiting for threads."); psFree(models); Index: branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionStamps.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/imcombine/pmSubtractionStamps.c (revision 33638) @@ -1208,5 +1208,13 @@ // XXX this is somewhat arbitrary... if (source->psfMagErr > 0.05) continue; - if (fabs(source->psfMag - source->apMag) > 0.5) continue; + if (isfinite(source->apMag)) { + if (fabs(source->psfMag - source->apMag) > 0.5) continue; + } else if (isfinite(source->apMagRaw)) { + if (fabs(source->psfMag - source->apMagRaw) > 0.5) continue; + } else { + // XXX: Should we carry on or drop this source? + // drop it for now + continue; + } if (source->modelPSF) { Index: branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_QGAUSS.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_QGAUSS.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_QGAUSS.c (revision 33638) @@ -402,4 +402,9 @@ assert (psf->params->n > PM_PAR_YPOS); assert (psf->params->n > PM_PAR_XPOS); + + if (! isfinite(Io)) { + fprintf(stderr, "non-finite Io passed to PM_MODEL_PARAMS_FROM_PSF\n"); + return false; + } PAR[PM_PAR_SKY] = 0.0; Index: branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_SERSIC.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_SERSIC.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/models/pmModel_SERSIC.c (revision 33638) @@ -192,4 +192,10 @@ psF32 z0 = PAR[PM_PAR_I0]*f1; psF32 f0 = PAR[PM_PAR_SKY] + z0; + + if (!isfinite(z0)) { + fprintf(stderr, "z0 is not finite for %f %f %f %f %f. Parameters: \n", X, Y, radius, z, f1); + fprintf(stderr, "%f %f %f %f %f %f %f %f\n", PAR[0], PAR[1], PAR[2], PAR[3], PAR[4], + PAR[5], PAR[6], PAR[7]); + } assert (isfinite(f2)); Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmFootprintCullPeaks.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmFootprintCullPeaks.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmFootprintCullPeaks.c (revision 33638) @@ -178,9 +178,9 @@ psArray *myFP = pmFootprintsFind(subImg, threshold, 5); if (!myFP) { - psWarning ("missing footprint?"); + psWarning ("missing footprint? threshold: %.f", threshold); continue; } if (!myFP->n) { - psWarning ("empty footprint?"); + psWarning ("empty footprint? threshold: %.f", threshold); psFree (myFP); continue; Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmModelUtils.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmModelUtils.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmModelUtils.c (revision 33638) @@ -144,8 +144,16 @@ *Io = source->peak->rawFlux; + +#ifndef ALLOW_NONFINITE_PEAK + // Gene says fail of peak !finite + if (!isfinite(*Io)) return false; +#else + // This is the way it used to be. Somtimes an infinite value Io made it's way down the pipeline + // causing assertion failures if (!isfinite(*Io) && !source->moments) return false; *Io = source->moments->Peak; if (!isfinite(*Io)) return false; +#endif return true; Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmPCMdata.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmPCMdata.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmPCMdata.c (revision 33638) @@ -136,4 +136,18 @@ sum += value; } + } + + if (!(sum > 0.0)) { + // Crazy PSF image print out some debugging information ... + fprintf(stderr, "invalid kernel sum %f found by pmPCMkernelFromPSF\n", sum); for (int j = psf->yMin; j <= psf->yMax; j++) { + fprintf(stderr, "Row %d\n", j); + for (int i = psf->xMin; i <= psf->xMax; i++) { + double value = source->psfImage->data.F32[y0 + j][x0 + i]; + fprintf(stderr, " %d %f\n", i, value); + } + } + fflush(stderr); + // ... but avoid the asssertion two lines down by escaping + goto escape; } assert (sum > 0.0); Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmPSFtryFitEXT.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmPSFtryFitEXT.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmPSFtryFitEXT.c (revision 33638) @@ -73,4 +73,12 @@ continue; } + // If mask object does not exist, mark the source as bad. + // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes + // psphot to exit with a fault. + if (source->maskObj == NULL) { + psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y); + psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; + continue; + } source->modelEXT = pmSourceModelGuess (source, options->type); @@ -89,5 +97,5 @@ // clear object mask to define valid pixels - psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask + psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask // exclude the poor fits Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.c (revision 33638) @@ -57,4 +57,7 @@ #define BLANK_HEADERS "BLANK.HEADERS" // Name of metadata in camera configuration containing header names // for putting values into a blank PHU +static bool pmReadoutReadXSRC(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXFIT(pmFPAfile *file, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); +static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char * exttype, psMetadata *hduHeader, psString xfitname, psArray *sources, long *sourceIndex); // lookup the EXTNAME values used for table data and image header segments @@ -961,5 +964,27 @@ psString dataname = NULL; psString deteffname = NULL; - if (!pmSourceIOextnames(&headname, &dataname, &deteffname, NULL, NULL, NULL, file, view)) { + psString xsrcname = NULL; + psString xfitname = NULL; + psString xradname = NULL; + + // determine the output table format. Assume if we need to output extendend source + // parameters that they may exist in the input. + // XXX: Perhaps we should use different recipe values. + // I.E. EXTENDED_SOURCE_ANALYSIS_READ or something like that + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, "PSPHOT"); + if (!status) { + psError(PS_ERR_UNKNOWN, true, "missing recipe PSPHOT in config data"); + return false; + } + // if this is not TRUE, the output files only contain the psf measurements. + bool XSRC_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_ANALYSIS"); + bool XFIT_OUTPUT = psMetadataLookupBool(&status, recipe, "EXTENDED_SOURCE_FITS"); + bool XRAD_OUTPUT = psMetadataLookupBool(&status, recipe, "RADIAL_APERTURES"); + + if (!pmSourceIOextnames(&headname, &dataname, &deteffname, + XSRC_OUTPUT ? &xsrcname : NULL, + XFIT_OUTPUT ? &xfitname : NULL, + XRAD_OUTPUT ? &xradname : NULL, + file, view)) { return false; } @@ -1039,4 +1064,50 @@ } + long *sourceIndex = NULL; + if (XSRC_OUTPUT || XFIT_OUTPUT || XRAD_OUTPUT) { + long seq_max = -1; + for (long i = sources->n -1; i >= 0; i--) { + pmSource *source = sources->data[i]; + if (source->seq < 0) { + // This can happen cmf files that have been corrupted + psError(PS_ERR_IO, true, "seq < 0 for source %ld: Suspect %s is corrupt", i, file->origname); + return false; + } + if (source->seq > seq_max) { + seq_max = source->seq; + } + } + sourceIndex = psAlloc((seq_max + 1) * sizeof(long)); + for (long i = 0; i < seq_max; i++) { + sourceIndex[i] = -1; + } + for (long i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + sourceIndex[source->seq] = i; + } + } + if (XSRC_OUTPUT && xsrcname) { + if (!pmReadoutReadXSRC(file, exttype, hdu->header, xsrcname, sources, sourceIndex)) { + // XXX: is this an error? + psErrorClear(); + } + psFree(xsrcname); + } + if (XFIT_OUTPUT && xfitname) { + if (!pmReadoutReadXFIT(file, exttype, hdu->header, xfitname, sources, sourceIndex)) { + // XXX: is this an error? + psErrorClear(); + } + psFree(xfitname); + } + if (XRAD_OUTPUT && xradname) { + if (!pmReadoutReadXRAD(file, readout, exttype, hdu->header, xradname, sources, sourceIndex)) { + // XXX: is this an error? + psErrorClear(); + } + psFree(xradname); + } + psFree(sourceIndex); + if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) { #if 0 @@ -1165,3 +1236,86 @@ } - +// XXX: We might be able to macroize this and reuse for the other types + +static bool pmReadoutReadXSRC(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString xsrcname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtName (file->fits, xsrcname)) { + psError(PS_ERR_UNKNOWN, false, "cannot find xsrc extension %s in %s", xsrcname, file->filename); + return false; + } + + psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", xsrcname); + return false; + } + + // XXX these are case-sensitive since the EXTYPE is case-sensitive + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { + if (!strcmp (exttype, "PS1_SV1")) { + status = pmSourcesRead_CMF_PS1_SV1_XSRC (file->fits, hduHeader, sources, sourceIndex); + } + } + psFree(tableHeader); + return status; +} + +static bool pmReadoutReadXFIT(pmFPAfile *file, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtName (file->fits, extname)) { + psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename); + return false; + } + + psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", extname); + return false; + } + + // XXX these are case-sensitive since the EXTYPE is case-sensitive + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { + if (!strcmp (exttype, "PS1_SV1")) { + status = pmSourcesRead_CMF_PS1_SV1_XFIT (file->fits, hduHeader, sources, sourceIndex); + } + } + psFree(tableHeader); + return status; +} +static bool pmReadoutReadXRAD(pmFPAfile *file, pmReadout *readout, char *exttype, psMetadata *hduHeader, psString extname, psArray *sources, long *sourceIndex) +{ + if (!psFitsMoveExtName (file->fits, extname)) { + psError(PS_ERR_UNKNOWN, false, "cannot find extension %s in %s", extname, file->filename); + return false; + } + + psMetadata *tableHeader = psFitsReadHeader(NULL, file->fits); // The FITS header + if (!tableHeader) psAbort("cannot read table header"); + + char *xtension = psMetadataLookupStr (NULL, tableHeader, "XTENSION"); + if (!xtension) psAbort("cannot read table type"); + if (strcmp (xtension, "BINTABLE")) { + psWarning ("no binary table in extension %s, skipping\n", extname); + return false; + } + + // XXX these are case-sensitive since the EXTYPE is case-sensitive + bool status = false; + if (file->type == PM_FPA_FILE_CMF) { + if (!strcmp (exttype, "PS1_SV1")) { + status = pmSourcesRead_CMF_PS1_SV1_XRAD (file->fits, readout, hduHeader, sources, sourceIndex); + } + } + psFree(tableHeader); + return status; +} Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.h =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.h (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO.h (revision 33638) @@ -93,4 +93,7 @@ psArray *pmSourcesRead_CMF_PS1_V4 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_SV1 (psFits *fits, psMetadata *header); +bool pmSourcesRead_CMF_PS1_SV1_XSRC (psFits *fits, psMetadata *header, psArray *sources, long *); +bool pmSourcesRead_CMF_PS1_SV1_XFIT (psFits *fits, psMetadata *header, psArray *sources, long *); +bool pmSourcesRead_CMF_PS1_SV1_XRAD (psFits *fits, pmReadout *readout, psMetadata *header, psArray *sources, long *); psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header); psArray *pmSourcesRead_CMF_PS1_DV2 (psFits *fits, psMetadata *header); Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourceIO_CMF_PS1_SV1.c (revision 33638) @@ -539,4 +539,98 @@ } +bool pmSourcesRead_CMF_PS1_SV1_XSRC(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XSRC Table contains no entries\n"); + return false; + } + + // petrosian mags are not saved, we need to calculate fluxes. For this we need exptime and zero point + float zeropt = psMetadataLookupF32(&status, hduHeader, "FPA.ZP"); + float exptime = psMetadataLookupF32(&status, hduHeader, "EXPTIME"); + float magOffset = zeropt + 2.5*log10(exptime); + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + pmSource *source = NULL; +#ifndef ASSUME_SORTED + long j = seq < sources->n ? seq : sources->n - 1; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } +#else + long j = sourceIndex[seq]; + psAssert(j >= 0 && j < sources->n, "invalid sourceIndex"); + source = sources->data[j]; +#endif + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + + if (!source->extpars) { + source->extpars = pmSourceExtendedParsAlloc (); + } + pmSourceExtendedPars *extpars = source->extpars; + + // Assume that X_EXT Y_EXT and sigmas match the psf src so skip + + // We don't have enough information to calculate the major and minor axis. Set major to 1. Should we scale this by + // psf size or something? + extpars->axes.major = 1.0; + extpars->axes.minor = extpars->axes.major * psMetadataLookupF32(&status, row, "F25_ARATIO"); + extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA"); + + float mag = psMetadataLookupF32(&status, row, "PETRO_MAG"); + float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR"); + if (isfinite(mag)) { + extpars->petrosianFlux = pow(10., (magOffset - mag) / 2.5); + if (isfinite(magErr)) { + extpars->petrosianFluxErr = extpars->petrosianFlux / magErr; + } + } + + extpars->petrosianRadius = psMetadataLookupF32(&status, row, "PETRO_RADIUS"); + extpars->petrosianRadiusErr= psMetadataLookupF32(&status, row, "PETRO_RADIUS_ERR"); + extpars->petrosianR50 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50"); + extpars->petrosianR50Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_50_ERR"); + extpars->petrosianR90 = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90"); + extpars->petrosianR90Err = psMetadataLookupF32(&status, row, "PETRO_RADIUS_90_ERR"); + extpars->petrosianFill = psMetadataLookupF32(&status, row, "PETRO_FILL"); + + psVector *radSB = psMetadataLookupVector(&status, row, "PROF_SB"); + psVector *radFlux = psMetadataLookupVector(&status, row, "PROF_FLUX"); + psVector *radFill = psMetadataLookupVector(&status, row, "PROF_FILL"); + + if (radSB && radSB->n > 0) { + extpars->radProfile = pmSourceRadialProfileAlloc(); + extpars->radProfile->binSB = psMemIncrRefCounter(radSB); + extpars->radProfile->binSum = psMemIncrRefCounter(radFlux); + extpars->radProfile->binFill = psMemIncrRefCounter(radFill); + } + + psFree(row); + } + + return true; +} + // XXX this layout is still the same as PS1_DEV_1 bool pmSourcesWrite_CMF_PS1_SV1_XFIT(psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, char *extname) @@ -684,4 +778,104 @@ psFree (outhead); psFree (table); + return true; +} + +bool pmSourcesRead_CMF_PS1_SV1_XFIT(psFits *fits, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XFIT Table contains no entries\n"); + return false; + } + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (!source->modelFits) { + // XXX: where to find the number of models to expect? + source->modelFits = psArrayAllocEmpty(5); + } + psString modelName = psMetadataLookupStr(&status, row, "MODEL_TYPE"); + if (!modelName) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model name for row %ld\n", i); + psFree(row); + return false; + } + pmModelType modelType = pmModelClassGetType(modelName); + if (modelType < 0) { + psError(PS_ERR_UNKNOWN, true, "Failed to find model type for %s\n", modelName); + psFree(row); + return false; + } + pmModel *model = pmModelAlloc(modelType); + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + PAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT"); + PAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT"); + dPAR[PM_PAR_XPOS] = psMetadataLookupF32(&status, row, "X_EXT_SIG"); + dPAR[PM_PAR_YPOS] = psMetadataLookupF32(&status, row, "Y_EXT_SIG"); + + model->mag = psMetadataLookupF32(&status, row, "EXT_INST_MAG"); + model->magErr = psMetadataLookupF32(&status, row, "EXT_INST_MAG_SIG"); + + psEllipseAxes axes; + axes.major = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ"); + axes.minor = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN"); + axes.theta = psMetadataLookupF32(&status, row, "EXT_THETA"); + if (!pmPSF_AxesToModel(PAR, axes, modelType)) { + // Do we need to fail here or can this happen? + psError(PS_ERR_UNKNOWN, false, "Failed to convert psf axes to model"); + psFree(model); + psFree(row); + return false; + } + // XXX: clean this up + if (model->params->n > 7) { + PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07"); + } + // read the covariance matrix + int nparams = model->params->n; + psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32); + for (int y = 0; y < nparams; y++) { + for (int x = 0; x < nparams; x++) { + char name[64]; + snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x); + covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name); + } + } + model->covar = covar; + + psArrayAdd(source->modelFits, 1, model); + psFree(model); + + psFree(row); + } + return true; } @@ -831,2 +1025,94 @@ return true; } + +bool pmSourcesRead_CMF_PS1_SV1_XRAD(psFits *fits, pmReadout *readout, psMetadata *hduHeader, psArray *sources, long *sourceIndex) +{ + PS_ASSERT_PTR_NON_NULL(fits, false); + PS_ASSERT_PTR_NON_NULL(sources, false); + + bool status; + long numSources = psFitsTableSize(fits); // Number of sources in table + if (numSources == 0) { + psError(psErrorCodeLast(), false, "XRAD Table contains no entries\n"); + return false; + } + + long seq_first = -1; + long seq_last = -1; + psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32); + long max_entries = -1; + long num_entries = -1; + + for (long i = 0; i < numSources; i++) { + psMetadata *row = psFitsReadTableRow(fits, i); // Table row + if (!row) { + psError(psErrorCodeLast(), false, "Unable to read row %ld of sources", i); + psFree(row); + return false; + } + // Find the source with this sequence number. + // XXX: I am assuming that sources is sorted in order of seq. + long seq = psMetadataLookupU32 (&status, row, "IPP_IDET"); + long j = seq < sources->n ? seq : sources->n - 1; + pmSource *source = NULL; + for (; j >= 0; j--) { + source = sources->data[j]; + if (source->seq == seq) { + break; + } + } + if (!source) { + psError(PS_ERR_UNKNOWN, false, "Failed to find source for row %ld sequence number %ld\n", i, seq); + psFree(row); + return false; + } + if (seq_first == -1) { + seq_first = seq; + } + if (seq == seq_first) { + psF32 value = psMetadataLookupF32(&status, row, "PSF_FWHM"); + psVectorAppend(fwhmValues, value); + } + if (seq == seq_last) { + num_entries++; + } else { + num_entries = 1; + seq_last = seq; + } + if (num_entries > max_entries) { + max_entries = num_entries; + } + + if (!source->radialAper) { + // XXX: where to find the number of models to expect? + source->radialAper = psArrayAllocEmpty(5); + } + pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc(); + + radialAper->flux = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX")); + radialAper->fluxStdev = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_STDEV")); + radialAper->fluxErr = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FLUX_ERR")); + radialAper->fill = psMemIncrRefCounter(psMetadataLookupVector(&status, row, "APER_FILL")); + + psArrayAdd(source->radialAper, 1, radialAper); + + psFree(radialAper); + psFree(row); + } + + // check for consistency between the length of fwhmValues and the maximum number of entries for each row + if (fwhmValues->n != max_entries) { + psError(PS_ERR_PROGRAMMING, true, "number of PSF_FWHM values found %ld does not match expected number: %ld\n", + fwhmValues->n, max_entries); + psAssert(0, "fixme"); + } + + if (!readout->analysis) { + readout->analysis = psMetadataAlloc(); + } + + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues); + psFree(fwhmValues); + + return true; +} Index: branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourcePhotometry.c (revision 33096) +++ branches/eam_branches/ipp-20111122/psModules/src/objects/pmSourcePhotometry.c (revision 33638) @@ -569,5 +569,4 @@ } } - return (true); }