Changeset 37827

Timestamp:

Jan 12, 2015, 12:35:53 PM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ipp-pv3-20140717-merge/psLib

Files:

• share/tai_utc.dat (modified) (1 diff)
• share/tai_utc.raw (modified) (1 diff)
• src/imageops/psImageMapFit.c (modified) (9 diffs)
• src/imageops/psImageMapFit.h (modified) (1 diff)
• src/math/psEllipse.c (modified) (1 diff)
• src/math/psEllipse.h (modified) (1 diff)
• src/math/psMinimizePolyFit.c (modified) (5 diffs, 1 prop)
• src/mathtypes/psImage.c (modified) (1 diff)

branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.dat

@@ -51 +51 @@
         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

branches/eam_branches/ipp-pv3-20140717-merge/psLib/share/tai_utc.raw

@@ -38 +38 @@
  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

branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.c

@@ -18 +18 @@
 
 #include <stdio.h>
+#include <stdlib.h>
 #include "psError.h"
 #include "psAbort.h"
@@ -35 +36 @@
 #include "psImageStructManip.h"
 #include "psImageMap.h"
+#include "psSparse.h"
 // #include "psImagePixelInterpolate.h"
 // #include "psImageUnbin.h"
@@ -118 +120 @@
     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 +264 @@
             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 +275 @@
                     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 +285 @@
     }
     // 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 +330 @@
         return true;
     }
-
+    
     // set bad values to NaN
     for (int i = 0; i < Nx*Ny; i++) {
@@ -341 +351 @@
     psFree (B);
     psFree (Empty);
-
     *pGoodFit = true;
     return true;
@@ -402 +411 @@
         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);

branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/imageops/psImageMapFit.h

@@ -20 +20 @@
 // 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,

branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.c

@@ -134 +134 @@
     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;

branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psEllipse.h

@@ -65 +65 @@
                                    );
 
+/// Convert shape to axes representation and compute errors
+psEllipseAxes psEllipseShapeToAxesWithErrors(psEllipseShape shape, ///< Shape of ellipse
+                                   psEllipseShape errors, ///<errors on shape params
+                                   double maxAR, ///< Maximum allowed axis ratio
+                                   psEllipseAxes *pShapeErrors ///< propagated errors on axes
+                                   );
+
 /// Convert axes to polarization representation
 psEllipsePol psEllipseAxesToPol(psEllipseAxes axes ///< Axes of ellipse

branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/math/psMinimizePolyFit.c

@@ -700 +700 @@
 
     // Define values that may be used by PS_POLYNOMIAL_ORD form.
+    // these scaling values put the dynamic range of the data into some vaguely sensible location
     bool scale = false;
+    bool median_zero = false;
     double median = 0.0;
     double sigma = 1.0;
-    psVector *sorted = NULL;
-    psVector *z64 = NULL;
-    
+
     switch (poly->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 +729 @@
           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 +748 @@
 #endif
 
-        
         result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
         psFree(z64); // Done with this.
@@ -758 +772 @@
 #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 +804 @@
         psFree(ZcoeffErr);
 
-      }
+      } // End scaling block.
         
         break;

branches/eam_branches/ipp-pv3-20140717-merge/psLib/src/mathtypes/psImage.c

@@ -91 +91 @@
     }
 
-    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));