IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 23, 2004, 8:31:49 AM (22 years ago)
Author:
desonia
Message:

migrated more psError calls to new style (i.e., psErrorMsg).

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/fft/psVectorFFT.c

    r1625 r1864  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2004-08-25 21:10:08 $
     7 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2004-09-23 18:31:49 $
    99 *
    1010 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2323#include "psImageExtraction.h"
    2424
     25#include "psDataManipErrors.h"
     26
    2527#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
    2628
    … …  
    4244
    4345    if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {
    44         psError(__func__, "Input image must be a 32-bit float or complex image (type=%d)", type);
     46        char* typeStr;
     47        PS_TYPE_NAME(typeStr,type);
     48        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
     49                   PS_ERR_BAD_PARAMETER_TYPE, true,
     50                   PS_ERRORTEXT_psVectorFFT_IMAGE_TYPE_UNSUPPORTED,
     51                   typeStr);
    4552        psFree(out);
    4653        return NULL;
    … …  
    4855
    4956    if ((type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE)) {
    50         psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
     57        char* typeStr;
     58        PS_TYPE_NAME(typeStr,type);
     59        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
     60                   PS_ERR_BAD_PARAMETER_TYPE, true,
     61                   PS_ERRORTEXT_psVectorFFT_REVERSE_NOT_COMPLEX,
     62                   typeStr);
    5163        psFree(out);
    5264        return NULL;
    … …  
    5567
    5668    if ((type != PS_TYPE_F32) && (direction == PS_FFT_FORWARD)) {
    57         psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
     69        char* typeStr;
     70        PS_TYPE_NAME(typeStr,type);
     71        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
     72                   PS_ERR_BAD_PARAMETER_TYPE, true,
     73                   PS_ERRORTEXT_psVectorFFT_FORWARD_NOT_REAL,
     74                   typeStr);
    5875        psFree(out);
    5976        return NULL;
    … …  
    94111    /* check if a plan exists now */
    95112    if (plan == NULL) {
    96         psError(__func__, "Failed to create FFTW plan.");
     113        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorFFT",
     114                   PS_ERR_BAD_PARAMETER_TYPE, true,
     115                   PS_ERRORTEXT_psVectorFFT_FFTW_PLAN_NULL);
    97116        psFree(out);
    98117        return NULL;
    … …  
    142161            outVec[i] = crealf(inVec[i]);
    143162        }
    144     } else {
    145         psError(__func__, "Can not extract real component from given vector type (%d).", type);
     163    } else if (type == PS_TYPE_C64) {
     164        psF64* outVec;
     165        psC64* inVec = in->data.C64;
     166
     167        out = psVectorRecycle(out, numElements, PS_TYPE_F64);
     168        out->n = numElements;
     169        outVec = out->data.F64;
     170
     171        for (unsigned int i = 0; i < numElements; i++) {
     172            outVec[i] = creal(inVec[i]);
     173        }
     174    } else {
     175        char* typeStr;
     176        PS_TYPE_NAME(typeStr,type);
     177        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorReal",
     178                   PS_ERR_BAD_PARAMETER_TYPE, true,
     179                   PS_ERRORTEXT_psVectorFFT_TYPE_UNSUPPORTED,
     180                   typeStr);
    146181        psFree(out);
    147182        return NULL;
    … …  
    186221            outVec[i] = cimagf(inVec[i]);
    187222        }
    188     } else {
    189         psError(__func__, "Can not extract imaginary component from given vector type (%d).", type);
     223    } else if (type == PS_TYPE_C64) {
     224        psF64* outVec;
     225        psC64* inVec = in->data.C64;
     226
     227        out = psVectorRecycle(out, numElements, PS_TYPE_F64);
     228        out->n = numElements;
     229        outVec = out->data.F64;
     230
     231        for (unsigned int i = 0; i < numElements; i++) {
     232            outVec[i] = cimag(inVec[i]);
     233        }
     234    } else {
     235        char* typeStr;
     236        PS_TYPE_NAME(typeStr,type);
     237        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorImaginary",
     238                   PS_ERR_BAD_PARAMETER_TYPE, true,
     239                   PS_ERRORTEXT_psVectorFFT_TYPE_UNSUPPORTED,
     240                   typeStr);
    190241        psFree(out);
    191242        return NULL;
    … …  
    213264
    214265    if (imag->type.type != type) {
    215         psError(__func__, "The inputs to psVectorComplex must be the same type.");
    216         psFree(out);
    217         return NULL;
    218     }
    219 
    220     if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
    221         psError(__func__, "The inputs to psVectorComplex can not be complex.");
     266        char* typeStrReal;
     267        char* typeStrImag;
     268        PS_TYPE_NAME(typeStrReal,type);
     269        PS_TYPE_NAME(typeStrImag,imag->type.type);
     270        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorComplex",
     271                   PS_ERR_BAD_PARAMETER_TYPE, true,
     272                   PS_ERRORTEXT_psVectorFFT_REAL_IMAG_TYPE_MISMATCH,
     273                   typeStrReal,typeStrImag);
    222274        psFree(out);
    223275        return NULL;
    … …  
    236288            outVec[i] = realVec[i] + I * imagVec[i];
    237289        }
    238     } else {
    239         psError(__func__, "Can not merge real and imaginary portions for given vector type (%d).", type);
     290    } else if (type == PS_TYPE_F64) {
     291        psC64* outVec;
     292        psF64* realVec = real->data.F64;
     293        psF64* imagVec = imag->data.F64;
     294
     295        out = psVectorRecycle(out, numElements, PS_TYPE_C64);
     296        out->n = numElements;
     297        outVec = out->data.C64;
     298
     299        for (unsigned int i = 0; i < numElements; i++) {
     300            outVec[i] = realVec[i] + I * imagVec[i];
     301        }
     302    } else {
     303        char* typeStr;
     304        PS_TYPE_NAME(typeStr,type);
     305        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorComplex",
     306                   PS_ERR_BAD_PARAMETER_TYPE, true,
     307                   PS_ERRORTEXT_psVectorFFT_NONREAL_NOTSUPPORTED,
     308                   typeStr);
    240309        psFree(out);
    241310        return NULL;
    … …  
    281350            outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]);
    282351        }
    283     } else {
    284         psError(__func__, "Can not compute complex conjugate for given vector type (%d).", type);
     352    } else if (type == PS_TYPE_C64) {
     353        psC64* outVec;
     354        psC64* inVec = in->data.C64;
     355
     356        out = psVectorRecycle(out, numElements, PS_TYPE_C64);
     357        out->n = numElements;
     358        outVec = out->data.C64;
     359
     360        for (unsigned int i = 0; i < numElements; i++) {
     361            outVec[i] = creal(inVec[i]) - I * cimag(inVec[i]);
     362        }
     363    } else {
     364        char* typeStr;
     365        PS_TYPE_NAME(typeStr,type);
     366        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorConjugate",
     367                   PS_ERR_BAD_PARAMETER_TYPE, true,
     368                   PS_ERRORTEXT_psVectorFFT_NONCOMPLEX_NOTSUPPORTED,
     369                   typeStr);
    285370        psFree(out);
    286371        return NULL;
    … …  
    309394    outNumElements = inHalfNumElements + 1;
    310395
    311     /* if not a complex number, this is not implemented */
    312     if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
    313         psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
    314         psFree(out);
    315         return NULL;
    316     }
    317 
    318396    if (type == PS_TYPE_C32) {
    319397        psF32* outVec;
    … …  
    340418        inAbs1 = cabsf(inVec[inHalfNumElements]);
    341419        outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
    342     } else {
    343         psError(__func__, "Can not power spectrum for given vector type (%d).", type);
    344         psFree(out);
    345         return NULL;
    346     }
    347 
    348     return out;
    349 
    350 }
     420    } else if (type == PS_TYPE_C64) {
     421        psF64* outVec;
     422        psC64* inVec = in->data.C64;
     423        psF64 inAbs1;
     424        psF64 inAbs2;
     425
     426        out = psVectorRecycle(out, outNumElements, PS_TYPE_F64);
     427        out->n = outNumElements;
     428        outVec = out->data.F64;
     429
     430        // from ADD: P_0 = |C_0|^2/N^2
     431        inAbs1 = cabs(inVec[0]);
     432        outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
     433
     434        // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
     435        for (unsigned int i = 1; i < inHalfNumElements; i++) {
     436            inAbs1 = cabs(inVec[i]);
     437            inAbs2 = cabs(inVec[inNumElements - i]);
     438            outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
     439        }
     440
     441        // from ADD: P_N/2 = |C_N/2|^2/N^2
     442        inAbs1 = cabs(inVec[inHalfNumElements]);
     443        outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
     444    } else {
     445        char* typeStr;
     446        PS_TYPE_NAME(typeStr,type);
     447        psErrorMsg(PS_ERRORNAME_DOMAIN "psVectorPowerSpectrum",
     448                   PS_ERR_BAD_PARAMETER_TYPE, true,
     449                   PS_ERRORTEXT_psVectorFFT_NONCOMPLEX_NOTSUPPORTED,
     450                   typeStr);
     451        psFree(out);
     452        return NULL;
     453    }
     454
     455    return out;
     456
     457}
Note: See TracChangeset for help on using the changeset viewer.