Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 5623) +++ /trunk/psLib/src/astro/psCoord.c (revision 5624) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.93 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-23 23:54:43 $ +* @version $Revision: 1.94 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-30 02:00:00 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -61,12 +61,10 @@ should rename this. -To derive this transformation, start with the simple 2-by-2 matrix inversion -based on discriminants. This will invert the - (x_2, y_2) = matrix(a b c d) * vector (x, y) -where there are no constant terms. Then you substitute - x_2 = x_1 - e - y_2 = y_1 - f -for (x_2, y_2) to get the desired inverse. +XXX: Use the ADD version which is based on determinants. *****************************************************************************/ + +// XXX EAM : below is the code using the standard matrix representation. +// note that this inversion requires x->nX == 1, y->nY == 1 and +// x->nY <= 1, y->nX <= 1 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) { @@ -130,4 +128,32 @@ // printf("HMMM: out->y: (%f %f %f)\n", out->y->coeff[0][0], out->y->coeff[1][0], out->y->coeff[0][1]); + + // unless the cross terms are available, set these matrix elements to 0 + psF64 r12 = 0.0; + if (transform->x->nY == 1) { + r12 = transform->x->coeff[0][1]; + } + psF64 r21 = 0.0; + if (transform->y->nX == 1) { + r21 = transform->y->coeff[1][0]; + } + psF64 r11 = transform->x->coeff[1][0]; + psF64 r22 = transform->y->coeff[0][1]; + psF64 xo = transform->x->coeff[0][0]; + psF64 yo = transform->y->coeff[0][0]; + + psF64 invDet = 1.0 / (r11 * r22 - r12 * r21); + + // apply the results back to the polynomials + out->x->coeff[0][0] = -invDet * (r22 * xo - r12 * yo); + out->y->coeff[0][0] = -invDet * (r11 * yo - r21 * xo); + out->x->coeff[1][0] = +invDet * r22; + out->y->coeff[0][1] = +invDet * r11; + if (transform->x->nY == 1) { + out->x->coeff[0][1] = -invDet * r12; + } + if (transform->y->nX == 1) { + out->y->coeff[1][0] = -invDet * r21; + } return(out); } @@ -327,6 +353,5 @@ XXX: Private Function. -piNormalize(): take an input angle in radians and convert it to the range -0:2*PI. +piNormalize(): take an input angle in radians and convert it to the range 0:2*PI. *****************************************************************************/ psF32 piNormalize(psF32 angle) @@ -380,6 +405,9 @@ PS_ASSERT_PTR_NON_NULL(projection, NULL); - psF64 theta = 0.0; - psF64 phi = 0.0; + psF64 phi, theta; + psF64 sinDp, cosDp, sinAlpha, cosAlpha, sinDelta, cosDelta; + psF64 sinTheta, cosPhiCT, sinPhiCT, zeta; + + bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN); // Allocate return value @@ -391,29 +419,43 @@ } - // Convert to projection spherical coordinate system - theta = asin( sin(coord->d)*sin(projection->D) + - cos(coord->d)*cos(projection->D)*cos(coord->r-projection->R)); - phi = atan2( -1.0*cos(coord->d)*sin(coord->r-projection->R), - sin(coord->d)*cos(projection->D) - cos(coord->d)*sin(projection->D)*cos(coord->r-projection->R) ); + if (zenithal) { + sinDp = sin(projection->D); + cosDp = cos(projection->D); + sinAlpha = sin(coord->r-projection->R); + cosAlpha = cos(coord->r-projection->R); + sinDelta = sin(coord->d); + cosDelta = cos(coord->d); + + sinTheta = sinDelta*sinDp + cosDelta*cosDp*cosAlpha; + cosPhiCT = sinDelta*cosDp - cosDelta*sinDp*cosAlpha; + sinPhiCT = -cosDelta*sinAlpha; + } else { + phi = coord->r - projection->R; + theta = coord->d - projection->D; + } // Perform the specified projection - // Gnomonic projection - if (projection->type == PS_PROJ_TAN) { - out->x = (cos(theta)*sin(phi))/sin(theta); - out->y = (-1.0*cos(theta)*cos(phi))/sin(theta); + switch (projection->type) { + case PS_PROJ_TAN: + // Gnomonic projection + out->x = +sinPhiCT / sinTheta; + out->y = -cosPhiCT / sinTheta; + break; + case PS_PROJ_SIN: // Othrographic projection - } else if (projection->type == PS_PROJ_SIN) { - out->x = cos(theta)*sin(phi); - out->y = -1.0*cos(theta)*cos(phi); + out->x = +sinPhiCT; + out->y = -cosPhiCT; + break; + case PS_PROJ_AIT: // Hammer-Aitoff projection - } else if ( projection->type == PS_PROJ_AIT) { - psF64 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); + zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); out->x = 2.0*zeta*cos(theta)*sin(phi/2.0); out->y = zeta*sin(theta); + break; + case PS_PROJ_PAR: // Parabolic projection - } else if ( projection->type == PS_PROJ_PAR) { out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0); out->y = M_PI*sin(theta/3.0); - } else { + default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, @@ -424,6 +466,6 @@ // Apply plate scales - out->x *= projection->Xs; - out->y *= projection->Ys; + out->x /= projection->Xs; + out->y /= projection->Ys; // Return output @@ -444,4 +486,10 @@ PS_ASSERT_PTR_NON_NULL(coord, NULL); PS_ASSERT_PTR_NON_NULL(projection, NULL); + + psF64 rho = 0.0; + psF64 sinTheta = 0.0; + psF64 cosTheta = 0.0; + psF64 sinPhi = 0.0; + psF64 cosPhi = 0.0; psF64 theta = 0.0; @@ -458,30 +506,42 @@ // Remove plate scales - // XXX: Verify this. EAM suggested we do a multiply, however that does - // not make sense if we also do the multiply in psProject(). - psF64 x = coord->x/projection->Xs; - psF64 y = coord->y/projection->Ys; + psF64 x = coord->x*projection->Xs; + psF64 y = coord->y*projection->Ys; + psF64 R = sqrt(x*x + y*y); + + bool zenithal = (projection->type == PS_PROJ_TAN) ||(projection->type == PS_PROJ_SIN); // Perform inverse projection - // Gnonomic deprojection - if ( projection->type == PS_PROJ_TAN) { - phi = atan(-1.0*x/y); - theta = atan(1.0/sqrt(x*x+y*y)); + switch (projection->type) { + case PS_PROJ_TAN: + // Gnonomic deprojection + rho = sqrt (1 + R*R); + sinTheta = 1 / rho; + cosTheta = R / rho; + sinPhi = (R == 0) ? 0.0 : +x / R; + cosPhi = (R == 0) ? 1.0 : -y / R; + break; + case PS_PROJ_SIN: // Orhtographic deprojection - } else if ( projection->type == PS_PROJ_SIN) { - phi = atan((-1.0*x)/y); - theta = atan( sqrt(1.0-(x*x+y*y)) / sqrt(x*x+y*y)); + sinTheta = sqrt (1 - R*R); + cosTheta = R; + sinPhi = (R == 0) ? 0.0 : +x / R; + cosPhi = (R == 0) ? 1.0 : -y / R; + break; + case PS_PROJ_AIT: // Hammer-Aitoff deprojection - } else if ( projection->type == PS_PROJ_AIT) { - psF64 z = sqrt(1.0 - ((x/4.0)*(x/4.0)) - ((y/2.0)*(y/2.0))); - phi = 2.0*atan((z*x) / (2.0*(2.0*z*z-1.0)) ); - theta = asin(y*z); + // XXX EAM : need range check on z^2 : must be > 0 + // XXX EAM : old code, ADD, and elixir code are discrepant re x/4, y/2 + rho = sqrt(1.0 - PS_SQR(x/4.0) - PS_SQR(y/2.0)); + phi = 2.0*atan2((2.0*rho*rho-1.0), x*rho); + theta = asin(y*rho); + break; + case PS_PROJ_PAR: // Parabolic deprojection - } else if ( projection->type == PS_PROJ_PAR) { - psF64 rho = y/M_PI; + rho = y/M_PI; phi = x/(1.0 - 4.0*rho*rho); theta = 3.0*asin(rho); - // Invalid deprojection type - } else { + break; + default: psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psCoord_PROJECTION_TYPE_UNKNOWN, @@ -491,10 +551,19 @@ } - // Convert from projection spherical coordinates - out->d = asin( sin(theta)*sin(projection->D) + - cos(theta)*cos(projection->D)*cos(phi) ); - out->r = projection->R + atan2( -1.0*cos(theta)*sin(phi), - sin(theta)*cos(projection->D) - - cos(theta)*sin(projection->D)*cos(phi) ); + if (zenithal) { + psF64 sinDp = sin(projection->D); + psF64 cosDp = cos(projection->D); + + // Convert from projection spherical coordinates + psF64 delta = asin(sinTheta*sinDp + cosTheta*cosDp*cosPhi); + psF64 sinAlphaF = -cosTheta*sinPhi; + psF64 cosAlphaF = -cosTheta*cosPhi*sinDp + sinTheta*cosDp; + + out->d = delta; + out->r = atan2(sinAlphaF, cosAlphaF) + projection->R; + } else { + out->r = phi + projection->R; + out->d = theta + projection->D; + } // Return sphere coordinate Index: /trunk/psLib/src/astro/psCoord.h =================================================================== --- /trunk/psLib/src/astro/psCoord.h (revision 5623) +++ /trunk/psLib/src/astro/psCoord.h (revision 5624) @@ -1,19 +1,19 @@ /** @file psCoord.h -* -* @brief Contains basic coordinate transformation definitions and operations -* -* This file defines the basic types for astronomical coordinate -* transformation -* -* @ingroup CoordinateTransform -* -* -* @author GLG, MHPCC -* -* @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-18 19:39:29 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -*/ + * + * @brief Contains basic coordinate transformation definitions and operations + * + * This file defines the basic types for astronomical coordinate + * transformation + * + * @ingroup CoordinateTransform + * + * + * @author GLG, MHPCC + * + * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:00:00 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii + */ #ifndef PS_COORD_H @@ -323,4 +323,7 @@ psPlaneTransform *transform ///< transform to invert ); +psPlaneTransform *p_psPlaneTransformLinearInvert_MHPCC( + psPlaneTransform *transform ///< transform to invert +); Index: /trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- /trunk/psLib/src/imageops/psImageConvolve.c (revision 5623) +++ /trunk/psLib/src/imageops/psImageConvolve.c (revision 5624) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-22 20:15:35 $ + * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:00:07 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -487,5 +487,5 @@ } -void psImageSmooth (psImage *image, +bool psImageSmooth (psImage *image, double sigma, double Nsigma) @@ -569,6 +569,8 @@ PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED, typeStr); - } - } + return false; + } + } + return true; } Index: /trunk/psLib/src/imageops/psImageConvolve.h =================================================================== --- /trunk/psLib/src/imageops/psImageConvolve.h (revision 5623) +++ /trunk/psLib/src/imageops/psImageConvolve.h (revision 5624) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-31 02:07:11 $ + * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:00:07 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -141,6 +141,8 @@ * Applies a circularly symmetric Gaussian smoothing first in x and then in y * directions with just a vector. This process is 2N faster than 2D convolutions (in general). + * + * @return bool TRUE if successful, otherwise FALSE */ -void psImageSmooth( +bool psImageSmooth( psImage *image, ///< the image to be smoothed double sigma, ///< the width of the smoothing kernel in pixels Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 5623) +++ /trunk/psLib/src/math/psMinimize.c (revision 5624) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.145 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-16 23:06:19 $ + * @version $Revision: 1.146 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:00:09 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1441,4 +1441,69 @@ } +// here is the definition for BuildSums4D. ifdef'ed away until it is used +// by psPolynomial4DFit.. +# if (0) + /****************************************************************************** + BuildSums4D(sums, x, y, z, t, nXterm, nYterm, nZterm, nTterm). equiv to + BuildSums2D(). The result is returned as a double **** + *****************************************************************************/ + static double ****BuildSums4D( + psF64 ****sums, + psF64 x, + psF64 y, + psF64 z, + psF64 t, + psS32 nXterm, + psS32 nYterm, + psS32 nZterm, + psS32 nTterm) +{ + psS32 nXsum = 0; + psS32 nYsum = 0; + psS32 nZsum = 0; + psS32 nTsum = 0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; + psF64 tSum = 1.0; + + nXsum = 2*nXterm; + nYsum = 2*nYterm; + nZsum = 2*nZterm; + nTsum = 2*nTterm; + if (sums == NULL) { + sums = (psF64 ****) psAlloc (nXsum*sizeof(psF64)); + for (int i = 0; i < nXsum; i++) { + sums[i] = (psF64 ***) psAlloc (nYsum*sizeof(psF64)); + for (int j = 0; j < nYsum; j++) { + sums[i][j] = (psF64 **) psAlloc (nZsum*sizeof(psF64)); + for (int k = 0; k < nZsum; k++) { + sums[i][j][k] = (psF64 *) psAlloc (nTsum*sizeof(psF64)); + } + } + } + } + // careful with this function: there is no size checking and realloc for reuse + + tSum = 1.0; + for (int m = 0; m < nTsum; m++) { + zSum = tSum; + for (int k = 0; k < nZsum; k++) { + ySum = zSum; + for (int j = 0; j < nYsum; j++) { + xSum = ySum; + for (int i = 0; i < nXsum; i++) { + sums[i][j][k][m] = xSum; + xSum *= x; + } + ySum *= y; + } + zSum *= z; + } + tSum *= t; + } + return (sums); +} +# endif /* BuildSums4D */ /****************************************************************************** @@ -1447,4 +1512,6 @@ in the output vector. This routine works on single-precision polynomials with double precision data. + +XXX EAM : this function is now deprecated: psPolynomial2DEvalVector handles F32 and F64 *****************************************************************************/ psVector *Polynomial2DEvalVectorD( @@ -1658,7 +1725,10 @@ // xSums look like: 1, x, x^2, ... x^(2n+1) // Build the B and A data structs. + // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects + // XXX EAM : this function is only valid for data vectors of F64 for (int k = 0; k < f->n; k++) { - if ((mask != NULL) && mask->data.U8[k]) + if ((mask != NULL) && (mask->data.U8[k] && maskValue)) { continue; + } if (x != NULL) { xSums = BuildSums1D(xSums, x->data.F64[k], nTerm); @@ -1670,6 +1740,6 @@ wt = 1.0; } else { - // this should filter fErr == 0 values - wt = 1.0 / PS_SQR(fErr->data.F64[k]); + // this filters fErr == 0 values + wt = (fErr->data.F64[k] == 0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]); } for (int i = 0; i < nTerm; i++) { @@ -1686,29 +1756,11 @@ } - // GaussJordan version - if (0) { - // does the solution in place - psGaussJordan (A, B); - - // the first nTerm entries in B correspond directly to the desired - // polynomial coefficients. this is only true for the 1D case - for (int k = 0; k < nTerm; k++) { - myPoly->coeff[k] = B->data.F64[k]; - } - } else { - // LUD version of the fit - psImage *ALUD = NULL; - psVector* outPerm = NULL; - psVector* coeffs = NULL; - - ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); - ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - for (int k = 0; k < nTerm; k++) { - myPoly->coeff[k] = coeffs->data.F64[k]; - } - psFree(ALUD); - psFree(coeffs); - psFree(outPerm); + // does the solution in place + psGaussJordan (A, B); + + // the first nTerm entries in B correspond directly to the desired + // polynomial coefficients. this is only true for the 1D case + for (int k = 0; k < nTerm; k++) { + myPoly->coeff[k] = B->data.F64[k]; } @@ -1819,5 +1871,5 @@ } - +// This function accepts F32 and F64 input vectors. psPolynomial1D *psVectorClipFitPolynomial1D( psPolynomial1D *poly, @@ -1827,8 +1879,8 @@ const psVector *f, const psVector *fErr, - const psVector *x) + const psVector *xIn) { // Internal pointers for possibly NULL vectors. - psVector *x32 = NULL; + psVector *x = NULL; PS_ASSERT_POLY_NON_NULL(poly, NULL); @@ -1836,24 +1888,93 @@ PS_ASSERT_PTR_NON_NULL(stats, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); if (mask != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); } - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + if (xIn != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, xIn, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(xIn, NULL); } if (fErr != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } - - psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + } + + // assign sequence vector if xIn is NULL + if (xIn == NULL) { + x = psVectorCreate (NULL, 0, f->n, 1, f->type.type); + } else { + x = (psVector *) xIn; + } + + // clipping range defined by min and max and/or clipSigma + float minClipSigma; + float maxClipSigma; + if (isfinite(stats->max)) { + maxClipSigma = fabs(stats->clipSigma); + } else { + maxClipSigma = fabs(stats->max); + } + if (isfinite(stats->min)) { + minClipSigma = fabs(stats->clipSigma); + } else { + minClipSigma = fabs(stats->min); + } + psVector *fit = NULL; + psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); + + // eventual expansion: user supplies one of various stats option pairs, + // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to + // evaluate the clipping sigma + // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used + stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + for (int N = 0; N < stats->clipIter; N++) { + int Nkeep = 0; + + poly = psVectorFitPolynomial1D (poly, mask, maskValue, f, fErr, x); + fit = psPolynomial1DEvalVector (poly, x); + resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); + + stats = psVectorStats (stats, resid, NULL, mask, maskValue); + float minClipValue = -minClipSigma*stats->sampleStdev; + float maxClipValue = +maxClipSigma*stats->sampleStdev; + + // set mask if pts are not valid + // we are masking out any point which is out of range + // recovery is not allowed with this scheme + for (int i = 0; i < resid->n; i++) { + if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { + continue; + } + if (resid->data.F64[i] - stats->sampleMedian > maxClipValue) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + if (resid->data.F64[i] - stats->sampleMedian < minClipValue) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + Nkeep ++; + } + + psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", + Nkeep, x->n); + + psFree (fit); + } // Free psVectors that were created for NULL arguments. - if (x == NULL) { - psFree(x32); - } - return(NULL); + if (xIn == NULL) { + psFree(x); + } + // Free other local temporary variables + psFree (resid); + + return (poly); } @@ -1929,6 +2050,8 @@ // Build the B and A data structs. for (int k = 0; k < x->n; k++) { - if ((mask != NULL) && mask->data.U8[k]) + if ((mask != NULL) && (mask->data.U8[k] & maskValue)) { continue; + } + Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm); @@ -1936,8 +2059,6 @@ wt = 1.0; } else { - // XXX: this should probably by fErr^2 !! - // this should filter fErr == 0 values - // XXX: Why isn't this fErr^2? - wt = 1.0 / fErr->data.F64[k]; + // this filters fErr == 0 values + wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]); } @@ -1962,11 +2083,7 @@ // does the solution in place - // XXX: Check return codes! psGaussJordan (A, B); - // XXX: Check return codes! - // ALUD = psMatrixLUD(ALUD, &outPerm, A); - // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - + // select the appropriate solution entries for (int n = 0; n < nXterm; n++) { for (int m = 0; m < nYterm; m++) { @@ -1985,4 +2102,6 @@ +// XXX EAM : I have implemented a single function to handle the mask/nomask cases +// this function can be deprecated psPolynomial2D* RobustFit2D_nomask( psPolynomial2D* poly, @@ -2251,5 +2370,5 @@ PS_ASSERT_PTR_NON_NULL(stats, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); if (mask != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); @@ -2258,21 +2377,75 @@ PS_ASSERT_VECTOR_NON_NULL(x, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); if (fErr != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } - - if (mask == NULL) { - // XXX: Change argument order. - poly = RobustFit2D_nomask(poly, x, y, f, fErr); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + } + + // clipping range defined by min and max and/or clipSigma + float minClipSigma; + float maxClipSigma; + if (isfinite(stats->max)) { + maxClipSigma = fabs(stats->max); } else { - // XXX: Use maskValue. - // XXX: Change argument order. - poly = RobustFit2D(poly, mask, x, y, f, fErr); - } + maxClipSigma = fabs(stats->clipSigma); + } + if (isfinite(stats->min)) { + minClipSigma = fabs(stats->min); + } else { + minClipSigma = fabs(stats->clipSigma); + } + psVector *fit = NULL; + psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64); + + // eventual expansion: user supplies one of various stats option pairs, + // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to + // evaluate the clipping sigma + // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used + stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + + for (int N = 0; N < stats->clipIter; N++) { + int Nkeep = 0; + + poly = psVectorFitPolynomial2D (poly, mask, maskValue, f, fErr, x, y); + fit = psPolynomial2DEvalVector (poly, x, y); + resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); + + stats = psVectorStats (stats, resid, NULL, mask, maskValue); + float minClipValue = -minClipSigma*stats->sampleStdev; + float maxClipValue = +maxClipSigma*stats->sampleStdev; + + // set mask if pts are not valid + // we are masking out any point which is out of range + // recovery is not allowed with this scheme + for (int i = 0; i < resid->n; i++) { + if ((mask != NULL) && (mask->data.U8[i] & maskValue)) { + continue; + } + if (resid->data.F64[i] - stats->sampleMedian > maxClipValue) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + if (resid->data.F64[i] - stats->sampleMedian < minClipValue) { + if (mask != NULL) { + mask->data.U8[i] |= 0x01; + } + continue; + } + Nkeep ++; + } + + psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n", + Nkeep, x->n); + + psFree (fit); + } + // Free local temporary variables + psFree (resid); if (poly == NULL) { Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 5623) +++ /trunk/psLib/src/math/psPolynomial.c (revision 5624) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.132 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-22 21:40:40 $ +* @version $Revision: 1.133 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-11-30 02:00:09 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -725,10 +725,9 @@ } -psPolynomial4D* psPolynomial4DAlloc( - unsigned int nX, - unsigned int nY, - unsigned int nZ, - unsigned int nT, - psPolynomialType type) +psPolynomial4D* psPolynomial4DAlloc( unsigned int nX, + unsigned int nY, + unsigned int nZ, + unsigned int nT, + psPolynomialType type) { PS_ASSERT_INT_NONNEGATIVE(nX, NULL); @@ -785,26 +784,4 @@ } - - - - - - - - - - - - - - - - - - - - - - psF64 psPolynomial1DEval(const psPolynomial1D* poly, psF64 x) @@ -824,4 +801,5 @@ } +// this function must accept F32 and F64 input x vectors psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly, const psVector *x) @@ -829,13 +807,25 @@ PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); psVector *tmp; - tmp = psVectorAlloc(x->n, PS_TYPE_F64); - for (unsigned int i=0;in;i++) { - tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); - } - + switch (x->type.type) { + case PS_TYPE_F64: + tmp = psVectorAlloc(x->n, PS_TYPE_F64); + for (unsigned int i=0;in;i++) { + tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); + } + break; + case PS_TYPE_F32: + tmp = psVectorAlloc(x->n, PS_TYPE_F32); + for (unsigned int i=0;in;i++) { + tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[i]); + } + break; + default: + psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n"); + return (NULL); + } return(tmp); } @@ -859,4 +849,6 @@ } +// this function must support input data types of F32 and F64 +// all input vectors data types must match (all F32 or all F64) psVector *psPolynomial2DEvalVector(const psPolynomial2D *poly, const psVector *x, @@ -866,7 +858,7 @@ PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); psVector *tmp; @@ -878,12 +870,37 @@ } - // Create output vector to return - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial at the specified points - for (unsigned int i=0; idata.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); - } - + switch (x->type.type) { + case PS_TYPE_F32: + if (y->type.type != x->type.type) { + psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors"); + return (NULL); + } + + // Create output vector to return + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + + // Evaluate the polynomial at the specified points + for (unsigned int i=0; idata.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]); + } + break; + case PS_TYPE_F64: + if (y->type.type != x->type.type) { + psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors"); + return (NULL); + } + + // Create output vector to return + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + + // Evaluate the polynomial at the specified points + for (unsigned int i=0; idata.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); + } + break; + default: + psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n"); + return (NULL); + } // Return output vector return(tmp); Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 5623) +++ /trunk/psLib/src/math/psStats.c (revision 5624) @@ -7,19 +7,19 @@ * on those data structures. * -* @author GLG, MHPCC -* -* XXX: The following stats members are never used, or set in this code. -* stats->robustN50 -* stats->clippedNvalues -* stats->binsize -* -* -* -* -* @version $Revision: 1.154 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-11-22 21:40:40 $ -* -* Copyright 2004 Maui High Performance Computing Center, University of Hawaii -*/ + * @author GLG, MHPCC + * + * XXX: The following stats members are never used, or set in this code. + * stats->robustN50 + * stats->clippedNvalues + * stats->binsize + * + * + * + * + * @version $Revision: 1.155 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:00:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ #include @@ -1131,4 +1131,27 @@ statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN); p_psMemSetPersistent(statsTmp, true); + } else { + // XXX EAM : initialize structure if already allocated + statsTmp->sampleMean = NAN; + statsTmp->sampleMedian = NAN; + statsTmp->sampleStdev = NAN; + statsTmp->sampleUQ = NAN; + statsTmp->sampleLQ = NAN; + statsTmp->robustMean = NAN; + statsTmp->robustMedian = NAN; + statsTmp->robustMode = NAN; + statsTmp->robustStdev = NAN; + statsTmp->robustUQ = NAN; + statsTmp->robustLQ = NAN; + statsTmp->robustN50 = -1; // XXX: This is never used + statsTmp->robustNfit = -1; + statsTmp->clippedMean = NAN; + statsTmp->clippedStdev = NAN; + statsTmp->clippedNvalues = -1; // XXX: This is never used + statsTmp->clipSigma = 3.0; + statsTmp->clipIter = 3; + statsTmp->min = NAN; + statsTmp->max = NAN; + statsTmp->binsize = NAN; // XXX: This is never used } Index: /trunk/psLib/src/mathtypes/psImage.c =================================================================== --- /trunk/psLib/src/mathtypes/psImage.c (revision 5623) +++ /trunk/psLib/src/mathtypes/psImage.c (revision 5624) @@ -9,6 +9,6 @@ * @author Ross Harman, MHPCC * - * @version $Revision: 1.90 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-11-22 19:58:16 $ + * @version $Revision: 1.91 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-11-30 02:00:10 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii Index: /trunk/psLib/test/sys/verified/tst_psMemory.stderr =================================================================== --- /trunk/psLib/test/sys/verified/tst_psMemory.stderr (revision 5623) +++ /trunk/psLib/test/sys/verified/tst_psMemory.stderr (revision 5624) @@ -169,12 +169,4 @@ |E|memCheckTypes (FILE:LINENO) psMemCheckBitSet failed in memCheckType. -|E|psPolynomial4DAlloc (FILE:LINENO) - Error: nX is 0 or less. -|E|psPolynomial4DAlloc (FILE:LINENO) - Error: nX is 0 or less. -|E|psPolynomial2DAlloc (FILE:LINENO) - Error: nX is 0 or less. -|E|psPolynomial2DAlloc (FILE:LINENO) - Error: nX is 0 or less. ---> TESTPOINT PASSED (psMemory{psMemCheckType} | tst_psMemory.c)