IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3702 for trunk/psLib/src/image


Ignore:
Timestamp:
Apr 14, 2005, 2:12:09 PM (21 years ago)
Author:
desonia
Message:

fixed image FFT and image convolution. The FFTW call had the wrong order
of dimensions, and the scaling was incorrect.

Location:
trunk/psLib/src/image
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r3682 r3702  
    99 *  @author Ross Harman, MHPCC
    1010 *
    11  *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-04-07 20:27:41 $
     11 *  @version $Revision: 1.64 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-04-15 00:12:08 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    197197}
    198198
    199 psImage* psImageCopy(psImage* restrict output,
     199psImage* psImageCopy(psImage* output,
    200200                     const psImage* input,
    201201                     psElemType type)
    … …  
    350350    return output;
    351351}
     352
     353bool p_psImageCopyToRawBuffer(void* buffer,
     354                              const psImage* input,
     355                              psElemType type)
     356{
     357    psElemType inDatatype;
     358    psS32 numRows;
     359    psS32 numCols;
     360
     361    if (input == NULL || input->data.V == NULL) {
     362        psError(PS_ERR_BAD_PARAMETER_NULL, true,
     363                PS_ERRORTEXT_psImage_IMAGE_NULL);
     364        return false;
     365    }
     366
     367    if (input->type.dimen != PS_DIMEN_IMAGE) {
     368        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     369                PS_ERRORTEXT_psImage_NOT_AN_IMAGE);
     370        return false;
     371    }
     372
     373    inDatatype = input->type.type;
     374    numRows = input->numRows;
     375    numCols = input->numCols;
     376
     377    // cover the trival case of copy of the same
     378    // datatype.
     379    if (type == inDatatype) {
     380        int rowSize = PSELEMTYPE_SIZEOF(inDatatype)*numCols;
     381        for (psS32 row=0;row<numRows;row++) {
     382            memcpy(&((psS8*)buffer)[row*rowSize], input->data.V[row], rowSize);
     383        }
     384        return true;
     385    }
     386
     387    #define PSIMAGE_BUFFER_COPY(INTYPE,OUTTYPE) { \
     388        ps##INTYPE *in; \
     389        ps##OUTTYPE *out = buffer; \
     390        for(psS32 row=0;row<numRows;row++) { \
     391            in = input->data.INTYPE[row]; \
     392            for (psS32 col=0;col<numCols;col++) { \
     393                *(out++) = *(in++); \
     394            } \
     395        } \
     396    }
     397
     398    #define PSIMAGE_BUFFER_COPY_CASE(OUT,OUTTYPE) { \
     399        switch (inDatatype) { \
     400        case PS_TYPE_S8: \
     401            PSIMAGE_BUFFER_COPY(S8,OUTTYPE); \
     402            break; \
     403        case PS_TYPE_S16: \
     404            PSIMAGE_BUFFER_COPY(S16,OUTTYPE); \
     405            break; \
     406        case PS_TYPE_S32: \
     407            PSIMAGE_BUFFER_COPY(S32,OUTTYPE); \
     408            break; \
     409        case PS_TYPE_S64: \
     410            PSIMAGE_BUFFER_COPY(S64,OUTTYPE); \
     411            break; \
     412        case PS_TYPE_U8: \
     413            PSIMAGE_BUFFER_COPY(U8,OUTTYPE); \
     414            break; \
     415        case PS_TYPE_U16: \
     416            PSIMAGE_BUFFER_COPY(U16,OUTTYPE); \
     417            break; \
     418        case PS_TYPE_U32: \
     419            PSIMAGE_BUFFER_COPY(U32,OUTTYPE); \
     420            break; \
     421        case PS_TYPE_U64: \
     422            PSIMAGE_BUFFER_COPY(U64,OUTTYPE); \
     423            break; \
     424        case PS_TYPE_F32: \
     425            PSIMAGE_BUFFER_COPY(F32,OUTTYPE); \
     426            break; \
     427        case PS_TYPE_F64: \
     428            PSIMAGE_BUFFER_COPY(F64,OUTTYPE); \
     429            break; \
     430        case PS_TYPE_C32: \
     431            PSIMAGE_BUFFER_COPY(C32,OUTTYPE); \
     432            break; \
     433        case PS_TYPE_C64: \
     434            PSIMAGE_BUFFER_COPY(C64,OUTTYPE); \
     435            break; \
     436        default: \
     437            break; \
     438        } \
     439    }
     440
     441    switch (type) {
     442    case PS_TYPE_S8:
     443        PSIMAGE_BUFFER_COPY_CASE(output, S8);
     444        break;
     445    case PS_TYPE_S16:
     446        PSIMAGE_BUFFER_COPY_CASE(output, S16);
     447        break;
     448    case PS_TYPE_S32:
     449        PSIMAGE_BUFFER_COPY_CASE(output, S32);
     450        break;
     451    case PS_TYPE_S64:
     452        PSIMAGE_BUFFER_COPY_CASE(output, S64);
     453        break;
     454    case PS_TYPE_U8:
     455        PSIMAGE_BUFFER_COPY_CASE(output, U8);
     456        break;
     457    case PS_TYPE_U16:
     458        PSIMAGE_BUFFER_COPY_CASE(output, U16);
     459        break;
     460    case PS_TYPE_U32:
     461        PSIMAGE_BUFFER_COPY_CASE(output, U32);
     462        break;
     463    case PS_TYPE_U64:
     464        PSIMAGE_BUFFER_COPY_CASE(output, U64);
     465        break;
     466    case PS_TYPE_F32:
     467        PSIMAGE_BUFFER_COPY_CASE(output, F32);
     468        break;
     469    case PS_TYPE_F64:
     470        PSIMAGE_BUFFER_COPY_CASE(output, F64);
     471        break;
     472    case PS_TYPE_C32:
     473        PSIMAGE_BUFFER_COPY_CASE(output, C32);
     474        break;
     475    case PS_TYPE_C64:
     476        PSIMAGE_BUFFER_COPY_CASE(output, C64);
     477        break;
     478    default: {
     479            char* typeStr;
     480            PS_TYPE_NAME(typeStr,type);
     481            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
     482                    PS_ERRORTEXT_psImage_IMAGE_TYPE_UNSUPPORTED,
     483                    typeStr);
     484            break;
     485        }
     486    }
     487    return true;
     488}
     489
    352490
    353491psS32 psImageFreeChildren(psImage* image)

  • trunk/psLib/src/image/psImage.h

    r3446 r3702  
    1111 *  @author Ross Harman, MHPCC
    1212 *
    13  *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-03-18 02:35:14 $
     13 *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-04-15 00:12:08 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    167167);
    168168
     169bool p_psImageCopyToRawBuffer(
     170    void* buffer,
     171    const psImage* input,
     172    psElemType type
     173);
     174
    169175/** Frees all children of a psImage.
    170176 *

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

    r3682 r3702  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-04-07 20:27:41 $
     7 *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-04-15 00:12:08 $
    99 *
    1010 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    398398        }
    399399
    400         psImageWriteSection(paddedImage,0,0,0,NULL,0,"paddedImage.fits");
    401 
    402400        // pad the kernel to the same size of paddedImage
    403401        psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL);
    … …  
    421419        }
    422420
    423         psImageWriteSection(paddedKernel,0,0,0,NULL,0,"paddedKernel.fits");
    424 
    425421        psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD);
    426422        if (kernelFourier == NULL) {
    … …  
    440436
    441437        // convolution in fourier domain is just a pixel-wise multiplication
    442         psBinaryOp(inFourier,inFourier,"*",kernelFourier);
     438        for (int row = 0; row < paddedRows; row++) {
     439            psC32* inRow = inFourier->data.C32[row];
     440            psC32* kRow = kernelFourier->data.C32[row];
     441            for (int col = 0; col < paddedCols; col++) {
     442                inRow[col] *= kRow[col];
     443            }
     444        }
    443445
    444446        psImage* complexOut = psImageFFT(NULL, inFourier,
    … …  
    451453        }
    452454
    453         {
    454             psImage* tmp = psImageReal(NULL,complexOut);
    455             psImageWriteSection(tmp,0,0,0,NULL,0,"complexOut.fits");
    456             psFree(tmp);
    457         }
    458 
    459455        // subset out the padded area now.
    460456        psImage* complexOutSansPad = psImageSubset(complexOut,
    … …  
    463459
    464460        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
    465         float factor = 1.0f/numCols/numRows;
     461        float factor = 1.0f/(float)paddedCols/(float)paddedRows;
    466462        for (psS32 row = 0; row < numRows; row++) {
    467463            psF32* outRow = out->data.F32[row];

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

    r3671 r3702  
    55 *  @author Robert DeSonia, MHPCC
    66 *
    7  *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
    8  *  @date $Date: 2005-04-06 01:12:58 $
     7 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
     8 *  @date $Date: 2005-04-15 00:12:08 $
    99 *
    1010 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    7474    // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place.
    7575    psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD;
    76     out = psImageCopy(out, in, PS_TYPE_C32);
    77     plan = fftwf_plan_dft_2d(numCols, numRows,
    78                              out->data.C32[0],
    79                              out->data.C32[0],
     76
     77    fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex));
     78    p_psImageCopyToRawBuffer(outBuffer, in, PS_TYPE_C32);
     79
     80    plan = fftwf_plan_dft_2d(numRows, numCols,
     81                             outBuffer,
     82                             outBuffer,
    8083                             sign,
    8184                             PS_FFTW_PLAN_RIGOR);
    … …  
    101104        // becomes an issue, the use of fftwf_plan_dft_r2c should be considered
    102105        // as well as fftwf_plan_dft_c2r.
    103         psImage* realOut = psImageCopy(NULL,out,PS_TYPE_F32);
    104         psFree(out);
    105         out = realOut;
    106     }
     106        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
     107        int index = 0;
     108        for (int row=0; row < numRows; row++) {
     109            psF32* outRow = out->data.F32[row];
     110            for (int col=0; col < numCols; col++) {
     111                outRow[col] = crealf(outBuffer[index++]); // take just the real part
     112            }
     113        }
     114    } else {
     115        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
     116        int index = 0;
     117        for (int row=0; row < numRows; row++) {
     118            psC32* outRow = out->data.C32[row];
     119            for (int col=0; col < numCols; col++) {
     120                outRow[col] = outBuffer[index++]; // take just the real part
     121            }
     122        }
     123        //        memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex));
     124    }
     125
     126    fftwf_free(outBuffer);
    107127
    108128    return out;
Note: See TracChangeset for help on using the changeset viewer.