Index: trunk/psModules/src/pmObjects.c
===================================================================
--- trunk/psModules/src/pmObjects.c	(revision 3078)
+++ trunk/psModules/src/pmObjects.c	(revision 3086)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-01-24 22:57:52 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-01-25 01:22:12 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22,19 +22,157 @@
 above the given threshold.  Returns a vector of type PS_TYPE_U32 containing
 the location (x value) of all peaks.
+ 
+XXX: What types should be supported?  Only F32 us implemented.
+ 
+XXX: We currently step through the input vector twice; once to determine the
+size of the output vector, then to set the values of the output vector.
+Depending upon actual use, this may need to be optimized.
  *****************************************************************************/
 psVector *pmFindVectorPeeks(const psVector *vector,
                             psF32 threshold)
 {
-    return(NULL);
+    PS_VECTOR_CHECK_NULL(vector, NULL);
+    PS_VECTOR_CHECK_TYPE(vector, PS_TYPE_F32, NULL);
+    int count = 0;
+    int n = vector->n;
+
+    if ((vector->data.F32[0] > vector->data.F32[1]) &&
+            (vector->data.F32[0] > threshold)) {
+        count++;
+    }
+
+    for (psU32 i = 1; i < n-1 ; i++) {
+        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
+                (vector->data.F32[i] > vector->data.F32[i+1]) &&
+                (vector->data.F32[i] > threshold)) {
+            count++;
+        }
+    }
+    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
+            (vector->data.F32[n-1] > threshold)) {
+        count++;
+    }
+
+    psVector *tmpVector = psVectorAlloc(count, PS_TYPE_U32);
+    count = 0;
+    if ((vector->data.F32[0] > vector->data.F32[1]) &&
+            (vector->data.F32[0] > threshold)) {
+        tmpVector->data.U32[count++] = 0;
+    }
+    for (psU32 i = 1; i < (n-1) ; i++) {
+        if ((vector->data.F32[i] > vector->data.F32[i-1]) &&
+                (vector->data.F32[i] > vector->data.F32[i+1]) &&
+                (vector->data.F32[i] > threshold)) {
+            tmpVector->data.U32[count++] = i;
+        }
+    }
+    if ((vector->data.F32[n-1] > vector->data.F32[n-2]) &&
+            (vector->data.F32[n-1] > threshold)) {
+        tmpVector->data.U32[count++] = n-1;
+    }
+
+    return(tmpVector);
+}
+
+/******************************************************************************
+XXX: Is there a better way to do this?
+ *****************************************************************************/
+psVector *p_psGetRowVectorFromImage(psImage *image,
+                                    psU32 row)
+{
+    PS_IMAGE_CHECK_NULL(image, NULL);
+    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
+
+    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
+    for (psU32 col = 0; col < image->numCols ; col++) {
+        tmpVector->data.F32[col] = image->data.F32[row][col];
+    }
+    return(tmpVector);
 }
 
 /******************************************************************************
 pmFindImagePeeks(image, threshold): Find all local peaks in the given psImage
-above the given threshold.  Returns a psList containing the location (x/y
-value) of all peaks.
+above the given threshold.  Returns a psList containing location (x/y value)
+of all peaks.
  *****************************************************************************/
 psList *pmFindImagePeeks(const psImage *image,
                          psF32 threshold)
 {
+    PS_IMAGE_CHECK_NULL(image, NULL);
+    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
+    psVector *tmpRow = NULL;
+    psU32 col = 0;
+    psU32 row = 0;
+    psVector *rowX = NULL;
+
+    tmpRow = p_psGetRowVectorFromImage((psImage *) image, 0);
+    psVector *row1 = pmFindVectorPeeks(tmpRow, threshold);
+    tmpRow = p_psGetRowVectorFromImage((psImage *) image, 1);
+    psVector *row2 = pmFindVectorPeeks(tmpRow, threshold);
+    tmpRow = p_psGetRowVectorFromImage((psImage *) image, 2);
+    psVector *row3 = pmFindVectorPeeks(tmpRow, threshold);
+
+    row = 0;
+    for (psU32 i = 0 ; i < row1->n ; i++ ) {
+        col = row1->data.U32[i];
+        if (col == 0) {
+            if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
+                // Add peak at location (row, col)
+            }
+        } else if (col < (image->numCols - 1)) {
+            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col+2])) {
+                // Add peak at location (row, col)
+            }
+
+        } else if (col == (image->numCols - 1)) {
+            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
+                // Add peak at location (row, col)
+            }
+
+        } else {
+            printf("XXX: ERROR!\n");
+        }
+    }
+
+    row = image->numRows - 1;
+    for (psU32 i = 0 ; i < row1->n ; i++ ) {
+        col = row1->data.U32[i];
+        if (col == 0) {
+            if ( (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col+1])) {
+                // Add peak at location (row, col)
+            }
+        } else if (col < (image->numCols - 1)) {
+            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row][col+1]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col+2])) {
+                // Add peak at location (row, col)
+            }
+
+        } else if (col == (image->numCols - 1)) {
+            if ( (image->data.F32[row][col] >= image->data.F32[row][col-1]) &&
+                    (image->data.F32[row][col] > image->data.F32[row+1][col]) &&
+                    (image->data.F32[row][col] >= image->data.F32[row+1][col-1])) {
+                // Add peak at location (row, col)
+            }
+
+        } else {
+            printf("XXX: ERROR!\n");
+        }
+    }
+
+
+
     return(NULL);
 }
