IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2605


Ignore:
Timestamp:
Dec 3, 2004, 12:29:43 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/image/psImageStats.c

    r2602 r2605  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-03 19:43:43 $
     11 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-03 22:29:43 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    174174    return (chebPolys);
    175175}
     176
     177/*****************************************************************************
     178psImageFitPolynomial(): This routine takes as input a 2-D image and produces
     179as output the coefficients of the Chebyshev polynomials which match that
     180input image.
     181  Input:
     182  Output:
     183  Internal Data Structures:
     184    chebPolys[i][j]
     185    sums[i][j]: This will contain the sum of
     186                input->data.F32[x][y] *
     187                psPolynomial1DEval(
     188chebPolys[i],
     189(float) x) *
     190                psPolynomial1DEval(
     191chebPolys[j],
     192(float) y,
     193);
     194        over all pixels (x,y) in the image.
     195  *****************************************************************************/
     196psPolynomial2D* psImageFitPolynomialORIG(psPolynomial2D* coeffs,
     197        const psImage* input)
     198{
     199    PS_IMAGE_CHECK_NULL(input, NULL);
     200    PS_IMAGE_CHECK_EMPTY(input, NULL);
     201    if ((input->type.type != PS_TYPE_S8) &&
     202            (input->type.type != PS_TYPE_U16) &&
     203            (input->type.type != PS_TYPE_F32) &&
     204            (input->type.type != PS_TYPE_F64)) {
     205        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n");
     206    }
     207    PS_POLY_CHECK_NULL(coeffs, NULL);
     208    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
     209    psS32 x = 0;
     210    psS32 y = 0;
     211    psS32 i = 0;
     212    psS32 j = 0;
     213    double **sums = NULL;
     214    psPolynomial1D* *chebPolys = NULL;
     215    psS32 maxChebyPoly = 0;
     216    double *cScalingFactors = NULL;
     217    double *rScalingFactors = NULL;
     218
     219    // Create the sums[][] data structure.  This
     220    // will hold the LHS of
     221    // equation
     222    // 29 in the ADD: sums[k][l] = SUM {
     223    // image(x,y) * Tk(x) * Tl(y) }
     224    sums = (double **)psAlloc(coeffs->nX * sizeof(double *));
     225    for (i = 0; i < coeffs->nX; i++) {
     226        sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double));
     227    }
     228    // We scale the pixel positions to values
     229    // between -1.0 and 1.0
     230    rScalingFactors = calcScaleFactors(input->numRows);
     231    cScalingFactors = calcScaleFactors(input->numCols);
     232
     233    // Determine how many Chebyshev polynomials
     234    // are needed, then create them.
     235    maxChebyPoly = coeffs->nX;
     236    if (coeffs->nY > coeffs->nX) {
     237        maxChebyPoly = coeffs->nY;
     238    }
     239    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
     240
     241    // Compute the sums[][] data structure.
     242    for (i = 0; i < coeffs->nX; i++) {
     243        for (j = 0; j < coeffs->nY; j++) {
     244            sums[i][j] = 0.0;
     245            for (x = 0; x < input->numRows; x++) {
     246                for (y = 0; y < input->numCols; y++) {
     247                    double pixel;
     248                    if (input->type.type == PS_TYPE_S8) {
     249                        pixel = (double) input->data.S8[x][y];
     250                    } else if (input->type.type == PS_TYPE_U16) {
     251                        pixel = (double) input->data.U16[x][y];
     252                    } else if (input->type.type == PS_TYPE_F32) {
     253                        pixel = (double) input->data.F32[x][y];
     254                    } else if (input->type.type == PS_TYPE_F64) {
     255                        pixel = input->data.F64[x][y];
     256                    }
     257                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) *
     258                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
     259                }
     260            }
     261        }
     262    }
     263
     264    for (i = 0; i < coeffs->nX; i++) {
     265        for (j = 0; j < coeffs->nY; j++) {
     266            coeffs->coeff[i][j] = sums[i][j];
     267            coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols);
     268
     269            if ((i != 0) && (j != 0)) {
     270                coeffs->coeff[i][j] *= 4.0;
     271            } else if ((i == 0) && (j == 0)) {
     272                coeffs->coeff[i][j] *= 1.0;
     273            } else {
     274                coeffs->coeff[i][j] *= 2.0;
     275            }
     276        }
     277    }
     278
     279    // Free the Chebyshev polynomials that were
     280    // created in this routine.
     281    for (i = 0; i < maxChebyPoly; i++) {
     282        psFree(chebPolys[i]);
     283    }
     284    psFree(chebPolys);
     285
     286    // Free some data
     287    for (i = 0; i < coeffs->nX; i++) {
     288        psFree(sums[i]);
     289    }
     290    psFree(sums);
     291    psFree(cScalingFactors);
     292    psFree(rScalingFactors);
     293
     294    return (coeffs);
     295}
     296
     297
     298
    176299
    177300/*****************************************************************************
    … …  
    216339    double *cScalingFactors = NULL;
    217340    double *rScalingFactors = NULL;
     341    psImage *nodes = psImageAlloc(input->numCols, input->numRows, PS_TYPE_F64);
     342
     343    double min = -1.0;
     344    double max = 1.0;
     345    double bma = 0.5 * (max-min);  // 1
     346    double bpa = 0.5 * (max+min);  // 0
     347    // We must calculate the value of the image at the nodes where the
     348    // Chebyshev polynomials are 0.
     349    for (x = 0; x < input->numRows; x++) {
     350        double xTmp = cos(M_PI * (0.5 + ((float) x)) / ((float) input->numRows));
     351        double xNode = - ((xTmp + bma + bpa) - 1.0);
     352        double xOrig = ((float) input->numRows) * (xNode - min) / (max - min);
     353
     354        for (y = 0; y < input->numCols; y++) {
     355            double yTmp = cos(M_PI * (0.5 + ((float) y)) / ((float) input->numCols));
     356            double yNode = - ((yTmp + bma + bpa) - 1.0);
     357            double yOrig = ((float) input->numCols) * (yNode - min) / (max - min);
     358
     359            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yNode, xNode, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     360            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yTmp, xTmp, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     361            printf("HMMMM: (xOrig, yOrig) is (%f, %f)\n", xOrig, yOrig);
     362            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yOrig, xOrig, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     363            printf("HMMMM: interpolated pixel [%d][%d] is %f, should be %f\n", x, y, nodes->data.F64[x][y], input->data.F32[x][y]);
     364        }
     365    }
    218366
    219367    // Create the sums[][] data structure.  This
    … …  
    246394                for (y = 0; y < input->numCols; y++) {
    247395                    double pixel;
    248                     if (input->type.type == PS_TYPE_S8) {
    249                         pixel = (double) input->data.S8[x][y];
    250                     } else if (input->type.type == PS_TYPE_U16) {
    251                         pixel = (double) input->data.U16[x][y];
    252                     } else if (input->type.type == PS_TYPE_F32) {
    253                         pixel = (double) input->data.F32[x][y];
    254                     } else if (input->type.type == PS_TYPE_F64) {
    255                         pixel = input->data.F64[x][y];
    256                     }
     396                    /*
     397                                        if (input->type.type == PS_TYPE_S8) {
     398                                            pixel = (double) input->data.S8[x][y];
     399                                        } else if (input->type.type == PS_TYPE_U16) {
     400                                            pixel = (double) input->data.U16[x][y];
     401                                        } else if (input->type.type == PS_TYPE_F32) {
     402                                            pixel = (double) input->data.F32[x][y];
     403                                        } else if (input->type.type == PS_TYPE_F64) {
     404                                            pixel = input->data.F64[x][y];
     405                                        }
     406                    */
     407                    pixel = nodes->data.F64[x][y];
    257408                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) *
    258409                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
    … …  
    291442    psFree(cScalingFactors);
    292443    psFree(rScalingFactors);
     444    psFree(nodes);
    293445
    294446    return (coeffs);
    … …  
    342494                for (j = 0; j < coeffs->nY; j++) {
    343495                    polySum +=
    344                         psPolynomial1DEval(
    345                             chebPolys[i],
    346                             rScalingFactors[x]
    347                         ) *
    348                         psPolynomial1DEval(
    349                             chebPolys[j],
    350                             cScalingFactors[y]
    351 
    352                         ) *
     496                        psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *
     497                        psPolynomial1DEval(chebPolys[j], cScalingFactors[y]) *
    353498                        coeffs->coeff[i][j];
    354 
    355499                }
    356500            }

  • trunk/psLib/src/imageops/psImageStats.c

    r2602 r2605  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-12-03 19:43:43 $
     11 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-12-03 22:29:43 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    174174    return (chebPolys);
    175175}
     176
     177/*****************************************************************************
     178psImageFitPolynomial(): This routine takes as input a 2-D image and produces
     179as output the coefficients of the Chebyshev polynomials which match that
     180input image.
     181  Input:
     182  Output:
     183  Internal Data Structures:
     184    chebPolys[i][j]
     185    sums[i][j]: This will contain the sum of
     186                input->data.F32[x][y] *
     187                psPolynomial1DEval(
     188chebPolys[i],
     189(float) x) *
     190                psPolynomial1DEval(
     191chebPolys[j],
     192(float) y,
     193);
     194        over all pixels (x,y) in the image.
     195  *****************************************************************************/
     196psPolynomial2D* psImageFitPolynomialORIG(psPolynomial2D* coeffs,
     197        const psImage* input)
     198{
     199    PS_IMAGE_CHECK_NULL(input, NULL);
     200    PS_IMAGE_CHECK_EMPTY(input, NULL);
     201    if ((input->type.type != PS_TYPE_S8) &&
     202            (input->type.type != PS_TYPE_U16) &&
     203            (input->type.type != PS_TYPE_F32) &&
     204            (input->type.type != PS_TYPE_F64)) {
     205        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unallowable image type.\n");
     206    }
     207    PS_POLY_CHECK_NULL(coeffs, NULL);
     208    PS_POLY_CHECK_TYPE(coeffs, PS_POLYNOMIAL_CHEB, NULL);
     209    psS32 x = 0;
     210    psS32 y = 0;
     211    psS32 i = 0;
     212    psS32 j = 0;
     213    double **sums = NULL;
     214    psPolynomial1D* *chebPolys = NULL;
     215    psS32 maxChebyPoly = 0;
     216    double *cScalingFactors = NULL;
     217    double *rScalingFactors = NULL;
     218
     219    // Create the sums[][] data structure.  This
     220    // will hold the LHS of
     221    // equation
     222    // 29 in the ADD: sums[k][l] = SUM {
     223    // image(x,y) * Tk(x) * Tl(y) }
     224    sums = (double **)psAlloc(coeffs->nX * sizeof(double *));
     225    for (i = 0; i < coeffs->nX; i++) {
     226        sums[i] = (double *)psAlloc(coeffs->nY * sizeof(double));
     227    }
     228    // We scale the pixel positions to values
     229    // between -1.0 and 1.0
     230    rScalingFactors = calcScaleFactors(input->numRows);
     231    cScalingFactors = calcScaleFactors(input->numCols);
     232
     233    // Determine how many Chebyshev polynomials
     234    // are needed, then create them.
     235    maxChebyPoly = coeffs->nX;
     236    if (coeffs->nY > coeffs->nX) {
     237        maxChebyPoly = coeffs->nY;
     238    }
     239    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly);
     240
     241    // Compute the sums[][] data structure.
     242    for (i = 0; i < coeffs->nX; i++) {
     243        for (j = 0; j < coeffs->nY; j++) {
     244            sums[i][j] = 0.0;
     245            for (x = 0; x < input->numRows; x++) {
     246                for (y = 0; y < input->numCols; y++) {
     247                    double pixel;
     248                    if (input->type.type == PS_TYPE_S8) {
     249                        pixel = (double) input->data.S8[x][y];
     250                    } else if (input->type.type == PS_TYPE_U16) {
     251                        pixel = (double) input->data.U16[x][y];
     252                    } else if (input->type.type == PS_TYPE_F32) {
     253                        pixel = (double) input->data.F32[x][y];
     254                    } else if (input->type.type == PS_TYPE_F64) {
     255                        pixel = input->data.F64[x][y];
     256                    }
     257                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) *
     258                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
     259                }
     260            }
     261        }
     262    }
     263
     264    for (i = 0; i < coeffs->nX; i++) {
     265        for (j = 0; j < coeffs->nY; j++) {
     266            coeffs->coeff[i][j] = sums[i][j];
     267            coeffs->coeff[i][j] /= (double)(input->numRows * input->numCols);
     268
     269            if ((i != 0) && (j != 0)) {
     270                coeffs->coeff[i][j] *= 4.0;
     271            } else if ((i == 0) && (j == 0)) {
     272                coeffs->coeff[i][j] *= 1.0;
     273            } else {
     274                coeffs->coeff[i][j] *= 2.0;
     275            }
     276        }
     277    }
     278
     279    // Free the Chebyshev polynomials that were
     280    // created in this routine.
     281    for (i = 0; i < maxChebyPoly; i++) {
     282        psFree(chebPolys[i]);
     283    }
     284    psFree(chebPolys);
     285
     286    // Free some data
     287    for (i = 0; i < coeffs->nX; i++) {
     288        psFree(sums[i]);
     289    }
     290    psFree(sums);
     291    psFree(cScalingFactors);
     292    psFree(rScalingFactors);
     293
     294    return (coeffs);
     295}
     296
     297
     298
    176299
    177300/*****************************************************************************
    … …  
    216339    double *cScalingFactors = NULL;
    217340    double *rScalingFactors = NULL;
     341    psImage *nodes = psImageAlloc(input->numCols, input->numRows, PS_TYPE_F64);
     342
     343    double min = -1.0;
     344    double max = 1.0;
     345    double bma = 0.5 * (max-min);  // 1
     346    double bpa = 0.5 * (max+min);  // 0
     347    // We must calculate the value of the image at the nodes where the
     348    // Chebyshev polynomials are 0.
     349    for (x = 0; x < input->numRows; x++) {
     350        double xTmp = cos(M_PI * (0.5 + ((float) x)) / ((float) input->numRows));
     351        double xNode = - ((xTmp + bma + bpa) - 1.0);
     352        double xOrig = ((float) input->numRows) * (xNode - min) / (max - min);
     353
     354        for (y = 0; y < input->numCols; y++) {
     355            double yTmp = cos(M_PI * (0.5 + ((float) y)) / ((float) input->numCols));
     356            double yNode = - ((yTmp + bma + bpa) - 1.0);
     357            double yOrig = ((float) input->numCols) * (yNode - min) / (max - min);
     358
     359            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yNode, xNode, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     360            //            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yTmp, xTmp, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     361            printf("HMMMM: (xOrig, yOrig) is (%f, %f)\n", xOrig, yOrig);
     362            nodes->data.F64[x][y] = psImagePixelInterpolate(input, yOrig, xOrig, NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
     363            printf("HMMMM: interpolated pixel [%d][%d] is %f, should be %f\n", x, y, nodes->data.F64[x][y], input->data.F32[x][y]);
     364        }
     365    }
    218366
    219367    // Create the sums[][] data structure.  This
    … …  
    246394                for (y = 0; y < input->numCols; y++) {
    247395                    double pixel;
    248                     if (input->type.type == PS_TYPE_S8) {
    249                         pixel = (double) input->data.S8[x][y];
    250                     } else if (input->type.type == PS_TYPE_U16) {
    251                         pixel = (double) input->data.U16[x][y];
    252                     } else if (input->type.type == PS_TYPE_F32) {
    253                         pixel = (double) input->data.F32[x][y];
    254                     } else if (input->type.type == PS_TYPE_F64) {
    255                         pixel = input->data.F64[x][y];
    256                     }
     396                    /*
     397                                        if (input->type.type == PS_TYPE_S8) {
     398                                            pixel = (double) input->data.S8[x][y];
     399                                        } else if (input->type.type == PS_TYPE_U16) {
     400                                            pixel = (double) input->data.U16[x][y];
     401                                        } else if (input->type.type == PS_TYPE_F32) {
     402                                            pixel = (double) input->data.F32[x][y];
     403                                        } else if (input->type.type == PS_TYPE_F64) {
     404                                            pixel = input->data.F64[x][y];
     405                                        }
     406                    */
     407                    pixel = nodes->data.F64[x][y];
    257408                    sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i],rScalingFactors[x]) *
    258409                                  psPolynomial1DEval(chebPolys[j], cScalingFactors[y]);
    … …  
    291442    psFree(cScalingFactors);
    292443    psFree(rScalingFactors);
     444    psFree(nodes);
    293445
    294446    return (coeffs);
    … …  
    342494                for (j = 0; j < coeffs->nY; j++) {
    343495                    polySum +=
    344                         psPolynomial1DEval(
    345                             chebPolys[i],
    346                             rScalingFactors[x]
    347                         ) *
    348                         psPolynomial1DEval(
    349                             chebPolys[j],
    350                             cScalingFactors[y]
    351 
    352                         ) *
     496                        psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *
     497                        psPolynomial1DEval(chebPolys[j], cScalingFactors[y]) *
    353498                        coeffs->coeff[i][j];
    354 
    355499                }
    356500            }
Note: See TracChangeset for help on using the changeset viewer.