Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.dat =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.dat (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.dat (revision 37827) @@ -51,3 +51,4 @@ 2453736.5 33.0000000 41317.0 0.0000000 2454832.5 34.0000000 41317.0 0.0000000 - 2456109.5 35.0000000 41317.0 0.0000000 + 2456109.5 35.0000000 41317.0 0.0000000 + 2457204.5 36.0000000 41317.0 0.0000000 Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.raw =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.raw (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.raw (revision 37827) @@ -38,2 +38,3 @@ 2009 JAN 1 =JD 2454832.5 TAI-UTC= 34.0 S + (MJD - 41317.) X 0.0 S 2012 JUL 1 =JD 2456109.5 TAI-UTC= 35.0 S + (MJD - 41317.) X 0.0 S + 2015 JUL 1 =JD 2457204.5 TAI-UTC= 36.0 S + (MJD - 41317.) X 0.0 S Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.c =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.c (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.c (revision 37827) @@ -18,4 +18,5 @@ #include +#include #include "psError.h" #include "psAbort.h" @@ -35,4 +36,5 @@ #include "psImageStructManip.h" #include "psImageMap.h" +#include "psSparse.h" // #include "psImagePixelInterpolate.h" // #include "psImageUnbin.h" @@ -118,5 +120,5 @@ psImageInit (A, 0.0); psVectorInit (B, 0.0); - + // we are looping over the Nx,Ny image map elements; // the matrix equation contains Nx*Ny rows and columns @@ -262,5 +264,5 @@ int I = n + Nx * m; B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py; - + // insert these values into their corresponding locations in A, B // float Sum = 0.0; @@ -273,4 +275,5 @@ int J = (n + jn) + Nx * (m + jm); A->data.F32[J][I] = sA[jn][jm]; + // fprintf (stderr, "A %d %d (%d %d : %d %d): %f\n", I, J, n, m, n + jn, m + jm, sA[jn][jm]); // Sum += sA[jn][jm]; @@ -282,10 +285,17 @@ } // fprintf (stderr, "Total: %f\n", Total); + + double MaxPivot = 0.0; + for (int i = 0; i < Nx*Ny; i++) { + MaxPivot = PS_MAX(MaxPivot, fabs(A->data.F32[i][i])); + // fprintf (stderr, "piv, max: %f : %f\n", A->data.F32[i][i], MaxPivot); + } // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0 psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8); psVectorInit (Empty, 0); + double MinPivot = 0.025*MaxPivot; for (int i = 0; i < Nx*Ny; i++) { - if (A->data.F32[i][i] == 0.0) { + if (fabs(A->data.F32[i][i]) < MinPivot) { Empty->data.S8[i] = 1; for (int j = 0; j < Nx*Ny; j++) { @@ -320,5 +330,5 @@ return true; } - + // set bad values to NaN for (int i = 0; i < Nx*Ny; i++) { @@ -341,5 +351,4 @@ psFree (B); psFree (Empty); - *pGoodFit = true; return true; @@ -402,4 +411,432 @@ psS32 Nkeep = 0; if (!psImageMapFit(pGoodFit, map, mask, maskValue, x, y, f, df)) { + psError(PS_ERR_UNKNOWN, false, "Could not fit image map.\n"); + psFree(resid); + if (!inMask) psFree (mask); + return false; + } + if (!*pGoodFit) { + psWarning ("bad fit to image map, try something else"); + psFree(resid); + if (!inMask) psFree (mask); + return true; + } + + psVector *fit = psImageMapEvalVector(map, mask, maskValue, x, y); + if (fit == NULL) { + psError(PS_ERR_UNKNOWN, false, "Failure in psImageMapEvalVector().\n"); + psFree(resid); + if (!inMask) psFree (mask); + return false; + } + for (int i = 0 ; i < f->n ; i++) { + resid->data.F32[i] = (f->data.F32[i] - fit->data.F32[i]); + } + + if (!psVectorStats(stats, resid, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "Failure to compute statistics on the resid vector.\n"); + psFree(resid); + psFree(fit); + if (!inMask) psFree (mask); + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.imageops", 5, "Mean is %f\n", meanValue); + psTrace("psLib.imageops", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; + + // 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 (psS32 i = 0; i < resid->n; i++) { + // XXX this prevents recovery of previously masked values + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValue) { + continue; + } + + if ((resid->data.F32[i] - meanValue > maxClipValue) || (resid->data.F32[i] - meanValue < minClipValue)) { + psTrace("psLib.imageops", 6, "Masking element %d : %f vs %f : resid is %f\n", i, f->data.F32[i], fit->data.F32[i], resid->data.F32[i]); + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= 0x01; + continue; + } + Nkeep++; + } + + // We should probably exit this loop if no new elements were masked since the fit won't + // change. + psTrace("psLib.imageops", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; + psFree(fit); + } + + // Free local temporary variables + psFree(resid); + if (!inMask) psFree (mask); + *pGoodFit = true; // XXX probably don't need to set this (set by psImageMapFit) + return true; +} + +// CZW: 2014-10-09 +// Sparse versions of MapFit and MapFitClip that assume the matrices are not filled. +bool psImageMapFitSparse(bool *pGoodFit, psImageMap *map, const psVector *mask, psVectorMaskType maskValue, + const psVector *x, const psVector *y, const psVector *f, const psVector *df) +{ + // XXX Add Asserts + + *pGoodFit = false; + + // dimensions of the output map image + int Nx = map->binning->nXruff; + int Ny = map->binning->nYruff; + + + // no spatial information, just calculate mean & stdev + if ((Nx == 1) && (Ny == 1)) { + psStatsInit(map->stats); + + // the user has supplied one of various stats option pairs, + psStatsOptions mean = psStatsMeanOption(map->stats->options); + psStatsOptions stdev = psStatsStdevOption(map->stats->options); + if (!psStatsSingleOption(mean)) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!psStatsSingleOption(stdev)) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; + } + + // XXX does ROBUST_MEDIAN work with weight? + if (!psVectorStats(map->stats, f, NULL, mask, maskValue)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + + map->map->data.F32[0][0] = psStatsGetValue(map->stats, mean); + map->error->data.F32[0][0] = psStatsGetValue(map->stats, stdev); + if (isfinite(map->map->data.F32[0][0]) && isfinite( map->error->data.F32[0][0])) { + *pGoodFit = true; + } + return true; + } + + if (Nx == 1) { + bool status; + status = psImageMapFit1DinY (pGoodFit, map, mask, maskValue, x, y, f, df); + return status; + } + if (Ny == 1) { + bool status; + status = psImageMapFit1DinX (pGoodFit, map, mask, maskValue, x, y, f, df); + return status; + } + + // set up the redirection table so we can use sA[-1][-1], etc + // XXX psKernel does this for you --- PAP. + float SAm[3][3], *SAv[3], **sA; + + for (int i = 0; i < 3; i++) { + SAv[i] = SAm[i] + 1; + } + sA = SAv + 1; + + // elements of the matrix equation Ax = B; we are solving for the vector x + // psImage *A = psImageAlloc (Nx*Ny, Nx*Ny, PS_TYPE_F32); + // psVector *B = psVectorAlloc (Nx*Ny, PS_TYPE_F32); + + // psImageInit (A, 0.0); + // psVectorInit (B, 0.0); + + // CZW: call to psSparseAlloc + // It should match old A, and each element of that should only touch four others. + psSparse *Asparse = psSparseAlloc(Nx * Ny, 4 * Nx * Ny); + + // we are looping over the Nx,Ny image map elements; + // the matrix equation contains Nx*Ny rows and columns + + for (int n = 0; n < Nx; n++) { + for (int m = 0; m < Ny; m++) { + // define & init summing variables + float rx_rx_ry_ry = 0; + float rx_rx_dy_ry = 0; + float dx_rx_ry_ry = 0; + float dx_rx_dy_ry = 0; + float fi_rx_ry = 0; + float rx_rx_py_py = 0; + float rx_rx_qy_py = 0; + float dx_rx_py_py = 0; + float dx_rx_qy_py = 0; + float fi_rx_py = 0; + float px_px_ry_ry = 0; + float px_px_dy_ry = 0; + float qx_px_ry_ry = 0; + float qx_px_dy_ry = 0; + float fi_px_ry = 0; + float px_px_py_py = 0; + float px_px_qy_py = 0; + float qx_px_py_py = 0; + float qx_px_qy_py = 0; + float fi_px_py = 0; + + // generate the sums for the fitting matrix element I,J + // I = n + nX*m + // J = (n + jn) + nX*(m + jm) + for (int i = 0; i < x->n; i++) { + + if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValue)) continue; + + // base coordinate offset for this point (x,y) relative to this map element (n,m) + float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5); + float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5); + + // edge cases to include: + bool edgeX = false; + edgeX |= ((n == 1) && (dx < -1.0)); + edgeX |= ((n == Nx - 2) && (dx > +1.0)); + + bool edgeY = false; + edgeY |= ((m == 1) && (dy < -1.0)); + edgeY |= ((m == Ny - 2) && (dy > +1.0)); + + // skip points outside of 2x2 grid centered on n,m: + if (!edgeX && (fabs(dx) > 1.0)) continue; + if (!edgeY && (fabs(dy) > 1.0)) continue; + + // related offset values + float rx = 1.0 - dx; + float ry = 1.0 - dy; + float px = 1.0 + dx; + float py = 1.0 + dy; + float qx = -dx; + float qy = -dy; + + // data value & weight for this point + float fi = f->data.F32[i]; + if (!isfinite(fi)) continue; + + float wt = 1.0; + if (df != NULL) { + if (df->data.F32[i] == 0.0) { + wt = 0.0; + } else { + if (!isfinite(df->data.F32[i])) continue; + wt = 1.0 / PS_SQR(df->data.F32[i]); // XXX test for dz == NULL or dz_i = 0 + } + } + + // sum the appropriate elements for the different quadrants + + int Qx = (dx >= 0) ? 1 : 0; + if (n == 0) Qx = 1; + if (n == Nx - 1) Qx = 0; + + int Qy = (dy >= 0) ? 1 : 0; + if (m == 0) Qy = 1; + if (m == Ny - 1) Qy = 0; + + assert (isfinite(fi)); + assert (isfinite(wt)); + assert (isfinite(rx)); + assert (isfinite(ry)); + + // points at offset 1,1 + if ((Qx == 1) && (Qy == 1)) { + rx_rx_ry_ry += rx*rx*ry*ry*wt; + rx_rx_dy_ry += rx*rx*dy*ry*wt; + dx_rx_ry_ry += dx*rx*ry*ry*wt; + dx_rx_dy_ry += dx*rx*dy*ry*wt; + fi_rx_ry += fi*rx*ry*wt; + } + // points at offset 1,0 + if ((Qx == 1) && (Qy == 0)) { + rx_rx_py_py += rx*rx*py*py*wt; + rx_rx_qy_py += rx*rx*qy*py*wt; + dx_rx_py_py += dx*rx*py*py*wt; + dx_rx_qy_py += dx*rx*qy*py*wt; + fi_rx_py += fi*rx*py*wt; + } + // points at offset 0,1 + if ((Qx == 0) && (Qy == 1)) { + px_px_ry_ry += px*px*ry*ry*wt; + px_px_dy_ry += px*px*dy*ry*wt; + qx_px_ry_ry += qx*px*ry*ry*wt; + qx_px_dy_ry += qx*px*dy*ry*wt; + fi_px_ry += fi*px*ry*wt; + } + // points at offset 0,0 + if ((Qx == 0) && (Qy == 0)) { + px_px_py_py += px*px*py*py*wt; + px_px_qy_py += px*px*qy*py*wt; + qx_px_py_py += qx*px*py*py*wt; + qx_px_qy_py += qx*px*qy*py*wt; + fi_px_py += fi*px*py*wt; + } + } + + // the chi-square derivatives have elements of the form g(n+jn,m+jm)*A(jn,jm), + // jn,jm = -1 to +1. Convert the sums above into the correct coefficients + sA[-1][-1] = qx_px_qy_py; + sA[-1][ 0] = qx_px_ry_ry + qx_px_py_py; + sA[-1][+1] = qx_px_dy_ry; + sA[ 0][-1] = rx_rx_qy_py + px_px_qy_py; + sA[ 0][ 0] = rx_rx_ry_ry + px_px_ry_ry + rx_rx_py_py + px_px_py_py; + sA[ 0][+1] = rx_rx_dy_ry + px_px_dy_ry; + sA[+1][-1] = dx_rx_qy_py; + sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py; + sA[+1][+1] = dx_rx_dy_ry; + + // I[ 0][ 0] = index for this n,m element: + int I = n + Nx * m; + // B->data.F32[I] = fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py; + // CZW: call to psSparseVector Element + if (fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py == 0.0) { + psSparseVectorElement(Asparse, I, 1.0); + } + else { + psSparseVectorElement(Asparse, I, fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py); + } + + // printf("ADDING: %d %g \n",I, fi_rx_ry + fi_rx_py + fi_px_ry + fi_px_py); + // insert these values into their corresponding locations in A, B + for (int jn = -1; jn <= +1; jn++) { + if (n + jn < 0) continue; + if (n + jn >= Nx) continue; + for (int jm = -1; jm <= +1; jm++) { + if (m + jm < 0) continue; + if (m + jm >= Ny) continue; + int J = (n + jn) + Nx * (m + jm); + // printf("A: %d %d %g\n",J,I,sA[jn][jm]); + // A->data.F32[J][I] = sA[jn][jm]; + // CZW: call to psSparseMatrixElement + if (J < I) { continue; } + psSparseMatrixElement(Asparse,J,I,sA[jn][jm]); // Ensure J < I? + + } + } + } + } + + // test for empty diagonal elements (unconstained cells), mark, and set pivots to 1.0 + // CZW: I'm not totally sure how to check these in the sparse context. + // Iterate over all ii pairs, and manually check the structure? +#if (0) + psVector *Empty = psVectorAlloc (Nx*Ny, PS_TYPE_S8); + psVectorInit (Empty, 0); + for (int i = 0; i < Nx*Ny; i++) { + if (A->data.F32[i][i] == 0.0) { + Empty->data.S8[i] = 1; + for (int j = 0; j < Nx*Ny; j++) { + A->data.F32[i][j] = 0.0; + A->data.F32[j][i] = 0.0; + } + A->data.F32[i][i] = 1.0; + B->data.F32[i] = 0.0; + } + } +#endif + + // CZW: call to psSparseSolve + psVector *solution = psVectorAlloc(Nx*Ny, PS_TYPE_F32); + psSparseConstraint Constraint; + Constraint.paramDelta = 1e-3; + Constraint.paramMin = -1e5; + Constraint.paramMax = 1e5; + solution = psSparseSolve(solution, Constraint, Asparse, 1000); + if (!solution) { + psFree(solution); + psFree(Asparse); + return(false); + } + +#if (0) + // CZW: This is a continuation of the above information + // set bad values to NaN + for (int i = 0; i < Nx*Ny; i++) { + if (Empty->data.S8[i]) { + B->data.F32[i] = NAN; + A->data.F32[i][i] = 0; + } + } +#endif + + for (int n = 0; n < Nx; n++) { + for (int m = 0; m < Ny; m++) { + int I = n + Nx * m; + map->map->data.F32[m][n] = solution->data.F32[I]; + map->error->data.F32[m][n] = NAN; // sqrt(A->data.F32[I][I]); // CZW: fix this to be a real error. + } + } + + // psFree (A); + // psFree (B); + // psFree (Empty); + // CZW: free things + psFree(solution); + psFree(Asparse); + + *pGoodFit = true; + return true; +} + +// measure residuals on each pass and clip outliers based on stats +bool psImageMapClipFitSparse(bool *pGoodFit, psImageMap *map, psStats *stats, psVector *inMask, psVectorMaskType maskValue, + const psVector *x, const psVector *y, const psVector *f, const psVector *df) +{ + // XXX add in full PS_ASSERTS + psAssert(map, "impossible"); + psAssert(stats, "impossible"); + psAssert(x, "impossible"); + psAssert(y, "impossible"); + psAssert(f, "impossible"); + + *pGoodFit = false; + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + psStatsOptions meanOption = psStatsMeanOption(stats->options); + psStatsOptions stdevOption = psStatsStdevOption(stats->options); + if (!psStatsSingleOption(meanOption)) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!psStatsSingleOption(stdevOption)) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; + } + + // clipping range defined by min and max and/or clipSigma + psF32 minClipSigma; + psF32 maxClipSigma; + if (isfinite(stats->max)) { + maxClipSigma = fabs(stats->max); + } else { + maxClipSigma = fabs(stats->clipSigma); + } + if (isfinite(stats->min)) { + minClipSigma = fabs(stats->min); + } else { + minClipSigma = fabs(stats->clipSigma); + } + + psVector *mask = inMask; + if (!inMask) { + mask = psVectorAlloc (x->n, PS_TYPE_VECTOR_MASK); + psVectorInit (mask, 0); + } + + // vector to store residuals + psVector *resid = psVectorAlloc(f->n, PS_TYPE_F32); + + psTrace("psLib.imageops", 4, "stats->clipIter is %d\n", stats->clipIter); + psTrace("psLib.imageops", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); + + for (psS32 N = 0; N < stats->clipIter; N++) { + psTrace("psLib.imageops", 6, "Loop iteration %d. Calling psImageMapFit()\n", N); + psS32 Nkeep = 0; + if (!psImageMapFitSparse(pGoodFit, map, mask, maskValue, x, y, f, df)) { psError(PS_ERR_UNKNOWN, false, "Could not fit image map.\n"); psFree(resid); Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.h =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.h (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.h (revision 37827) @@ -20,4 +20,28 @@ // fit the image map to a set of points bool psImageMapClipFit(bool *pGoodFit, + psImageMap *map, + psStats *stats, + psVector *mask, // WARNING: Mask is modified! + psVectorMaskType maskValue, + const psVector *x, + const psVector *y, + const psVector *f, + const psVector *df + ); + + +// fit the image map to a set of points, using sparse matrix tools +bool psImageMapFitSparse(bool *pGoodFit, + psImageMap *map, + const psVector *mask, + psVectorMaskType maskValue, // + const psVector *x, + const psVector *y, + const psVector *f, + const psVector *df + ); + +// fit the image map to a set of points using sparse matrix tools +bool psImageMapClipFitSparse(bool *pGoodFit, psImageMap *map, psStats *stats, Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.c =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.c (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.c (revision 37827) @@ -134,4 +134,66 @@ assert (isfinite(axes.major)); assert (isfinite(axes.minor)); + + return axes; +} +// ellipse derotation (sx, sy, sxy) -> (major, minor, theta) +psEllipseAxes psEllipseShapeToAxesWithErrors(psEllipseShape shape, psEllipseShape shapeErrors, double maxAR, psEllipseAxes *pShapeErrors) +{ + psEllipseAxes axes; + psEllipseAxes errors; + + double f1 = 1.0 / PS_SQR(shape.sy) + 1.0 / PS_SQR(shape.sx); + double f2 = 1.0 / PS_SQR(shape.sy) - 1.0 / PS_SQR(shape.sx); + + double f32 = PS_SQR(f2) + 4*PS_SQR(shape.sxy); + double f3 = sqrt(f32); + + double df1 = - 2.0 * shapeErrors.sy / (PS_SQR(shape.sy) * shape.sy) + - 2.0 * shapeErrors.sx / (PS_SQR(shape.sx) * shape.sx); + double df2 = - 2.0 * shapeErrors.sy / (PS_SQR(shape.sy) * shape.sy) + + 2.0 * shapeErrors.sx / (PS_SQR(shape.sx) * shape.sx); + + // df3 = 0.5 * df32 / sqrt(f32) = 0.5 * df32 / f3 + double df3 = 0.5 * (2.0*f2*df2 + 2 * 4 * shape.sxy * shapeErrors.sxy) / f3; + + double f13 = f1 + f3; + axes.minor = sqrt (2.0 / f13); + + // dminor = -0.5 * sqrt (2) * df13 * f13**-3/2 = - 0.5 * df13 * minor / f13 + errors.minor = - 0.5 * (df1 + df3) * axes.minor / f13; + + axes.theta = -0.5 * atan2 (+2.0*shape.sxy, f2); + + // according to wikipedia the derivitive of atan2(y, x) is + // + // dAtan2(y, x) = - y * dx / (x**2 + y**2)) + x * dy / (x**2 + y**2) (for x > 0 and y!=0) + // where + // y = 2 * sxy dy = 2 * dsxy x = f2 so dx = df2 + + // dtheta = -0.5 * dAtan2(y, x) + + errors.theta = -0.5 * ( -2.0 * shape.sxy * df2 + f2 * 2 *shapeErrors.sxy ) / + (PS_SQR(f2) + 4 * PS_SQR(shape.sxy)) ; + + // long, thin objects are likely to have a poorly measured major axis + // the angle and minor axis are likely to be ok. + // restrict the axis ratio + double rAR2 = (f1 - f3) / (f1 + f3); + if (rAR2 < 1.0/PS_SQR(maxAR)) { + axes.major = axes.minor * maxAR; + errors.major = errors.minor * maxAR; + } else { + axes.major = sqrt (2.0 / (f1 - f3)); + + // dmajor = -0.5 * (df1 - df3) * sqrt(2) * (f1 - f3)**-3/2 + // = -0.5 * (df2 - df3) * major / (f1 - f3) + errors.major = -0.5 * axes.major * (df1 - df3) / (f1 - f3); + } + + assert (isfinite(axes.theta)); + assert (isfinite(axes.major)); + assert (isfinite(axes.minor)); + + *pShapeErrors = errors; return axes; Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.h =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.h (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.h (revision 37827) @@ -65,4 +65,11 @@ ); +/// Convert shape to axes representation and compute errors +psEllipseAxes psEllipseShapeToAxesWithErrors(psEllipseShape shape, ///< Shape of ellipse + psEllipseShape errors, ///type) { case PS_POLYNOMIAL_ORD: if ((f64->n < 10000)&&(poly->nX > 1)) { scale = true; - sorted = psVectorSort(NULL,x64); + + // generate a subset of the unmasked values: + psVector *tmp = psVectorAllocEmpty (x64->n, PS_TYPE_F64); + for (int itmp = 0; itmp < x64->n; itmp++) { + if (mask && (mask->data.PS_TYPE_VECTOR_MASK_DATA[itmp] && maskValue)) continue; + psVectorAppend (tmp, x64->data.F64[itmp]); + } + psVector *sorted = psVectorSort(NULL,tmp); median = sorted->data.F64[sorted->n / 2]; // CZW: I'm not bothering to scale this because it doesn't really matter. sigma = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]); psFree(sorted); + psFree(tmp); if ((!isfinite(median))|| @@ -721,9 +729,16 @@ sigma = 1.0; } - if ((median == 0.0)&&(sigma == 1.0)) { - scale = false; + if (fabs(median) < 1e-10) { + // CZW 2014-10-21: This median is small and close to zero. This can cause issues with the + // scaling, + median_zero = true; + if ((sigma == 1.0)||(fabs(sigma) <= 1e-10)) { + // Don't bother scaling if sigma is unity (it's already scaled) or if the sigma calculation has gone wrong. + scale = false; + } } + // I can't see a way to not clobber x if it's already F64, so make a copy.x - z64 = psVectorCopy(NULL,x64,PS_TYPE_F64); + psVector *z64 = psVectorCopy(NULL,x64,PS_TYPE_F64); psBinaryOp(z64,z64,"-",psScalarAlloc(median,PS_TYPE_F64)); psBinaryOp(z64,z64,"/",psScalarAlloc(sigma,PS_TYPE_F64)); @@ -733,5 +748,4 @@ #endif - result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64); psFree(z64); // Done with this. @@ -758,22 +772,28 @@ #endif } + for (psS32 i = 0; i <= poly->nX; i++) { poly->coeff[i] = 0.0; poly->coeffErr[i] = 0.0; - - for (psS32 j = 0; j <= poly->nX; j++) { + if (median_zero) { // If the median is zero, the obtained solution just needs to be scaled by the sigma values. + poly->coeff[i] = Zcoeff[i] * pow(1.0 / sigma,i); + poly->coeffErr[i] = ZcoeffErr[i] * pow(1.0 / sigma,i); + } + else { // Otherwise, do the correct transformations by expanding the (x-m)/s terms. + for (psS32 j = 0; j <= poly->nX; j++) { #if (CZW) - printf(" %d %d %f %f %f %f => %f\n", - i,j,Zcoeff[j], - pow(1.0 / sigma,j) * pow(-1,j - i), - pow(median,j - i), - 1.0 * psBinomialCoeff(j,i), - Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i) - ); + printf(" %d %d %f %f %f %f => %f\n", + i,j,Zcoeff[j], + pow(1.0 / sigma,j) * pow(-1,j - i), + pow(median,j - i), + 1.0 * psBinomialCoeff(j,i), + Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i) + ); #endif - poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i); - poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2); + poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i); + poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2); + } + poly->coeffErr[i] = sqrt(poly->coeffErr[i]); } - poly->coeffErr[i] = sqrt(poly->coeffErr[i]); #if (CZW) printf("poly1d: unscaled parameters: %d %f %f\n", @@ -784,5 +804,5 @@ psFree(ZcoeffErr); - } + } // End scaling block. break; Index: /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/mathtypes/psImage.c =================================================================== --- /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/mathtypes/psImage.c (revision 37826) +++ /branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/mathtypes/psImage.c (revision 37827) @@ -91,5 +91,5 @@ } - long numBytes = numRows * numCols * elementSize; + size_t numBytes = (size_t) numRows * (size_t) numCols * (size_t) elementSize; psImage* image = (psImage* ) p_psAlloc(file, lineno, func, sizeof(psImage));