@@ -141,4 +279,14 @@
 /******************************************************************************
 XXX: Why only *x argument?
+ 
+May x[0] is x and x[1] is y?
+ 
+    params->data.F32[0] = So;
+    params->data.F32[1] = Zo;
+    params->data.F32[2] = Xo;
+    params->data.F32[3] = Yo;
+    params->data.F32[4] = sqrt(2.0) / SigmaX;
+    params->data.F32[5] = sqrt(2.0) / SigmaY;
+    params->data.F32[6] = Sxy;
  *****************************************************************************/
 psF32 pmMinLM_Gauss2D(psVector *deriv,
@@ -146,9 +294,32 @@
                       psVector *x)
 {
-    return(0.0);
-}
-
-/******************************************************************************
-XXX: Why only *x argument?
+    psF32 X = x->data.F32[0] - params->data.F32[2];
+    psF32 Y = x->data.F32[1] - params->data.F32[3];
+    psF32 px = params->data.F32[4]*X;
+    psF32 py = params->data.F32[5]*Y;
+    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
+    psF32 r = exp(-z);
+    psF32 f = params->data.F32[1]*r + params->data.F32[0];
+
+    psF32 q = params->data.F32[1]*r;
+    deriv->data.F32[0] = +1.0;
+    deriv->data.F32[1] = +r;
+    deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y);
+    deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X);
+    deriv->data.F32[4] = -2.0*q*px*X;
+    deriv->data.F32[5] = -2.0*q*py*Y;
+    deriv->data.F32[6] = -q*X*Y;
+
+    return(f);
+}
+
+/******************************************************************************
+    params->data.F32[0] = So;
+    params->data.F32[1] = Zo;
+    params->data.F32[2] = Xo;
+    params->data.F32[3] = Yo;
+    params->data.F32[4] = sqrt(2) / SigmaX;
+    params->data.F32[5] = sqrt(2) / SigmaY;
+    params->data.F32[6] = Sxy;
  *****************************************************************************/
 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv,
