IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1775


Ignore:
Timestamp:
Sep 9, 2004, 4:55:59 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib
Files:
3 added
11 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/dataManip/psFunctions.c

    r1765 r1775  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-09 23:39:36 $
     9 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-10 02:52:02 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    15931593    float width;
    15941594
     1595    if (fabs(max-min) < FLT_EPSILON) {
     1596        psLogMsg(__func__, PS_LOG_WARN, "min and max are the same.\n");
     1597        return(NULL);
     1598    }
     1599
    15951600    tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
    1596 
    15971601    tmp->n = numSplines;
    15981602
    15991603    tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    16001604    for (i=0;i<numSplines;i++) {
    1601         (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     1605        (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    16021606    }
    16031607
     
    16261630
    16271631    numSplines = bounds->n - 1;
     1632    tmp->n = numSplines;
    16281633
    16291634    tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    16301635    for (i=0;i<numSplines;i++) {
    1631         (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     1636        (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    16321637    }
    16331638
     
    16681673    mid = (max-min)/2;
    16691674
    1670     while (min != max) {
     1675    while (min+1 < max) {
     1676        psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max);
    16711677        mid = (max-min)/2;
    16721678
     
    17211727
    17221728
    1723 psSpline1D *psSpline1DEGen(const psSpline1D *spline,
    1724                            psVector *data)
     1729psSpline1D *psSpline1DGen(const psSpline1D *spline,
     1730                          psVector *data)
    17251731{
    17261732    int i;
     
    17341740    }
    17351741    // XXX: verify your math
    1736     if (1 == (spline->spline[0])->n) {
     1742    if (2 == (spline->spline[0])->n) {
    17371743        for (i=0;i<spline->n;i++) {
    17381744            slope = (data->data.F32[i+1] - data->data.F32[i]) /
     
    17421748            (spline->spline[0])->coeff[1] = slope;
    17431749        }
    1744 
    1745         psLogMsg(__func__, PS_LOG_WARN,
    1746                  "Linear splines have not been implemented.");
    1747         return(NULL);
    1748     } else if (3 == (spline->spline[0])->n) {
     1750        return((psSpline1D *) spline);
     1751    } else if (4 == (spline->spline[0])->n) {
    17491752        psLogMsg(__func__, PS_LOG_WARN,
    17501753                 "Cubic splines have not been implemented.");
     
    17531756        psLogMsg(__func__, PS_LOG_WARN,
    17541757                 "Don't know how to generate %d-order splines.",
    1755                  (spline->spline[0])->n);
     1758                 (spline->spline[0])->n-1);
    17561759        return(NULL);
    17571760    }
  • trunk/psLib/src/dataManip/psFunctions.h

    r1721 r1775  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-08 06:04:10 $
     14*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-10 02:52:02 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    402402                               const psSpline1D *spline);
    403403
     404psSpline1D *psSpline1DGen(const psSpline1D *spline,
     405                          psVector *data);
     406
     407
    404408/* \} */// End of MathGroup Functions
    405409
  • trunk/psLib/src/dataManip/psMatrix.c

    r1440 r1775  
    2121 *  @author Ross Harman, MHPCC
    2222 *   
    23  *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
    24  *  @date $Date: 2004-08-09 23:34:57 $
     23 *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
     24 *  @date $Date: 2004-09-10 02:52:02 $
    2525 *
    2626 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    205205}
    206206
    207 psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage, const psVector* inVector, const
    208                           psVector* inPerm)
     207psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage,
     208                          const psVector* inVector, const psVector* inPerm)
    209209{
    210210    int arraySize = 0;
  • trunk/psLib/src/dataManip/psMinimize.c

    r1734 r1775  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-08 23:32:03 $
     11 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-10 02:52:02 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    846846/******************************************************************************
    847847GUS
    848  
    849 XXX: bug report for the psMatrix covar argument.
    850848 *****************************************************************************/
    851849bool psMinimizeLM(psMinimization *min,
     
    856854                  psMinimizeLMFunc func)
    857855{
    858     /*
    859         psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
    860         psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    861         psVector *deriv;
    862         int i;
    863         int k;
    864         float tmp;
    865      
    866         tmp = func(deriv, params, coords);
    867         for (k=0;k<params->n;k++) {
    868             beta->data.F32[k] = 0.0;
     856    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
     857    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
     858    psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
     859    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     860    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     861    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     862    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
     863    psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
     864    int i;
     865    int j;
     866    int k;
     867    float newValue;
     868    float oldValue;
     869    float lamda = 1.0;
     870
     871    min->lastDelta = -HUGE;
     872    min->iter = 0;
     873    min->value = func(deriv, params, coords);
     874
     875    while ((min->lastDelta > min->tol) &&
     876            (min->iter < min->maxIter)) {
     877
     878        for (i=0;i<params->n;i++) {
     879            beta->data.F32[i] = deriv->data.F32[i];
     880        }
     881
     882        for (j=0;j<params->n;j++) {
     883            for (k=0;k<params->n;k++) {
     884                alpha->data.F32[j][k] = deriv->data.F32[j] * deriv->data.F32[j];
     885            }
     886        }
     887        for (j=0;j<params->n;j++) {
     888            for (k=0;k<params->n;k++) {
     889                if (j != k) {
     890                    A->data.F32[j][k] = alpha->data.F32[j][k];
     891                } else {
     892                    A->data.F32[j][k] = (1.0 + lamda) * alpha->data.F32[j][k];
     893                }
     894            }
     895        }
     896        aOut = psMatrixLUD(aOut, perm, A);
     897        newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
     898
     899        oldValue = min->value;
     900        newValue = func(newDeriv, newParams, coords);
     901        if (oldValue > newValue) {
     902            min->lastDelta = oldValue - newValue;
     903            min->value = newValue;
     904
    869905            for (i=0;i<params->n;i++) {
    870                 beta->data.F32[k]+= deriv
    871             }
    872     */
     906                params->data.F32[i] = newParams->data.F32[i];
     907                deriv->data.F32[i] = newDeriv->data.F32[i];
     908            }
     909            min->value = func(deriv, params, coords);
     910
     911            lamda*= 0.1;
     912        } else {
     913            lamda*= 10.0;
     914        }
     915        min->iter++;
     916    }
     917    psFree(beta);
     918    psFree(perm);
     919    psFree(newParams);
     920    psFree(alpha);
     921    psFree(A);
     922    psFree(aOut);
     923    psFree(deriv);
     924    psFree(newDeriv);
     925
     926    if ((min->iter < min->maxIter) ||
     927            (min->lastDelta <= min->tol)) {
     928        return(true);
     929    }
     930
    873931    return(false);
    874932}
  • trunk/psLib/src/math/psMatrix.c

    r1440 r1775  
    2121 *  @author Ross Harman, MHPCC
    2222 *   
    23  *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
    24  *  @date $Date: 2004-08-09 23:34:57 $
     23 *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
     24 *  @date $Date: 2004-09-10 02:52:02 $
    2525 *
    2626 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    205205}
    206206
    207 psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage, const psVector* inVector, const
    208                           psVector* inPerm)
     207psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage,
     208                          const psVector* inVector, const psVector* inPerm)
    209209{
    210210    int arraySize = 0;
  • trunk/psLib/src/math/psMinimize.c

    r1734 r1775  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-08 23:32:03 $
     11 *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-10 02:52:02 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    846846/******************************************************************************
    847847GUS
    848  
    849 XXX: bug report for the psMatrix covar argument.
    850848 *****************************************************************************/
    851849bool psMinimizeLM(psMinimization *min,
     
    856854                  psMinimizeLMFunc func)
    857855{
    858     /*
    859         psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
    860         psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    861         psVector *deriv;
    862         int i;
    863         int k;
    864         float tmp;
    865      
    866         tmp = func(deriv, params, coords);
    867         for (k=0;k<params->n;k++) {
    868             beta->data.F32[k] = 0.0;
     856    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
     857    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
     858    psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
     859    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     860    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     861    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     862    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
     863    psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
     864    int i;
     865    int j;
     866    int k;
     867    float newValue;
     868    float oldValue;
     869    float lamda = 1.0;
     870
     871    min->lastDelta = -HUGE;
     872    min->iter = 0;
     873    min->value = func(deriv, params, coords);
     874
     875    while ((min->lastDelta > min->tol) &&
     876            (min->iter < min->maxIter)) {
     877
     878        for (i=0;i<params->n;i++) {
     879            beta->data.F32[i] = deriv->data.F32[i];
     880        }
     881
     882        for (j=0;j<params->n;j++) {
     883            for (k=0;k<params->n;k++) {
     884                alpha->data.F32[j][k] = deriv->data.F32[j] * deriv->data.F32[j];
     885            }
     886        }
     887        for (j=0;j<params->n;j++) {
     888            for (k=0;k<params->n;k++) {
     889                if (j != k) {
     890                    A->data.F32[j][k] = alpha->data.F32[j][k];
     891                } else {
     892                    A->data.F32[j][k] = (1.0 + lamda) * alpha->data.F32[j][k];
     893                }
     894            }
     895        }
     896        aOut = psMatrixLUD(aOut, perm, A);
     897        newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
     898
     899        oldValue = min->value;
     900        newValue = func(newDeriv, newParams, coords);
     901        if (oldValue > newValue) {
     902            min->lastDelta = oldValue - newValue;
     903            min->value = newValue;
     904
    869905            for (i=0;i<params->n;i++) {
    870                 beta->data.F32[k]+= deriv
    871             }
    872     */
     906                params->data.F32[i] = newParams->data.F32[i];
     907                deriv->data.F32[i] = newDeriv->data.F32[i];
     908            }
     909            min->value = func(deriv, params, coords);
     910
     911            lamda*= 0.1;
     912        } else {
     913            lamda*= 10.0;
     914        }
     915        min->iter++;
     916    }
     917    psFree(beta);
     918    psFree(perm);
     919    psFree(newParams);
     920    psFree(alpha);
     921    psFree(A);
     922    psFree(aOut);
     923    psFree(deriv);
     924    psFree(newDeriv);
     925
     926    if ((min->iter < min->maxIter) ||
     927            (min->lastDelta <= min->tol)) {
     928        return(true);
     929    }
     930
    873931    return(false);
    874932}
  • trunk/psLib/src/math/psPolynomial.c

    r1765 r1775  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-09 23:39:36 $
     9 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-10 02:52:02 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    15931593    float width;
    15941594
     1595    if (fabs(max-min) < FLT_EPSILON) {
     1596        psLogMsg(__func__, PS_LOG_WARN, "min and max are the same.\n");
     1597        return(NULL);
     1598    }
     1599
    15951600    tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
    1596 
    15971601    tmp->n = numSplines;
    15981602
    15991603    tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    16001604    for (i=0;i<numSplines;i++) {
    1601         (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     1605        (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    16021606    }
    16031607
     
    16261630
    16271631    numSplines = bounds->n - 1;
     1632    tmp->n = numSplines;
    16281633
    16291634    tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    16301635    for (i=0;i<numSplines;i++) {
    1631         (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     1636        (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    16321637    }
    16331638
     
    16681673    mid = (max-min)/2;
    16691674
    1670     while (min != max) {
     1675    while (min+1 < max) {
     1676        psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max);
    16711677        mid = (max-min)/2;
    16721678
     
    17211727
    17221728
    1723 psSpline1D *psSpline1DEGen(const psSpline1D *spline,
    1724                            psVector *data)
     1729psSpline1D *psSpline1DGen(const psSpline1D *spline,
     1730                          psVector *data)
    17251731{
    17261732    int i;
     
    17341740    }
    17351741    // XXX: verify your math
    1736     if (1 == (spline->spline[0])->n) {
     1742    if (2 == (spline->spline[0])->n) {
    17371743        for (i=0;i<spline->n;i++) {
    17381744            slope = (data->data.F32[i+1] - data->data.F32[i]) /
     
    17421748            (spline->spline[0])->coeff[1] = slope;
    17431749        }
    1744 
    1745         psLogMsg(__func__, PS_LOG_WARN,
    1746                  "Linear splines have not been implemented.");
    1747         return(NULL);
    1748     } else if (3 == (spline->spline[0])->n) {
     1750        return((psSpline1D *) spline);
     1751    } else if (4 == (spline->spline[0])->n) {
    17491752        psLogMsg(__func__, PS_LOG_WARN,
    17501753                 "Cubic splines have not been implemented.");
     
    17531756        psLogMsg(__func__, PS_LOG_WARN,
    17541757                 "Don't know how to generate %d-order splines.",
    1755                  (spline->spline[0])->n);
     1758                 (spline->spline[0])->n-1);
    17561759        return(NULL);
    17571760    }
  • trunk/psLib/src/math/psPolynomial.h

    r1722 r1775  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-08 06:04:10 $
     14*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-10 02:52:02 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    402402                               const psSpline1D *spline);
    403403
     404psSpline1D *psSpline1DGen(const psSpline1D *spline,
     405                          psVector *data);
     406
     407
    404408/* \} */// End of MathGroup Functions
    405409
  • trunk/psLib/src/math/psSpline.c

    r1765 r1775  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-09 23:39:36 $
     9 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-10 02:52:02 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    15931593    float width;
    15941594
     1595    if (fabs(max-min) < FLT_EPSILON) {
     1596        psLogMsg(__func__, PS_LOG_WARN, "min and max are the same.\n");
     1597        return(NULL);
     1598    }
     1599
    15951600    tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
    1596 
    15971601    tmp->n = numSplines;
    15981602
    15991603    tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    16001604    for (i=0;i<numSplines;i++) {
    1601         (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     1605        (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    16021606    }
    16031607
     
    16261630
    16271631    numSplines = bounds->n - 1;
     1632    tmp->n = numSplines;
    16281633
    16291634    tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    16301635    for (i=0;i<numSplines;i++) {
    1631         (tmp->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
     1636        (tmp->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
    16321637    }
    16331638
     
    16681673    mid = (max-min)/2;
    16691674
    1670     while (min != max) {
     1675    while (min+1 < max) {
     1676        psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max);
    16711677        mid = (max-min)/2;
    16721678
     
    17211727
    17221728
    1723 psSpline1D *psSpline1DEGen(const psSpline1D *spline,
    1724                            psVector *data)
     1729psSpline1D *psSpline1DGen(const psSpline1D *spline,
     1730                          psVector *data)
    17251731{
    17261732    int i;
     
    17341740    }
    17351741    // XXX: verify your math
    1736     if (1 == (spline->spline[0])->n) {
     1742    if (2 == (spline->spline[0])->n) {
    17371743        for (i=0;i<spline->n;i++) {
    17381744            slope = (data->data.F32[i+1] - data->data.F32[i]) /
     
    17421748            (spline->spline[0])->coeff[1] = slope;
    17431749        }
    1744 
    1745         psLogMsg(__func__, PS_LOG_WARN,
    1746                  "Linear splines have not been implemented.");
    1747         return(NULL);
    1748     } else if (3 == (spline->spline[0])->n) {
     1750        return((psSpline1D *) spline);
     1751    } else if (4 == (spline->spline[0])->n) {
    17491752        psLogMsg(__func__, PS_LOG_WARN,
    17501753                 "Cubic splines have not been implemented.");
     
    17531756        psLogMsg(__func__, PS_LOG_WARN,
    17541757                 "Don't know how to generate %d-order splines.",
    1755                  (spline->spline[0])->n);
     1758                 (spline->spline[0])->n-1);
    17561759        return(NULL);
    17571760    }
  • trunk/psLib/src/math/psSpline.h

    r1722 r1775  
    1212*  @author George Gusciora, MHPCC
    1313*
    14 *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
    15 *  @date $Date: 2004-09-08 06:04:10 $
     14*  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
     15*  @date $Date: 2004-09-10 02:52:02 $
    1616*
    1717*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    402402                               const psSpline1D *spline);
    403403
     404psSpline1D *psSpline1DGen(const psSpline1D *spline,
     405                          psVector *data);
     406
     407
    404408/* \} */// End of MathGroup Functions
    405409
  • trunk/psLib/test/dataManip/Makefile

    r1308 r1775  
    33##  Makefile:   test/sysUtils
    44##
    5 ##  $Revision: 1.35 $  $Name: not supported by cvs2svn $
    6 ##  $Date: 2004-07-28 02:50:38 $
     5##  $Revision: 1.36 $  $Name: not supported by cvs2svn $
     6##  $Date: 2004-09-10 02:55:20 $
    77##
    88##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     
    2222tst_psFunc00 \
    2323tst_psFunc01 \
     24tst_psFunc02 \
     25tst_psFunc03 \
     26tst_psFunc04 \
    2427tst_psHist00 \
    2528tst_psHist01 \
Note: See TracChangeset for help on using the changeset viewer.