IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1781


Ignore:
Timestamp:
Sep 10, 2004, 12:26:09 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psFunctions.c

    r1775 r1781  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-10 02:52:02 $
     9 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-10 22:25:59 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    441441    float xSum = 1.0;
    442442
     443    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     444            "---- Calling p_psOrdPolynomial1DEval(%f)\n", x);
     445    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     446            "Polynomial order is %d\n", myPoly->n);
     447    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
     448        psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     449                "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]);
     450    }
     451
     452
    443453    if (NULL == myPoly) {
    444454        psAbort(__func__, "psPolynomial1DEval(): myPoly is NULL\n");
    … …  
    455465
    456466    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
    457         if (myPoly->mask[loop_x] != 0) {
     467        if (myPoly->mask[loop_x] == 0) {
    458468            polySum += xSum * myPoly->coeff[loop_x];
    459469            xSum *= x;
    … …  
    535545        ySum = xSum;
    536546        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    537             if (myPoly->mask[loop_x][loop_y] != 0) {
     547            if (myPoly->mask[loop_x][loop_y] == 0) {
    538548                polySum += ySum * myPoly->coeff[loop_x][loop_y];
    539549                ySum *= y;
    … …  
    564574    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
    565575        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    566             if (myPoly->mask[loop_x][loop_y] != 0) {
     576            if (myPoly->mask[loop_x][loop_y] == 0) {
    567577                polySum += myPoly->coeff[loop_x][loop_y] *
    568578                           psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    639649            zSum = ySum;
    640650            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    641                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     651                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    642652                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
    643653                    zSum *= z;
    … …  
    675685        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    676686            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    677                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     687                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    678688                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
    679689                               psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    768778                zSum = ySum;
    769779                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    770                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     780                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    771781                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
    772782                        zSum *= z;
    … …  
    811821            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    812822                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    813                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     823                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    814824                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
    815825                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
    … …  
    11291139
    11301140    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
    1131         if (myPoly->mask[loop_x] != 0) {
     1141        if (myPoly->mask[loop_x] == 0) {
    11321142            polySum += xSum * myPoly->coeff[loop_x];
    11331143            xSum *= x;
    … …  
    12101220        ySum = xSum;
    12111221        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    1212             if (myPoly->mask[loop_x][loop_y] != 0) {
     1222            if (myPoly->mask[loop_x][loop_y] == 0) {
    12131223                polySum += ySum * myPoly->coeff[loop_x][loop_y];
    12141224                ySum *= y;
    … …  
    12391249    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
    12401250        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    1241             if (myPoly->mask[loop_x][loop_y] != 0) {
     1251            if (myPoly->mask[loop_x][loop_y] == 0) {
    12421252                polySum += myPoly->coeff[loop_x][loop_y] *
    12431253                           psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    13131323            zSum = ySum;
    13141324            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1315                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1325                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    13161326                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
    13171327                    zSum *= z;
    … …  
    13491359        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    13501360            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1351                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1361                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    13521362                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
    13531363                               psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    14441454                zSum = ySum;
    14451455                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1446                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1456                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    14471457                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
    14481458                        zSum *= z;
    … …  
    14871497            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    14881498                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1489                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1499                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    14901500                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
    14911501                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
    … …  
    16081618    tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
    16091619    width = (max - min) / ((float) numSplines);
    1610     tmpDomain = min;
    1611     for (i=0;i<numSplines+1;i++) {
     1620
     1621    (tmp->domains)[0] = min;
     1622    tmpDomain = min+width;
     1623    for (i=1;i<numSplines+1;i++) {
    16121624        (tmp->domains)[i] = tmpDomain;
    16131625        tmpDomain+= width;
    16141626    }
     1627    (tmp->domains)[numSplines] = max;
    16151628
    16161629    return(tmp);
    … …  
    16611674    int mid;
    16621675
     1676    psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1677            "---- Calling VectorBinDisect(%f)\n", x);
     1678
    16631679    if ((x < bins[0]) ||
    16641680            (x > bins[numBins-1])) {
    … …  
    16701686
    16711687    min = 0;
    1672     max = numBins-1;
    1673     mid = (max-min)/2;
    1674 
    1675     while (min+1 < max) {
    1676         psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max);
    1677         mid = (max-min)/2;
    1678 
    1679         if (x < bins[mid]) {
    1680             max = mid;
     1688    max = numBins-2;
     1689    mid = ((max+1)-min)/2;
     1690
     1691    while (min != max) {
     1692        psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1693                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
     1694                min, mid, max, x, bins[mid]);
     1695
     1696        if (x == bins[mid]) {
     1697            psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1698                    "---- Exiting VectorBinDisect(): bin %d\n", min);
     1699            return(min);
     1700        } else if (x < bins[mid]) {
     1701            max = mid-1;
    16811702        } else {
    16821703            min = mid;
    16831704        }
    1684     }
    1685 
     1705        mid = ((max+1)+min)/2;
     1706    }
     1707
     1708    psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1709            "---- Exiting VectorBinDisect(): bin %d\n", min);
    16861710    return(min);
    16871711}
    … …  
    17331757    float slope;
    17341758
     1759    psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1760            "---- Calling psSpline1DGen()\n");
     1761
    17351762    if (data->n != (1 + spline->n)) {
    17361763        psLogMsg(__func__, PS_LOG_WARN,
    … …  
    17441771            slope = (data->data.F32[i+1] - data->data.F32[i]) /
    17451772                    (spline->domains[i+1] - spline->domains[i]);
    1746             (spline->spline[0])->coeff[0] = data->data.F32[i] -
     1773            (spline->spline[i])->coeff[0] = data->data.F32[i] -
    17471774                                            (slope * spline->domains[i]);
    1748             (spline->spline[0])->coeff[1] = slope;
    1749         }
     1775            //printf("HMMM (%f, %f, %f)\n", data->data.F32[i], slope, spline->domains[i]);
     1776
     1777            (spline->spline[i])->coeff[1] = slope;
     1778            psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1779                    "---- spline %d coeffs are (%f, %f)\n", i,
     1780                    (spline->spline[i])->coeff[0],
     1781                    (spline->spline[i])->coeff[1]);
     1782        }
     1783        psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1784                "---- Exiting psSpline1DGen()\n");
    17501785        return((psSpline1D *) spline);
    17511786    } else if (4 == (spline->spline[0])->n) {

  • trunk/psLib/src/dataManip/psMinimize.c

    r1775 r1781  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-10 02:52:02 $
     11 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-10 22:25:59 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    845845}
    846846/******************************************************************************
     847This routine will take an procedure which calculates an arbitrary function
     848and it's derivative and minimize it.
     849 
    847850GUS
     851 
     852XXX: Do this:
     853 After checking that all entries in the paramMask are 1 or 0, when
     854 forming the A matrix from alpha, try this:
     855 
     856     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     857 
     858 
    848859 *****************************************************************************/
    849860bool psMinimizeLM(psMinimization *min,
    … …  
    857868    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
    858869    psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
     870    psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
    859871    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    860872    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    … …  
    869881    float lamda = 1.0;
    870882
     883    for (i=0;i<params->n;i++) {
     884        origParams->data.F32[i] = params->data.F32[i];
     885    }
     886
    871887    min->lastDelta = -HUGE;
    872888    min->iter = 0;
    873889    min->value = func(deriv, params, coords);
     890    for (i=0;i<params->n;i++) {
     891        if ((paramMask != NULL) &&
     892                (paramMask->data.U8[i] != 0)) {
     893            deriv->data.F32[i] = 0.0;
     894        }
     895    }
    874896
    875897    while ((min->lastDelta > min->tol) &&
    … …  
    897919        newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
    898920
     921        for (i=0;i<params->n;i++) {
     922            if ((paramMask != NULL) &&
     923                    (paramMask->data.U8[i] != 0)) {
     924                newParams->data.F32[i] = origParams->data.F32[i];
     925            }
     926        }
     927
    899928        oldValue = min->value;
    900929        newValue = func(newDeriv, newParams, coords);
     930        for (i=0;i<params->n;i++) {
     931            if ((paramMask != NULL) &&
     932                    (paramMask->data.U8[i] != 0)) {
     933                deriv->data.F32[i] = 0.0;
     934            }
     935        }
     936
    901937        if (oldValue > newValue) {
    902938            min->lastDelta = oldValue - newValue;
    903939            min->value = newValue;
    904940
     941            // No need to check the paramMask here since we already did so
     942            // before the last function evaluation.
    905943            for (i=0;i<params->n;i++) {
    906944                params->data.F32[i] = newParams->data.F32[i];
    … …  
    918956    psFree(perm);
    919957    psFree(newParams);
     958    psFree(origParams);
    920959    psFree(alpha);
    921960    psFree(A);
    … …  
    17691808 
    17701809The basic idea is calculate chi-squared for a set of params/coords/errors.
    1771 This functions uses globale variables to receive the function pointer, the
     1810This functions uses global variables to receive the function pointer, the
    17721811data values, and the data errors.
    17731812 *****************************************************************************/

  • trunk/psLib/src/math/psMinimize.c

    r1775 r1781  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-10 02:52:02 $
     11 *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-10 22:25:59 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    845845}
    846846/******************************************************************************
     847This routine will take an procedure which calculates an arbitrary function
     848and it's derivative and minimize it.
     849 
    847850GUS
     851 
     852XXX: Do this:
     853 After checking that all entries in the paramMask are 1 or 0, when
     854 forming the A matrix from alpha, try this:
     855 
     856     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
     857 
     858 
    848859 *****************************************************************************/
    849860bool psMinimizeLM(psMinimization *min,
    … …  
    857868    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
    858869    psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
     870    psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
    859871    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    860872    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    … …  
    869881    float lamda = 1.0;
    870882
     883    for (i=0;i<params->n;i++) {
     884        origParams->data.F32[i] = params->data.F32[i];
     885    }
     886
    871887    min->lastDelta = -HUGE;
    872888    min->iter = 0;
    873889    min->value = func(deriv, params, coords);
     890    for (i=0;i<params->n;i++) {
     891        if ((paramMask != NULL) &&
     892                (paramMask->data.U8[i] != 0)) {
     893            deriv->data.F32[i] = 0.0;
     894        }
     895    }
    874896
    875897    while ((min->lastDelta > min->tol) &&
    … …  
    897919        newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
    898920
     921        for (i=0;i<params->n;i++) {
     922            if ((paramMask != NULL) &&
     923                    (paramMask->data.U8[i] != 0)) {
     924                newParams->data.F32[i] = origParams->data.F32[i];
     925            }
     926        }
     927
    899928        oldValue = min->value;
    900929        newValue = func(newDeriv, newParams, coords);
     930        for (i=0;i<params->n;i++) {
     931            if ((paramMask != NULL) &&
     932                    (paramMask->data.U8[i] != 0)) {
     933                deriv->data.F32[i] = 0.0;
     934            }
     935        }
     936
    901937        if (oldValue > newValue) {
    902938            min->lastDelta = oldValue - newValue;
    903939            min->value = newValue;
    904940
     941            // No need to check the paramMask here since we already did so
     942            // before the last function evaluation.
    905943            for (i=0;i<params->n;i++) {
    906944                params->data.F32[i] = newParams->data.F32[i];
    … …  
    918956    psFree(perm);
    919957    psFree(newParams);
     958    psFree(origParams);
    920959    psFree(alpha);
    921960    psFree(A);
    … …  
    17691808 
    17701809The basic idea is calculate chi-squared for a set of params/coords/errors.
    1771 This functions uses globale variables to receive the function pointer, the
     1810This functions uses global variables to receive the function pointer, the
    17721811data values, and the data errors.
    17731812 *****************************************************************************/

  • trunk/psLib/src/math/psPolynomial.c

    r1775 r1781  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-10 02:52:02 $
     9 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-10 22:25:59 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    441441    float xSum = 1.0;
    442442
     443    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     444            "---- Calling p_psOrdPolynomial1DEval(%f)\n", x);
     445    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     446            "Polynomial order is %d\n", myPoly->n);
     447    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
     448        psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     449                "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]);
     450    }
     451
     452
    443453    if (NULL == myPoly) {
    444454        psAbort(__func__, "psPolynomial1DEval(): myPoly is NULL\n");
    … …  
    455465
    456466    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
    457         if (myPoly->mask[loop_x] != 0) {
     467        if (myPoly->mask[loop_x] == 0) {
    458468            polySum += xSum * myPoly->coeff[loop_x];
    459469            xSum *= x;
    … …  
    535545        ySum = xSum;
    536546        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    537             if (myPoly->mask[loop_x][loop_y] != 0) {
     547            if (myPoly->mask[loop_x][loop_y] == 0) {
    538548                polySum += ySum * myPoly->coeff[loop_x][loop_y];
    539549                ySum *= y;
    … …  
    564574    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
    565575        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    566             if (myPoly->mask[loop_x][loop_y] != 0) {
     576            if (myPoly->mask[loop_x][loop_y] == 0) {
    567577                polySum += myPoly->coeff[loop_x][loop_y] *
    568578                           psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    639649            zSum = ySum;
    640650            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    641                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     651                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    642652                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
    643653                    zSum *= z;
    … …  
    675685        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    676686            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    677                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     687                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    678688                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
    679689                               psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    768778                zSum = ySum;
    769779                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    770                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     780                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    771781                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
    772782                        zSum *= z;
    … …  
    811821            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    812822                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    813                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     823                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    814824                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
    815825                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
    … …  
    11291139
    11301140    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
    1131         if (myPoly->mask[loop_x] != 0) {
     1141        if (myPoly->mask[loop_x] == 0) {
    11321142            polySum += xSum * myPoly->coeff[loop_x];
    11331143            xSum *= x;
    … …  
    12101220        ySum = xSum;
    12111221        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    1212             if (myPoly->mask[loop_x][loop_y] != 0) {
     1222            if (myPoly->mask[loop_x][loop_y] == 0) {
    12131223                polySum += ySum * myPoly->coeff[loop_x][loop_y];
    12141224                ySum *= y;
    … …  
    12391249    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
    12401250        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    1241             if (myPoly->mask[loop_x][loop_y] != 0) {
     1251            if (myPoly->mask[loop_x][loop_y] == 0) {
    12421252                polySum += myPoly->coeff[loop_x][loop_y] *
    12431253                           psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    13131323            zSum = ySum;
    13141324            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1315                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1325                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    13161326                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
    13171327                    zSum *= z;
    … …  
    13491359        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    13501360            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1351                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1361                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    13521362                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
    13531363                               psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    14441454                zSum = ySum;
    14451455                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1446                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1456                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    14471457                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
    14481458                        zSum *= z;
    … …  
    14871497            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    14881498                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1489                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1499                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    14901500                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
    14911501                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
    … …  
    16081618    tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
    16091619    width = (max - min) / ((float) numSplines);
    1610     tmpDomain = min;
    1611     for (i=0;i<numSplines+1;i++) {
     1620
     1621    (tmp->domains)[0] = min;
     1622    tmpDomain = min+width;
     1623    for (i=1;i<numSplines+1;i++) {
    16121624        (tmp->domains)[i] = tmpDomain;
    16131625        tmpDomain+= width;
    16141626    }
     1627    (tmp->domains)[numSplines] = max;
    16151628
    16161629    return(tmp);
    … …  
    16611674    int mid;
    16621675
     1676    psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1677            "---- Calling VectorBinDisect(%f)\n", x);
     1678
    16631679    if ((x < bins[0]) ||
    16641680            (x > bins[numBins-1])) {
    … …  
    16701686
    16711687    min = 0;
    1672     max = numBins-1;
    1673     mid = (max-min)/2;
    1674 
    1675     while (min+1 < max) {
    1676         psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max);
    1677         mid = (max-min)/2;
    1678 
    1679         if (x < bins[mid]) {
    1680             max = mid;
     1688    max = numBins-2;
     1689    mid = ((max+1)-min)/2;
     1690
     1691    while (min != max) {
     1692        psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1693                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
     1694                min, mid, max, x, bins[mid]);
     1695
     1696        if (x == bins[mid]) {
     1697            psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1698                    "---- Exiting VectorBinDisect(): bin %d\n", min);
     1699            return(min);
     1700        } else if (x < bins[mid]) {
     1701            max = mid-1;
    16811702        } else {
    16821703            min = mid;
    16831704        }
    1684     }
    1685 
     1705        mid = ((max+1)+min)/2;
     1706    }
     1707
     1708    psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1709            "---- Exiting VectorBinDisect(): bin %d\n", min);
    16861710    return(min);
    16871711}
    … …  
    17331757    float slope;
    17341758
     1759    psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1760            "---- Calling psSpline1DGen()\n");
     1761
    17351762    if (data->n != (1 + spline->n)) {
    17361763        psLogMsg(__func__, PS_LOG_WARN,
    … …  
    17441771            slope = (data->data.F32[i+1] - data->data.F32[i]) /
    17451772                    (spline->domains[i+1] - spline->domains[i]);
    1746             (spline->spline[0])->coeff[0] = data->data.F32[i] -
     1773            (spline->spline[i])->coeff[0] = data->data.F32[i] -
    17471774                                            (slope * spline->domains[i]);
    1748             (spline->spline[0])->coeff[1] = slope;
    1749         }
     1775            //printf("HMMM (%f, %f, %f)\n", data->data.F32[i], slope, spline->domains[i]);
     1776
     1777            (spline->spline[i])->coeff[1] = slope;
     1778            psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1779                    "---- spline %d coeffs are (%f, %f)\n", i,
     1780                    (spline->spline[i])->coeff[0],
     1781                    (spline->spline[i])->coeff[1]);
     1782        }
     1783        psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1784                "---- Exiting psSpline1DGen()\n");
    17501785        return((psSpline1D *) spline);
    17511786    } else if (4 == (spline->spline[0])->n) {

  • trunk/psLib/src/math/psSpline.c

    r1775 r1781  
    77 *  polynomials.  It also contains a Gaussian functions.
    88 *
    9  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    10  *  @date $Date: 2004-09-10 02:52:02 $
     9 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     10 *  @date $Date: 2004-09-10 22:25:59 $
    1111 *
    1212 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    441441    float xSum = 1.0;
    442442
     443    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     444            "---- Calling p_psOrdPolynomial1DEval(%f)\n", x);
     445    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     446            "Polynomial order is %d\n", myPoly->n);
     447    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
     448        psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
     449                "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]);
     450    }
     451
     452
    443453    if (NULL == myPoly) {
    444454        psAbort(__func__, "psPolynomial1DEval(): myPoly is NULL\n");
    … …  
    455465
    456466    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
    457         if (myPoly->mask[loop_x] != 0) {
     467        if (myPoly->mask[loop_x] == 0) {
    458468            polySum += xSum * myPoly->coeff[loop_x];
    459469            xSum *= x;
    … …  
    535545        ySum = xSum;
    536546        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    537             if (myPoly->mask[loop_x][loop_y] != 0) {
     547            if (myPoly->mask[loop_x][loop_y] == 0) {
    538548                polySum += ySum * myPoly->coeff[loop_x][loop_y];
    539549                ySum *= y;
    … …  
    564574    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
    565575        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    566             if (myPoly->mask[loop_x][loop_y] != 0) {
     576            if (myPoly->mask[loop_x][loop_y] == 0) {
    567577                polySum += myPoly->coeff[loop_x][loop_y] *
    568578                           psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    639649            zSum = ySum;
    640650            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    641                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     651                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    642652                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
    643653                    zSum *= z;
    … …  
    675685        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    676686            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    677                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     687                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    678688                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
    679689                               psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    768778                zSum = ySum;
    769779                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    770                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     780                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    771781                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
    772782                        zSum *= z;
    … …  
    811821            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    812822                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    813                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     823                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    814824                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
    815825                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
    … …  
    11291139
    11301140    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
    1131         if (myPoly->mask[loop_x] != 0) {
     1141        if (myPoly->mask[loop_x] == 0) {
    11321142            polySum += xSum * myPoly->coeff[loop_x];
    11331143            xSum *= x;
    … …  
    12101220        ySum = xSum;
    12111221        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    1212             if (myPoly->mask[loop_x][loop_y] != 0) {
     1222            if (myPoly->mask[loop_x][loop_y] == 0) {
    12131223                polySum += ySum * myPoly->coeff[loop_x][loop_y];
    12141224                ySum *= y;
    … …  
    12391249    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
    12401250        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    1241             if (myPoly->mask[loop_x][loop_y] != 0) {
     1251            if (myPoly->mask[loop_x][loop_y] == 0) {
    12421252                polySum += myPoly->coeff[loop_x][loop_y] *
    12431253                           psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    13131323            zSum = ySum;
    13141324            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1315                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1325                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    13161326                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
    13171327                    zSum *= z;
    … …  
    13491359        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    13501360            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1351                 if (myPoly->mask[loop_x][loop_y][loop_z] != 0) {
     1361                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
    13521362                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
    13531363                               psPolynomial1DEval(x, chebPolys[loop_x]) *
    … …  
    14441454                zSum = ySum;
    14451455                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1446                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1456                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    14471457                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
    14481458                        zSum *= z;
    … …  
    14871497            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
    14881498                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
    1489                     if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] != 0) {
     1499                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
    14901500                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
    14911501                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
    … …  
    16081618    tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
    16091619    width = (max - min) / ((float) numSplines);
    1610     tmpDomain = min;
    1611     for (i=0;i<numSplines+1;i++) {
     1620
     1621    (tmp->domains)[0] = min;
     1622    tmpDomain = min+width;
     1623    for (i=1;i<numSplines+1;i++) {
    16121624        (tmp->domains)[i] = tmpDomain;
    16131625        tmpDomain+= width;
    16141626    }
     1627    (tmp->domains)[numSplines] = max;
    16151628
    16161629    return(tmp);
    … …  
    16611674    int mid;
    16621675
     1676    psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1677            "---- Calling VectorBinDisect(%f)\n", x);
     1678
    16631679    if ((x < bins[0]) ||
    16641680            (x > bins[numBins-1])) {
    … …  
    16701686
    16711687    min = 0;
    1672     max = numBins-1;
    1673     mid = (max-min)/2;
    1674 
    1675     while (min+1 < max) {
    1676         psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4, "(min, max) is (%d, %d)\n", min, max);
    1677         mid = (max-min)/2;
    1678 
    1679         if (x < bins[mid]) {
    1680             max = mid;
     1688    max = numBins-2;
     1689    mid = ((max+1)-min)/2;
     1690
     1691    while (min != max) {
     1692        psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1693                "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
     1694                min, mid, max, x, bins[mid]);
     1695
     1696        if (x == bins[mid]) {
     1697            psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1698                    "---- Exiting VectorBinDisect(): bin %d\n", min);
     1699            return(min);
     1700        } else if (x < bins[mid]) {
     1701            max = mid-1;
    16811702        } else {
    16821703            min = mid;
    16831704        }
    1684     }
    1685 
     1705        mid = ((max+1)+min)/2;
     1706    }
     1707
     1708    psTrace(".psLib.dataManip.psFunctions.VectorBinDisect", 4,
     1709            "---- Exiting VectorBinDisect(): bin %d\n", min);
    16861710    return(min);
    16871711}
    … …  
    17331757    float slope;
    17341758
     1759    psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1760            "---- Calling psSpline1DGen()\n");
     1761
    17351762    if (data->n != (1 + spline->n)) {
    17361763        psLogMsg(__func__, PS_LOG_WARN,
    … …  
    17441771            slope = (data->data.F32[i+1] - data->data.F32[i]) /
    17451772                    (spline->domains[i+1] - spline->domains[i]);
    1746             (spline->spline[0])->coeff[0] = data->data.F32[i] -
     1773            (spline->spline[i])->coeff[0] = data->data.F32[i] -
    17471774                                            (slope * spline->domains[i]);
    1748             (spline->spline[0])->coeff[1] = slope;
    1749         }
     1775            //printf("HMMM (%f, %f, %f)\n", data->data.F32[i], slope, spline->domains[i]);
     1776
     1777            (spline->spline[i])->coeff[1] = slope;
     1778            psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1779                    "---- spline %d coeffs are (%f, %f)\n", i,
     1780                    (spline->spline[i])->coeff[0],
     1781                    (spline->spline[i])->coeff[1]);
     1782        }
     1783        psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4,
     1784                "---- Exiting psSpline1DGen()\n");
    17501785        return((psSpline1D *) spline);
    17511786    } else if (4 == (spline->spline[0])->n) {

  • trunk/psLib/test/dataManip/tst_psFunc04.c

    r1775 r1781  
    1111#include "psFunctions.h"
    1212
    13 #define N 10
     13#define N 8
    1414
    1515int main()
    … …  
    2323    psVector *data;
    2424
    25     psTraceSetLevel(".", 9);
     25    psTraceSetLevel(".", 0);
    2626    /****************************************************************************/
    2727    /****************************************************************************/
    … …  
    3333    tmpSpline = psSpline1DAlloc(N, 1, 1.0, 10.0);
    3434
     35
    3536    for (i=0;i<N+1;i++) {
    36         data->data.F32[i] = (float) tmpSpline->domains[i];
     37        printf("tmpSpline->domains[%d] is %f\n", i, tmpSpline->domains[i]);
     38    }
     39
     40    for (i=0;i<N+1;i++) {
     41        data->data.F32[i] = tmpSpline->domains[i];
    3742    }
    3843
Note: See TracChangeset for help on using the changeset viewer.