Index: /trunk/psModules/src/pmSubtractBias.c =================================================================== --- /trunk/psModules/src/pmSubtractBias.c (revision 2847) +++ /trunk/psModules/src/pmSubtractBias.c (revision 2848) @@ -6,6 +6,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-29 21:18:33 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-29 22:44:33 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -205,5 +205,5 @@ for (i=0;in) / ((float) n); - newVec->data.F32[i] = psSpline1DEval(x, mySpline); + newVec->data.F32[i] = psSpline1DEval(mySpline, x); } } else { @@ -234,4 +234,6 @@ return(in); } + + // Check for an unallowable pmFit. if ((fit != PM_OVERSCAN_NONE) && (fit != PM_OVERSCAN_ROWS) && @@ -241,4 +243,5 @@ return(in); } + // Check for an unallowable pmOverscanAxis. if ((overScanAxis != PM_OVERSCAN_NONE) && (overScanAxis != PM_OVERSCAN_ROWS) && @@ -259,6 +262,6 @@ double statValue; psImage *myOverscanImage = NULL; - psPolynomial1D *myPoly; - psSpline1D *mySpline; + psPolynomial1D *myPoly = NULL; + psSpline1D *mySpline = NULL; int nBin; // XXX: This comment isn't right? @@ -308,4 +311,9 @@ } + // This check is redundant with above code. + if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) { + psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d)\n", overScanAxis); + return(in); + } tmpOverscan = (psListElem *) overscans->head; @@ -313,126 +321,194 @@ myOverscanImage = (psImage *) tmpOverscan->data; - // XXX: Is there a better way to extract a psVector from a psImage without - // having to copy every element in that vector? - // XXX: If we get here, do we know that overScanAxis == PM_OVERSCAN_ROWS - // or PM_OVERSCAN_COL? - - if ((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS)) { - if (overScanAxis == PM_OVERSCAN_ROWS) { - if (myOverscanImage->numCols != (in->image)->numCols) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n", - myOverscanImage->numCols, in->image->numCols); - } - // We create a row vector and subtract this vector from image. - overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); - for (i=0;in;i++) { - overscanVector->data.F32[i] = 0.0; - } - tmpRow = psVectorAlloc((in->image)->numRows, PS_TYPE_F32); - - for (i=0;inumCols;i++) { - for (j=0;jnumRows;j++) { - tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i]; - } - stat = psVectorStats(stat, tmpRow, NULL, NULL, 0); - p_psGetStatValue(stat, &statValue); - overscanVector->data.F32[i] = statValue; - } - psFree(tmpRow); + if (overScanAxis == PM_OVERSCAN_ROWS) { + if (myOverscanImage->numCols != (in->image)->numCols) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n", + myOverscanImage->numCols, in->image->numCols); + } + + // We create a row vector and subtract this vector from image. + // XXX: Is there a better way to extract a psVector from a psImage without + // having to copy every element in that vector? + overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32); + for (i=0;in;i++) { + overscanVector->data.F32[i] = 0.0; + } + tmpRow = psVectorAlloc((in->image)->numRows, PS_TYPE_F32); + + // For each column of the input image, loop through every row, + // collect the pixel in that row, then performed the specified + // statistical op on those pixels. Store this in overscanVector. + for (i=0;inumCols;i++) { + for (j=0;jnumRows;j++) { + tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i]; + } + stat = psVectorStats(stat, tmpRow, NULL, NULL, 0); + p_psGetStatValue(stat, &statValue); + overscanVector->data.F32[i] = statValue; + } + psFree(tmpRow); + + // Scale the overscan vector to the size of the input image. + if (overscanVector->n != in->image->numCols) { if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { psVector *newVec = ScaleOverscanVector(overscanVector, - (in->image)->numCols, + in->image->numCols, fitSpec, fit); psFree(overscanVector); overscanVector = newVec; - } - } - - if (overScanAxis == PM_OVERSCAN_COLUMNS) { - // We create a column vector and subtract this vector from image. - if (myOverscanImage->numRows != (in->image)->numRows) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n", - myOverscanImage->numRows, in->image->numRows); - } - overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); - for (i=0;in;i++) { - overscanVector->data.F32[i] = 0.0; - } - tmpCol = psVectorAlloc((in->image)->numCols, PS_TYPE_F32); - - for (i=0;inumRows;i++) { - for (j=0;jnumCols;j++) { - tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j]; - } - stat = psVectorStats(stat, tmpCol, NULL, NULL, 0); - p_psGetStatValue(stat, &statValue); - overscanVector->data.F32[i] = statValue; - } - - psFree(tmpCol); + } else { + psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.\n"); + psFree(overscanVector); + return(in); + } + } + } + + if (overScanAxis == PM_OVERSCAN_COLUMNS) { + if (myOverscanImage->numRows != (in->image)->numRows) { + psLogMsg(__func__, PS_LOG_WARN, + "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n", + myOverscanImage->numRows, in->image->numRows); + } + + // We create a column vector and subtract this vector from image. + overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32); + for (i=0;in;i++) { + overscanVector->data.F32[i] = 0.0; + } + tmpCol = psVectorAlloc((in->image)->numCols, PS_TYPE_F32); + + // For each row of the input image, loop through every column, + // collect the pixel in that row, then performed the specified + // statistical op on those pixels. Store this in overscanVector. + for (i=0;inumRows;i++) { + for (j=0;jnumCols;j++) { + tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j]; + } + stat = psVectorStats(stat, tmpCol, NULL, NULL, 0); + p_psGetStatValue(stat, &statValue); + overscanVector->data.F32[i] = statValue; + } + psFree(tmpCol); + + // Scale the overscan vector to the size of the input image. + if (overscanVector->n != in->image->numRows) { if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { psVector *newVec = ScaleOverscanVector(overscanVector, - (in->image)->numRows, + in->image->numRows, fitSpec, fit); psFree(overscanVector); overscanVector = newVec; - } - } - - // - // Re-bin the overscan vector (change its length). - // - if ((nBinOrig > 0) && (nBinOrig < overscanVector->n)) { - numBins = 1+((overscanVector->n)/nBinOrig); - myBin = psVectorAlloc(numBins, PS_TYPE_F32); - binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); - - for (i=0;in > ((i*nBinOrig)+j)) { - binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; - } else { - // XXX: we get here if nBinOrig does not evenly divide - // the overscanVector vector. This is the last bin. Should - // we change the binVec->n to acknowledge that? - binVec->n = j; - } - } - stat = psVectorStats(stat, binVec, NULL, NULL, 0); - p_psGetStatValue(stat, &statValue); - myBin->data.F32[i] = statValue; - } - - // Change the effective size of overscanVector. - overscanVector->n = numBins; - for (i=0;idata.F32[i] = myBin->data.F32[i]; - } - psFree(binVec); - psFree(myBin); - nBin = nBinOrig; - } else { - nBin = 1; - } - - // - // Fit a polynomial to the overscan vector. - // - if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) { - if (fit == PM_FIT_POLYNOMIAL) { - myPoly = (psPolynomial1D *) fitSpec; - myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, NULL); - for (i=0;idata.F32[i] = psPolynomial1DEval(myPoly, (float) i); - } - - } else if (fit == PM_FIT_SPLINE) { - mySpline = (psSpline1D *) fitSpec; - // XXX: The spline functions do not work if the number of - // data points and the number of splines are not equal. - } - } + } else { + psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector. Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.\n"); + psFree(overscanVector); + return(in); + } + } + } + + // + // Re-bin the overscan vector (change its length). + // + if ((nBinOrig > 0) && (nBinOrig < overscanVector->n)) { + numBins = 1+((overscanVector->n)/nBinOrig); + myBin = psVectorAlloc(numBins, PS_TYPE_F32); + binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); + + for (i=0;in > ((i*nBinOrig)+j)) { + binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; + } else { + // XXX: we get here if nBinOrig does not evenly divide + // the overscanVector vector. This is the last bin. Should + // we change the binVec->n to acknowledge that? + binVec->n = j; + } + } + stat = psVectorStats(stat, binVec, NULL, NULL, 0); + p_psGetStatValue(stat, &statValue); + myBin->data.F32[i] = statValue; + } + + // Change the effective size of overscanVector. + overscanVector->n = numBins; + for (i=0;idata.F32[i] = myBin->data.F32[i]; + } + psFree(binVec); + psFree(myBin); + nBin = nBinOrig; + } else { + nBin = 1; + } + + // At this point the number of data points in overscanVector should be + // equal to the number of rows/columns (whatever is appropriate) in the + // image divided by numBins. + // + if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) { + // + // Fit a polynomial or spline to the overscan vector. + // + if (fit == PM_FIT_POLYNOMIAL) { + myPoly = (psPolynomial1D *) fitSpec; + myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, NULL); + if (myPoly == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.\n"); + psFree(overscanVector); + return(in); + } + } else if (fit == PM_FIT_SPLINE) { + mySpline = (psSpline1D *) fitSpec; + mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); + if (mySpline == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.\n"); + psFree(overscanVector); + return(in); + } + } + + // + // Subtract fitted overscan vector row-wise from the image. + // + if (overScanAxis == PM_OVERSCAN_ROWS) { + for (i=0;i<(in->image)->numCols;i++) { + psF32 tmpF32 = 0.0; + if (fit == PM_FIT_POLYNOMIAL) { + tmpF32 = psPolynomial1DEval(myPoly, ((float) i) / ((float) nBin)); + } else if (fit == PM_FIT_SPLINE) { + tmpF32 = psSpline1DEval(mySpline, ((float) i) / ((float) nBin)); + } + + for (j=0;j<(in->image)->numRows;j++) { + (in->image)->data.F32[j][i]-= tmpF32; + + } + } + } + + // + // Subtract fitted overscan vector column-wise from the image. + // + if (overScanAxis == PM_OVERSCAN_COLUMNS) { + for (i=0;i<(in->image)->numRows;i++) { + psF32 tmpF32 = 0.0; + if (fit == PM_FIT_POLYNOMIAL) { + tmpF32 = psPolynomial1DEval(myPoly, ((float) i) / ((float) nBin)); + } else if (fit == PM_FIT_SPLINE) { + tmpF32 = psSpline1DEval(mySpline, ((float) i) / ((float) nBin)); + } + + for (j=0;j<(in->image)->numCols;j++) { + (in->image)->data.F32[i][j]-= tmpF32; + } + } + } + } else { + // If we get here, then no polynomials were fit to the overscan + // vector. We simply subtract it, taking into account binning, + // from the image. // @@ -457,9 +533,10 @@ } } - psFree(overscanVector); - } + } + + psFree(overscanVector); + tmpOverscan = tmpOverscan->next; } - psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, Index: /trunk/psModules/src/pmSubtractSky.c =================================================================== --- /trunk/psModules/src/pmSubtractSky.c (revision 2847) +++ /trunk/psModules/src/pmSubtractSky.c (revision 2848) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-29 21:18:33 $ + * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-29 22:44:33 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -423,4 +423,5 @@ psTrace(".psModule.pmSubtractSky", 4, "---- pmSubtractSky() begin ----\n"); + if ((fit != PM_FIT_NONE) && (fit != PM_FIT_POLYNOMIAL) &&