IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1014


Ignore:
Timestamp:
Jun 12, 2004, 9:46:04 AM (22 years ago)
Author:
desonia
Message:

added additional tests.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/test/image/tst_psImageFFT.c

    r1011 r1014  
    66 *  @author Robert DeSonia, MHPCC
    77 *
    8  *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-06-12 02:17:42 $
     8 *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-06-12 19:46:04 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2828
    2929static int testImageFFT(void);
    30 #if 0
    3130static int testImageRealImaginary(void);
    3231static int testImageComplex(void);
    3332static int testImageConjugate(void);
    3433static int testImagePowerSpectrum(void);
    35 #endif
    3634
    3735testDescription tests[] = {
    3836                              {testImageFFT,"600-testImageFFT",0},
    39                               #if 0
    4037                              {testImageRealImaginary,"601-testImageRealImaginary",0},
    4138                              {testImageComplex,"602-testImageComplex",0},
    4239                              {testImageConjugate,"603-testImageConjugate",0},
    4340                              {testImagePowerSpectrum,"604-testImagePowerSpectrum",0},
    44                               #endif
    4541                              {NULL}
    4642                          };
    … …  
    7369
    7470    // 1. assign a image to a radial sinisoid
    75     //    GENIMAGE(img,m,n,F32, sinf(sqrt((n/2-row)*(n/2-row))/n*2*M_PI));
    76     GENIMAGE(img,m,n,F32, 1);
    77     psImageWriteSection(img,0,0,0,NULL,0,"test.fits");
     71    GENIMAGE(img,m,n,F32, sinf((32.0f-row)/16.0f*M_PI)+sinf((64.0f-col)/32.0f*M_PI));
    7872
    7973    // 2. perform a forward transform
    … …  
    8478    }
    8579
    86     img3 = psImageReal(img3,img2);
    87     psImageWriteSection(img2,0,0,0,NULL,1,"test.fits");
    88 
    89     #if 0
    9080    // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    91     for (unsigned int n = 0; n<100; n++) {
    92         if (n==1 || n==99) {
    93             if (fabsf(cabsf(img2->data.C32[n]) - 50.0f) > 0.1f) {
    94                 psError(__func__,"FFT didn't work for vector (n=%d)",n);
    95                 return 2;
    96             }
    97         } else {
    98             if (fabsf(cabsf(img2->data.C32[n])) > 0.1f) {
    99                 psError(__func__,"FFT didn't work for vector (n=%d)",n);
    100                 return 3;
     81    for (unsigned int row=0;row<n;row++) {
     82        psC32* img2Row = img2->data.C32[row];
     83        for (unsigned int col=0;col<m;col++) {
     84            psF32 mag = cabsf(img2Row[col])/m/n;
     85            if (mag > 0.1f) {
     86                // must be (0,1) or (0,n-1) or (1,0) or (n-1,0)
     87                if (! (col == 0 && (row == 1 || row == n-1))
     88                        && ! (row == 0 && (col==1 || col == n-1)) ) {
     89                    psError(__func__,"Result invalid at %d,%d (%.2f)",col,row,mag);
     90                }
     91            } else if ( (col == 0 && (row == 1 || row == n-1))
     92                        || (row == 0 && (col==1 || col == n-1)) ) {
     93                psError(__func__,"Result invalid at %d,%d (%.2f)",col,row,mag);
    10194            }
    10295        }
    … …  
    109102        return 4;
    110103    }
    111     for (unsigned int n = 0; n<100; n++) {
    112         psF32 val =sinf((psF32)n / 50.0f * M_PI);
    113         psF32 vecVal = crealf(img3->data.C32[n])/100;
    114         if (fabsf(vecVal - val) > 0.1f) {
    115             psError(__func__,"Reverse FFT didn't give me the original vector back (n=%d) (%.2f vs %.2f)",
    116                     n,vecVal,val);
    117             return 5;
    118         }
    119     }
    120     #endif
     104
     105    for (unsigned int row=0;row<n;row++) {
     106        psC32* img3Row = img3->data.C32[row];
     107        psF32* imgRow = img->data.F32[row];
     108        for (unsigned int col=0;col<m;col++) {
     109            psF32 pixel = creal(img3Row[col])/m/n;
     110            if (fabsf(pixel-imgRow[col]) > 0.1) {
     111                psError(__func__,"Reverse FFT didn't gime original image back (%d,%d %.2f vs %.2f)",
     112                        col,row,pixel,imgRow[col]);
     113                return 5;
     114            }
     115        }
     116    }
     117
    121118    psImageFree(img);
    122119    psImageFree(img2);
    … …  
    125122    return 0;
    126123}
    127 #if 0
     124
    128125int testImageRealImaginary(void)
    129126{
    … …  
    131128    psImage* img2 = NULL;
    132129    psImage* img3 = NULL;
     130    unsigned int m = 128;
     131    unsigned int n = 64;
    133132
    134133    /*
    … …  
    139138
    140139    // 1. create a C32 complex vector with distinctly different real and imaginary parts.
    141     img=psImageAlloc(100,PS_TYPE_C32);
    142     img->n = img->nalloc;
    143     for (unsigned int n = 0; n<100; n++) {
    144         img->data.C32[n] = n + I * (n*2);
    145     }
     140    GENIMAGE(img,m,n,C32, row + I * col);
    146141
    147142    // 2. call psImageReal and psImageImaginary
    … …  
    169164
    170165    // 3. compare results to the real/imaginary components of input
    171     for (unsigned int n = 0; n<100; n++) {
    172         psF32 r = n;
    173         psF32 i = (n*2);
    174         if (fabsf(img2->data.F32[n] -r) > FLT_EPSILON) {
    175             psError(__func__,"psImageReal didn't return the real portion at n=%d",
    176                     n);
    177             return 5;
    178         }
    179         if (fabsf(img3->data.F32[n] -i) > FLT_EPSILON) {
    180             psError(__func__,"psImageImaginary didn't return the real portion at n=%d",
    181                     n);
    182             return 6;
     166    for (unsigned int row=0;row<n;row++) {
     167        psF32* img2Row = img2->data.F32[row];
     168        psF32* img3Row = img3->data.F32[row];
     169        for (unsigned int col=0;col<m;col++) {
     170            if (fabsf(img2Row[col] - row) > FLT_EPSILON) {
     171                psError(__func__,"psImageReal didn't return the real portion at n=%d",
     172                        n);
     173                return 5;
     174            }
     175            if (fabsf(img3Row[col] - col) > FLT_EPSILON) {
     176                psError(__func__,"psImageImaginary didn't return the imag portion at n=%d",
     177                        n);
     178                return 6;
     179            }
    183180        }
    184181    }
    … …  
    196193    psImage* img2 = NULL;
    197194    psImage* img3 = NULL;
     195    unsigned int m = 128;
     196    unsigned int n = 64;
    198197
    199198    /*
    … …  
    214213
    215214    // 1. create two unique psF32 vectors of the same size
    216     img=psImageAlloc(100,PS_TYPE_F32);
    217     img2=psImageAlloc(100,PS_TYPE_F32);
    218     img->n = img->nalloc;
    219     img2->n = img2->nalloc;
    220     for (unsigned int n = 0; n<100; n++) {
    221         img->data.F32[n] = n;
    222         img2->data.F32[n] = (n*2);
    223     }
     215    GENIMAGE(img,m,n,F32,row);
     216    GENIMAGE(img2,m,n,F32,col);
    224217
    225218    // 2. call psImageComplex
    … …  
    228221    // 3. verify that the result is a psC32
    229222    if (img3->type.type != PS_TYPE_C32) {
    230         psError(__func__,"Vector Type from psImageComplex is not complex? (%d)",
     223        psError(__func__,"Image Type from psImageComplex is not complex? (%d)",
    231224                img3->type.type);
    232225        return 1;
    … …  
    235228    // 4. call psImageReal and psImageImaginary on step 2 results (not needed, just use crealf/cimagf)
    236229    // 5. compare step 4 results to input.
    237     for (unsigned int n = 0; n<100; n++) {
    238         if (fabsf(crealf(img3->data.C32[n]) - n) > FLT_EPSILON ||
    239                 fabsf(cimagf(img3->data.C32[n]) - (n*2)) > FLT_EPSILON) {
    240             psError(__func__,"psImageComplex result is invalid (n=%d, %.2f+%.2fi)",
    241                     n,crealf(img3->data.C32[n]),cimagf(img3->data.C32[n]));
    242             return 2;
    243         };
    244     }
    245 
    246 
    247     // 6. create a psF32 and a psF64 vector of the same size
    248     img2 = psImageRecycle(img2,PS_TYPE_F64, 100);
     230    for (unsigned int row=0;row<n;row++) {
     231        psC32* img3Row = img3->data.C32[row];
     232        for (unsigned int col=0;col<m;col++) {
     233            if (fabsf(crealf(img3Row[col]) - row) > FLT_EPSILON ||
     234                    fabsf(cimagf(img3Row[col]) - col) > FLT_EPSILON) {
     235                psError(__func__,"psImageComplex result is invalid (%d,%d, %.2f+%.2fi)",
     236                        col,row,crealf(img3Row[col]),cimagf(img3Row[col]));
     237                return 2;
     238            }
     239        }
     240    }
     241
     242    // 6. create a psF32 and a psF64 image of the same size
     243    img2 = psImageRecycle(img2,m,n,PS_TYPE_F64);
    249244
    250245    // 7. call psImageComplex
    … …  
    253248    // 8. verify that an appropriate error occurred. (this partially has to be done via inspection)
    254249    if (img3 != NULL) {
    255         psError(__func__,"psImageComplex returned a vector though input types mismatched.");
     250        psError(__func__,"psImageComplex returned a image though input types mismatched.");
    256251        return 3;
    257252    }
    258253
    259254    // 9. create two psf32 vectors of different sizes
    260     img2 = psImageRecycle(img2,PS_TYPE_F32,200);
     255    img2 = psImageRecycle(img2,m/2,n,PS_TYPE_F32);
    261256
    262257    // 10. call psImageComplex
     258    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error (size mismatch).");
    263259    img3 = psImageComplex(img3,img,img2);
    264260
    265     // 11. verify thet an appropriate error occurred. (actually, it isn't an error...)
    266     if (img3->n != 100) {
    267         psError(__func__,"psImageComplex returned a vector though input sizes mismatched.");
     261    // 11. verify thet an appropriate error occurred.
     262    if (img3 != NULL) {
     263        psError(__func__,"psImageComplex returned a image though input sizes mismatched.");
    268264        return 4;
    269265    }
    … …  
    280276    psImage* img = NULL;
    281277    psImage* img2 = NULL;
     278    unsigned int m = 128;
     279    unsigned int n = 64;
    282280
    283281    /*
    … …  
    289287
    290288    // 1. create a psC32 with unique real and imaginary values.
    291     img=psImageAlloc(100,PS_TYPE_C32);
    292     img->n = img->nalloc;
    293     for (unsigned int n = 0; n<100; n++) {
    294         img->data.C32[n] = n + I * (n*2);
    295     }
     289    GENIMAGE(img,m,n,C32, row + I * col);
    296290
    297291    // 2. call psImageConjugate
    … …  
    300294    // 3. verify result is psC32
    301295    if (img2->type.type != PS_TYPE_C32) {
    302         psError(__func__,"the psImageConjugate didn't return a C32 vector");
     296        psError(__func__,"the psImageConjugate didn't return a C32 image.");
    303297        return 1;
    304298    }
    305299
    306300    // 4. verify each value is conjugate of input (a+bi -> a-bi)
    307     for (unsigned int n = 0; n<100; n++) {
    308         if (fabsf(crealf(img->data.C32[n]) - crealf(img2->data.C32[n])) > FLT_EPSILON ||
    309                 fabsf(cimagf(img->data.C32[n]) + cimagf(img2->data.C32[n])) > FLT_EPSILON) {
    310             psError(__func__,"psImageComplex result is invalid (n=%d, %.2f+%.2fi)",
    311                     n,crealf(img2->data.C32[n]),cimagf(img2->data.C32[n]));
    312             return 2;
    313         };
     301    for (unsigned int row=0;row<n;row++) {
     302        psC32* img2Row = img2->data.C32[row];
     303        for (unsigned int col=0;col<m;col++) {
     304            if (fabsf(crealf(img2Row[col]) - row) > FLT_EPSILON ||
     305                    fabsf(cimagf(img2Row[col]) + col) > FLT_EPSILON) {
     306                psError(__func__,"psImageComplex result is invalid (%d,%d, %.2f+%.2fi)",
     307                        col,row,crealf(img2Row[col]),cimagf(img2Row[col]));
     308                return 2;
     309            }
     310        }
    314311    }
    315312
    … …  
    324321    psImage* img = NULL;
    325322    psImage* img2 = NULL;
    326     psF32 val;
     323    unsigned int m = 128;
     324    unsigned int n = 64;
    327325
    328326    /*
    … …  
    334332
    335333    // 1. create a psC32 vector with unique real and imaginary components
    336     img=psImageAlloc(100,PS_TYPE_C32);
    337     img->n = img->nalloc;
    338     for (unsigned int n = 0; n<100; n++) {
    339         img->data.C32[n] = n + I * sinf(((psF32)n) / 50.f * M_PI);
    340     }
     334    GENIMAGE(img,m,n,C32, row + I * col);
    341335
    342336    // 2. call psImagePowerSpectrum
    … …  
    350344
    351345    // 4. verify the values are the square of the absolute values of the original
    352     //   (ADD specifies something else)
    353     //   P_0 = |C_0|^2/N^2
    354     //   P_j = (|C_j|^2+|C_N-j|^2)/N^2
    355     //   P_N/2 = |C_N/2|^2/N^2
    356     //  where j = 1,2,...,(N/2-1)
    357 
    358     val = cabsf(img->data.C32[0])*cabsf(img->data.C32[0])/100/100;
    359     if (fabsf(img2->data.F32[0] - val) > FLT_EPSILON) {
    360         psError(__func__,"psImagePowerSpectrum result is invalid (n=0, %.2f %.2f)",
    361                 img2->data.F32[0],val);
    362         return 2;
    363     };
    364 
    365     for (unsigned int n = 1; n<50; n++) {
    366         val = ( cabsf(img->data.C32[n])*cabsf(img->data.C32[n])+
    367                 cabsf(img->data.C32[100-n])*cabsf(img->data.C32[100-n]) ) /100/100;
    368 
    369         if (fabsf(val - img2->data.F32[n]) > FLT_EPSILON) {
    370             psError(__func__,"psImagePowerSpectrum result is invalid (n=%d, %.2f %.2f)",
    371                     n,img2->data.F32[n],val);
    372             return 2;
    373         };
    374     }
    375 
    376     val = cabsf(img->data.C32[50])*cabsf(img->data.C32[50])/100/100;
    377     if (fabsf(img2->data.F32[50] - val) > FLT_EPSILON) {
    378         psError(__func__,"psImagePowerSpectrum result is invalid (n=50, %.2f %.2f)",
    379                 img2->data.F32[0],val);
    380         return 2;
    381     };
    382 
    383     psImageFree(img);
    384     psImageFree(img2);
    385 
    386     return 0;
    387 }
    388 #endif
     346    for (unsigned int row=0;row<n;row++) {
     347        psC32* imgRow = img->data.C32[row];
     348        psF32* img2Row = img2->data.F32[row];
     349        for (unsigned int col=0;col<m;col++) {
     350            psF32 power = cabs(imgRow[col]);
     351            power *= power/n/n/m/m;
     352
     353            if (fabsf(img2Row[col] - power) > 2.0f*FLT_EPSILON) {
     354                psError(__func__,"psImagePowerSpectrum result is invalid (%d,%d, %.2f vs %.2f)",
     355                        col,row,img2Row[col],power);
     356                return 2;
     357            }
     358        }
     359    }
     360
     361    psImageFree(img);
     362    psImageFree(img2);
     363
     364    return 0;
     365}
Note: See TracChangeset for help on using the changeset viewer.