@@ -156,8 +327,40 @@
                             psVector *x)
 {
-    return(0.0);
-}
-
-/******************************************************************************
+    psF32 X = x->data.F32[0] - params->data.F32[2];
+    psF32 Y = x->data.F32[1] - params->data.F32[3];
+    psF32 px = params->data.F32[4]*X;
+    psF32 py = params->data.F32[5]*Y;
+    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
+    psF32 t = 1 + z + 0.5*z*z;
+    psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */
+    psF32 f = params->data.F32[1]*r + params->data.F32[0];
+
+    //
+    // note difference from a pure gaussian: q = params->data.F32[1]*r
+    //
+
+    psF32 q = params->data.F32[1]*r*r*t;
+    deriv->data.F32[0] = +1.0;
+    deriv->data.F32[1] = +r;
+    deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
+    deriv->data.F32[3] = q *
+                         (2.0*py*params->data.F32[5] + params->data.F32[6]*X);
+    deriv->data.F32[4] = -2.0*q*px*X;
+    deriv->data.F32[5] = -2.0*q*py*Y;
+    deriv->data.F32[6] = -q*X*Y;
+
+    return(f);
+}
+
+/******************************************************************************
+    params->data.F32[0] = So;
+    params->data.F32[1] = Zo;
+    params->data.F32[2] = Xo;
+    params->data.F32[3] = Yo;
+    params->data.F32[4] = Sx;
+    params->data.F32[5] = Sy;
+    params->data.F32[6] = Sxy;
+    params->data.F32[7] = B2;
+    params->data.F32[8] = B3;
  *****************************************************************************/
 psF32 pmMinLM_Wauss2D(psVector *deriv,
@@ -165,8 +368,48 @@
                       psVector *x)
 {
-    return(0.0);
-}
-
-/******************************************************************************
+    psF32 X = x->data.F32[0] - params->data.F32[2];
+    psF32 Y = x->data.F32[1] - params->data.F32[2];
+    psF32 px = params->data.F32[4]*X;
+    psF32 py = params->data.F32[5]*Y;
+    psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y;
+    psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0);
+    psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */
+    psF32 f = params->data.F32[1]*r + params->data.F32[0];
+
+    //
+    // note difference from gaussian: q = params->data.F32[1]*r
+    //
+    psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0));
+    deriv->data.F32[0] = +1.0;
+    deriv->data.F32[1] = +r;
+    deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y);
+    deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X);
+    deriv->data.F32[4] = -2.0*q*px*X;
+    deriv->data.F32[5] = -2.0*q*py*Y;
+    deriv->data.F32[6] = -q*X*Y;
+    deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t;
+    deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0;
+    //
+    // The values of 100 dampen the swing of params->data.F32[7,8] */
+    //
+
+    return(f);
+}
+
+// XXX: What should these be?
+#define FFACTOR 1.0
+#define FSCALE 1.0
+/******************************************************************************
+    params->data.F32[0] = So;
+    params->data.F32[1] = Zo;
+    params->data.F32[2] = Xo;
+    params->data.F32[3] = Yo;
+    params->data.F32[4] = SxInner;
+    params->data.F32[5] = SyInner;
+    params->data.F32[6] = SxyInner;
+    params->data.F32[7] = SxOuter;
+    params->data.F32[8] = SyOuter;
+    params->data.F32[9] = SxyOuter;
+    params->data.F32[10] = N;
  *****************************************************************************/
 psF32 pmMinLM_TwistGauss2D(psVector *deriv,
@@ -174,8 +417,57 @@
                            psVector *x)
 {
-    return(0.0);
-}
-
-/******************************************************************************
+    psF32 X = x->data.F32[0] - params->data.F32[2];
+    psF32 Y = x->data.F32[1] - params->data.F32[3];
+    psF32 px1 = params->data.F32[4]*X;
+    psF32 py1 = params->data.F32[5]*Y;
+    psF32 px2 = params->data.F32[7]*X;
+    psF32 py2 = params->data.F32[8]*Y;
+    psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y;
+    psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y;
+    psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10]));
+
+
+    psF32 f = params->data.F32[5]*r + params->data.F32[6];
+    psF32 q1 = params->data.F32[5]*PS_SQR(r);
+    psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0));
+    deriv->data.F32[0] = +1.0;
+    deriv->data.F32[1] = +r;
+    deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y);
+    deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X);
+
+    //
+    // These fudge factors impede the growth of params->data.F32[4] beyond
+    // params->data.F32[7].
+    //
+    psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]);
+    psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
+    deriv->data.F32[4] = -2.0*q1*px1*X*f2;
+
+    //
+    // These fudge factors impede the growth of params->data.F32[5] beyond
+    // params->data.F32[8].
+    //
+    f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]);
+    f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0;
+    deriv->data.F32[5] = -2.0*q1*py1*Y*f2;
+    deriv->data.F32[6] = -q1*X*Y;
+    deriv->data.F32[7] = -2.0*q2*px2*X;
+    deriv->data.F32[8] = -2.0*q2*py2*Y;
+    deriv->data.F32[9] = -q2*X*Y;
+    deriv->data.F32[10] = -q1*log(z2);
+
+    return(f);
+}
+
+/******************************************************************************
+    float Sersic()
+    params->data.F32[0] = So;
+    params->data.F32[1] = Zo;
+    params->data.F32[2] = Xo;
+    params->data.F32[3] = Yo;
+    params->data.F32[4] = Sx;
+    params->data.F32[5] = Sy;
+    params->data.F32[6] = Sxy;
+    params->data.F32[7] = Nexp;
  *****************************************************************************/
 psF32 pmMinLM_Sersic(psVector *deriv,
@@ -187,4 +479,17 @@
 
 /******************************************************************************
+    float SersicBulge()
+    params->data.F32[0] So;
+    params->data.F32[1] Zo;
+    params->data.F32[2] Xo;
+    params->data.F32[3] Yo;
+    params->data.F32[4] SxInner;
+    params->data.F32[5] SyInner;
+    params->data.F32[6] SxyInner;
+    params->data.F32[7] Zd;
+    params->data.F32[8] SxOuter;
+    params->data.F32[9] SyOuter;
+    params->data.F32[10] = SxyOuter;
+    params->data.F32[11] = Nexp;
  *****************************************************************************/
 psF32 pmMinLM_SersicCore(psVector *deriv,