Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

vuFourierVolume_IO.cpp

Go to the documentation of this file.
00001 #include "vuFourierVolume_IO.h"
00002 #include "Volume/Fourier/Unimodal/3d/3d.h"
00003 
00004 #define MAX_FILTER_SIZE 6
00005 
00006 template <int S> bool vuFourierVolume_IO<S>::
00007 preprocessSpatialInput(byte *data, dword XSize, dword YSize, dword ZSize,
00008                        float s, float mult_pad, dword add_pad,
00009                        bool doWrapAndInit)
00010 {
00011   dword dim = computeDimensions(XSize, YSize, ZSize, mult_pad, add_pad);
00012 
00013   m_XSize  = m_YSize = m_ZSize = dim;
00014   m_Volume = new float[m_XSize * m_YSize * m_ZSize * 2 * S];
00015 
00016   // transporm bytes into complex numbers + do some padding
00017   readSpatial(data, m_Volume, m_XSize, m_YSize, m_ZSize, XSize, YSize, ZSize);
00018 
00019   // apply fft on spatial domain
00020   preprocess();
00021 
00022   if (doWrapAndInit) {
00023     setWrap(MAX_FILTER_SIZE/2);
00024 
00025     XSize   = m_XSize;
00026     YSize   = m_YSize;
00027     ZSize   = m_ZSize;
00028 
00029     m_XSize = XSize + m_Wrap * 2;
00030     m_YSize = YSize + m_Wrap * 2;
00031     m_ZSize = ZSize + m_Wrap * 2;
00032 
00033     float *dest = new float[m_XSize * m_YSize * m_ZSize * 2 * S];
00034   
00035     // make m_Volume larger (according to m_Wrap)
00036     padFourier(m_Volume, dest, m_XSize, m_YSize, m_ZSize, XSize, YSize, ZSize);
00037     CHECKNDELETE(m_Volume);
00038     m_Volume = dest;
00039 
00040     wrapAndInitialize(s);
00041   }
00042   else
00043     m_Wrap = 0;
00044 
00045   return true;
00046 }
00047 
00048 template <int S>
00049 void vuFourierVolume_IO<S>::preprocess()
00050 {
00051   cout << "shifting volume. " << flush;
00052   shift3D(m_Volume, m_XSize, m_YSize, m_ZSize);
00053 
00054   cout << "transforming volume. " << flush;
00055   transform3D(m_Volume, m_XSize, m_YSize, m_ZSize);
00056 
00057   cout << "shifting volume. " << flush;
00058   shift3D(m_Volume, m_XSize, m_YSize, m_ZSize);
00059   cout << "done\n" << flush;
00060 }
00061 
00062 template <int S> void vuFourierVolume_IO<S>::
00063 transform3D(float* vol, dword xx, dword yy, dword zz)
00064 {
00065   fftwnd_plan plan;
00066 
00067   plan = fftw3d_create_plan(xx, yy, zz, FFTW_FORWARD,
00068                             FFTW_ESTIMATE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00069 
00070   fftwnd(plan, S, (fftw_complex*)vol, S, 1, NULL, 0, 0);
00071 
00072   fftwnd_destroy_plan(plan);
00073   dword count  =  xx * yy * zz * 2 * S;
00074   float size_f = (float) xx* yy * zz;
00075   float *ptr   = vol;
00076   for (dword i = 0; i < count; ++i) *(ptr++) /= size_f;
00077 }
00078 
00079 template <int S> void vuFourierVolume_IO<S>::
00080 inverseTransform3D(float* vol, dword xx, dword yy, dword zz)
00081 {
00082   fftwnd_plan plan;
00083 
00084   plan = fftw3d_create_plan(xx, yy, zz, FFTW_BACKWARD,
00085                             FFTW_ESTIMATE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00086 
00087   fftwnd(plan, S, (fftw_complex*)vol, S, 1, NULL, 0, 0);
00088 
00089   fftwnd_destroy_plan(plan);
00090 }
00091 
00092 
00093 template <int S> void vuFourierVolume_IO<S>::
00094 shift3D(float* data,dword XSize,dword YSize,dword ZSize)
00095 {
00096   dword xstep  = XSize * 2 * S;
00097   dword ystep  = YSize * XSize * 2 * S;
00098   dword index1 = 0;
00099   int   start1 = 0;
00100   int   val1   = -1;
00101 
00102   dword index2;
00103   int   val2, start2;
00104 
00105   for(dword k=0; k<ZSize; k++) {
00106     index2 = index1;
00107     start2 = start1;
00108     val2   = val1;
00109     for(dword j=0; j<YSize; j++) {
00110       for(dword i=start2; i<xstep; i+=4*S) {
00111         float *ptr = &data[index2+i];
00112         for (int l=0; l<S; l++) {
00113           *(ptr++) *= -1.0f;
00114           *(ptr++) *= -1.0f;
00115         }
00116       }
00117       val2 *= -1;
00118       start2 += val2 * 2 * S;
00119       index2 += xstep;
00120     }
00121     val1 *= -1;
00122     start1 += val1*2*S;
00123     index1 += ystep;
00124   }
00125 }
00126 
00127 template <int S> bool vuFourierVolume_IO<S>::
00128 scaleAndWriteToFourierFile(const char *fileName, float scale)
00129 {
00130   if (scale !=1.0 && scale > 0 || true) {
00131     dword count = m_XSize * m_YSize * m_ZSize * 2 * S;
00132     float *ptr  = m_Volume;
00133     for (dword i=0; i<count; i++) *(ptr++) *= scale;
00134   }
00135 
00136   FILE *file = fopen(fileName, "wb");
00137   if (file == NULL) return false;
00138 
00139   dword XSize = m_XSize - 2 * m_Wrap;
00140   dword YSize = m_YSize - 2 * m_Wrap;
00141   dword ZSize = m_ZSize - 2 * m_Wrap;
00142 
00143   bool status = vu1712<S>::writeHeader(file, XSize, YSize, ZSize);
00144   if (status) {
00145     status = vu1712<S>::writeData(file, m_Volume, 
00146                                   m_XSize, m_YSize, m_ZSize,
00147                                   XSize, YSize, ZSize);
00148   }
00149   fclose(file);
00150   return status;
00151 }
00152 
00153 /*
00154    in  -> input data  (spatial domain volume)
00155    out -> output data (fourier domain volume)
00156 */
00157 template <int S>
00158 void vuFourierVolume_IO<S>::readSpatial(byte *in, float* out,
00159                                    dword XX, dword YY, dword ZZ,
00160                                    dword XXsmall, dword YYsmall, dword ZZsmall)
00161 {
00162   byte  *i_ptr = in;
00163   float *v_ptr = out;
00164 
00165   dword lowX = (XX - XXsmall) / 2;
00166   dword lowY = (YY - YYsmall) / 2;
00167   dword lowZ = (ZZ - ZZsmall) / 2;
00168 
00169   dword hiX = XX - lowX - XXsmall;
00170   dword hiY = YY - lowY - YYsmall;
00171   dword hiZ = ZZ - lowZ - ZZsmall;
00172 
00173   v_ptr = pad(v_ptr, lowZ * XX * YY * 2 * S);
00174   for (dword k = 0; k < ZZsmall; k++) {
00175     v_ptr = pad(v_ptr, lowY * XX * 2 * S);
00176     for (dword j = 0; j < YYsmall; j++) {
00177       v_ptr = pad(v_ptr, lowX * 2 * S);
00178       for (dword i = 0; i < XXsmall; i++) {
00179         for (int l=0; l<S; l++) {
00180           *(v_ptr++) = *(i_ptr++);
00181           *(v_ptr++) = 0.0f;
00182         }
00183       }
00184       v_ptr = pad(v_ptr, hiX * 2 * S);
00185     }
00186     v_ptr = pad(v_ptr, hiY * XX * 2 * S);
00187   }
00188   v_ptr = pad(v_ptr, hiZ * XX * YY * 2 * S);
00189 }
00190 
00191 // complex pad n spaces
00192 template <int S>
00193 float* vuFourierVolume_IO<S>::pad(float* v, dword n)
00194 {
00195   for (dword i = 0; i < n; ++i) *(v++) = 0.0f;
00196   return v;
00197 }
00198 
00199 template <int S>
00200 void vuFourierVolume_IO<S>::padFourier(float *in, float* out,
00201                                    dword XX, dword YY, dword ZZ,
00202                                    dword XXsmall, dword YYsmall, dword ZZsmall)
00203 {
00204   float *i_ptr = in;
00205   float *v_ptr = out;
00206 
00207   dword lowX = (XX - XXsmall) / 2;
00208   dword lowY = (YY - YYsmall) / 2;
00209   dword lowZ = (ZZ - ZZsmall) / 2;
00210 
00211   dword hiX = XX - lowX - XXsmall;
00212   dword hiY = YY - lowY - YYsmall;
00213   dword hiZ = ZZ - lowZ - ZZsmall;
00214 
00215   v_ptr = pad(v_ptr, lowZ * XX * YY * 2 * S);
00216   for (dword k = 0; k < ZZsmall; k++) {
00217     v_ptr = pad(v_ptr, lowY * XX * 2 * S);
00218     for (dword j = 0; j < YYsmall; j++) {
00219       v_ptr = pad(v_ptr, lowX * 2 * S);
00220       for (dword i = 0; i < XXsmall; i++) {
00221         for (int l=0; l<S; l++) {
00222           *(v_ptr++) = *(i_ptr++);
00223           *(v_ptr++) = *(i_ptr++);
00224         }
00225       }
00226       v_ptr = pad(v_ptr, hiX * 2 * S);
00227     }
00228     v_ptr = pad(v_ptr, hiY * XX * 2 * S);
00229   }
00230   v_ptr = pad(v_ptr, hiZ * XX * YY * 2 * S);
00231 }
00232 
00233 template <int S>
00234 bool vuFourierVolume_IO<S>::getSpatialDataFromVUF(byte *&spatVol,
00235                                                   dword &XSize,
00236                                                   dword &YSize,
00237                                                   dword &ZSize,
00238                                                   float scale,
00239                                                   vuString fileName)
00240 {
00241   CHECKNDELETE(spatVol);
00242 
00243   if (S != 1) {
00244     cerr << "vuFourierVolume_IO::getSpatialDataFromVUF(): ";
00245     cerr << "only S=1 is supported" << endl;
00246     return false;
00247   }
00248 
00249   float *freqVol = NULL;
00250   bool  status   = true;
00251   if (fileName.isEmpty()) return false;
00252   FILE *file  = fopen(fileName, "rb");
00253   if (file == NULL) return false;
00254 
00255   status = vu1712<S>::readHeader(file, XSize, YSize, ZSize);
00256 
00257   if (status) {
00258     freqVol = new float[XSize * YSize * ZSize * S * 2];
00259     status = vu1712<S>::readData(file, freqVol,
00260                                  XSize, YSize, ZSize,
00261                                  XSize, YSize, ZSize);
00262   }
00263   fclose(file);
00264 
00265   if (status) {
00266     vuFourierVolume_IO<S>::inverseTransform3D(freqVol, XSize, YSize, ZSize);
00267 
00268     dword count  = XSize * YSize * ZSize;
00269     float *s_ptr = freqVol;
00270     spatVol      = new byte[count];
00271     byte  *d_ptr = spatVol;
00272 
00273     if (scale == 0) scale = 1.0;
00274 
00275     for (dword i=0; i<count; i++) {
00276       register float tmp = *s_ptr * scale;
00277 
00278       if (tmp > 255)
00279         tmp = 255.0f;
00280       else if (tmp < 0)
00281         tmp = 0.0f;
00282 
00283       *d_ptr = (byte)tmp;
00284 
00285       d_ptr += 1;
00286       s_ptr += 2;
00287     }
00288   }
00289 
00290   CHECKNDELETE(freqVol);
00291   return status;
00292 }

Generated on Wed Dec 15 21:20:34 2004 for vuVolume by  doxygen 1.3.9.1