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

FVR/FFT.cpp

Go to the documentation of this file.
00001 #include "Transform.h"
00002 #include "Image_io.h"
00003 #include <stdio.h>
00004 #include <fftw.h>
00005 #include <iostream.h>
00006 
00007 namespace FVR_NS
00008 {
00009 
00010 #define WISDOM "fvr_fftw.wis"
00011 
00012 static fftwnd_plan g_plan_2D;
00013 static bool g_plan_exists = false;
00014 static dword g_x3d;
00015 static dword g_y3d;
00016 static dword g_z3d;
00017 
00018 void initTransforms(void)
00019 {
00020         FILE* wisdom_file = fopen(WISDOM, "a+");
00021         if(!wisdom_file) {
00022             cout<<"Could not load wisdom."<<endl;
00023             return;
00024         }
00025         rewind(wisdom_file);
00026         fftw_import_wisdom_from_file(wisdom_file);
00027         fclose(wisdom_file);
00028 }
00029 
00030 void destroyTransforms(void)
00031 {
00032         FILE* wisdom_file = fopen(WISDOM, "w");
00033         if(!wisdom_file) {
00034             cout << "Could not save wisdom."<<endl;
00035             return;
00036         }
00037         fftw_export_wisdom_to_file(wisdom_file);
00038         fclose(wisdom_file);
00039 }
00040 
00041 void initTransform2D(dword XSize, dword YSize)
00042 {
00043   g_plan_2D = fftw2d_create_plan(YSize, XSize, FFTW_BACKWARD, FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00044   g_plan_exists=true;
00045 }
00046 
00047 void transform2D(t_data* x)
00048 {
00049   if(g_plan_exists)
00050     fftwnd_one(g_plan_2D, (fftw_complex*) x, 0);
00051   else cout << "no plan" << endl;
00052 }
00053 
00054 void destroyTransform2D(void)
00055 {
00056     if(g_plan_exists) {
00057         fftwnd_destroy_plan(g_plan_2D);
00058         g_plan_exists = false;
00059     }
00060 }
00061 
00062 void initTransform3D(dword XSize, dword YSize, dword ZSize)
00063 {
00064         g_x3d = XSize;
00065         g_y3d = YSize;
00066         g_z3d = ZSize;
00067 }
00068 
00069 void transform3D(t_data* x)
00070 {
00071         fftwnd_plan p;
00072         p = fftw3d_create_plan(g_z3d, g_y3d, g_x3d, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00073         
00074         fftwnd_one(p, (fftw_complex*) x, NULL);
00075         fftwnd_destroy_plan(p);
00076         dword size_i = g_z3d * g_y3d * g_x3d * 2;
00077         t_data size_f = (t_data) g_z3d * g_y3d * g_x3d;
00078         for (dword i = 0; i < size_i; ++i)
00079                 x[i] /= size_f;
00080 }
00081 
00082 void destroyTransform3D(void)
00083 {
00084         g_x3d = 0;
00085         g_y3d = 0;
00086         g_z3d = 0;
00087 }
00088 
00089 void shift3D(t_data* x, dword XSize, dword YSize, dword ZSize)
00090 {
00091   dword xstep;
00092   dword ystep;
00093   dword i;
00094   dword j;
00095   dword k;
00096   dword index1;
00097   int start1 = 0;
00098   int val1   = -1;
00099   dword index2;
00100   int start2;
00101   int val2;
00102   
00103   xstep = XSize * 2;
00104   ystep = YSize * XSize * 2;
00105   
00106   index1 = 0;
00107   for(k = 0; k < ZSize; ++k)
00108     {
00109       index2 = index1;
00110       start2 = start1;
00111       val2   = val1;
00112       for(j = 0; j < YSize; ++j)
00113         {
00114           for(i = start2; i < xstep; i += 4)
00115             {
00116               x[index2+i]   *= -1.0f;
00117               x[index2+i+1] *= -1.0f;
00118             }
00119           val2 *= -1;
00120           start2 += val2*2;
00121           index2 += xstep;
00122         }
00123       val1 *= -1;
00124       start1 += val1*2;
00125       index1 += ystep;
00126     }
00127 }
00128 
00129 void shift2D(t_data* x, dword XSize, dword YSize)
00130 {
00131   dword xstep;
00132   dword i;
00133   dword j;
00134   dword index;
00135   int start = 0;
00136   int val = -1;
00137   
00138   xstep = XSize * 2;
00139   
00140   index = 0;
00141   for (j = 0; j < YSize; ++j) {
00142     for (i = start; i < xstep; i += 4) {
00143       x[index+i]   *= -1.0f;
00144       x[index+i+1] *= -1.0f;
00145     }
00146     val *= -1;
00147     start += val * 2;
00148     index += xstep;
00149   }
00150 }
00151 
00152 
00153 
00154 void shift_copy2D(t_data* dest, t_data* src, dword XSize, dword YSize)
00155 {
00156         dword i, j;
00157         dword xlast = XSize / 2;
00158         dword ylast = YSize / 2;
00159 
00160         t_data* src_ptr = src;
00161         t_data* dest_ptr = dest;
00162 
00163         for (j = 0; j < ylast; ++j) {
00164                 for (i = 0; i < xlast; ++i) {
00165                         *(dest_ptr++) = *(src_ptr++) * -1.0f;
00166                         *(dest_ptr++) = *(src_ptr++) * -1.0f;
00167                         *(dest_ptr++) = *(src_ptr++);
00168                         *(dest_ptr++) = *(src_ptr++);
00169                 }
00170                 for (i = 0; i < xlast; ++i) {
00171                         *(dest_ptr++) = *(src_ptr++);
00172                         *(dest_ptr++) = *(src_ptr++);
00173                         *(dest_ptr++) = *(src_ptr++) * -1.0f;
00174                         *(dest_ptr++) = *(src_ptr++) * -1.0f;
00175                 }
00176         }
00177 }
00178 
00179 inline int vcoord(int x, int y, int z)
00180 {
00181   return ((z * g_y3d + y) * g_x3d + x) * 2;
00182 }
00183 
00184 void doStuff(t_data* x, int xMax, int yMax, int zMax)
00185 {
00186   fftwnd_plan myplan;
00187   myplan = fftw3d_create_plan(xMax, yMax, zMax, FFTW_BACKWARD,
00188                          FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00189   //shift3D(x, xMax, yMax, zMax);
00190   fftwnd_one(myplan, (fftw_complex*)x, NULL);
00191   //  dword vsize = xMax * yMax * zMax;
00192   //  for(dword i = 0; i < vsize * 2; i++)
00193   //  x[i] /= (t_data)vsize;
00194   shift3D(x, xMax, yMax, zMax);
00195   fftwnd_destroy_plan(myplan);
00196 }
00197 
00198 t_data eval_f_c(t_data x)
00199 {
00200   if(!((x >=0 && x <= 1)))
00201     cout << "Invalid x value: " << x << endl << flush;
00202 
00203   return x;
00204   
00205   if(x <= .5)
00206     return x;
00207   else
00208     return 1.0 - x;
00209   
00210   //  return sin(x * M_PI - M_PI/2.0);
00211 }
00212 
00213 t_data eval_f_r(t_data x)
00214 {
00215   return .0f;
00216 }
00217 
00218 t_data rfilter(t_data r, t_data c, t_data x)
00219 {
00220   x = x/(t_data)g_x3d;// - M_PI;
00221   t_data cf = eval_f_c(x);
00222   t_data rf = eval_f_r(x);
00223   return r * rf - c * cf;
00224 }
00225 
00226 t_data cfilter(t_data r, t_data c, t_data x)
00227 {
00228   x = x/(t_data)g_x3d;// - M_PI;
00229   t_data cf = eval_f_c(x);
00230   t_data rf = eval_f_r(x);
00231   return rf * c + cf * r;
00232 }
00237 t_data roundThings(t_data value)
00238 {
00239   return value;//fabs(value)<0.000001?0.0f:value;
00240 }
00241 
00242 vuVector* computeGradient(t_data* volume)
00243 {
00244   //  fftwnd_plan p = fftw3d_create_plan(64, 64, 64, FFTW_FORWARD, FFTW_USE_WISDOM);
00245   t_data *tv = new t_data[g_x3d * g_y3d * g_z3d * 2];
00246 
00247   vuVector *v = new vuVector[g_x3d * g_y3d * g_z3d];
00248   cout << "Computing X Component of frequency moment..." << flush;
00249   //create Wx; X component of frequency moment
00250   //memcpy(&tv[0], &volume[0], g_x3d * g_y3d * g_z3d * 2 * sizeof(t_data));
00251 
00252   for (dword i = 0; i < g_x3d; i++)
00253     for(dword j = 0; j < g_y3d; j++)
00254       for(dword k = 0; k < g_z3d; k++)
00255         {
00256           t_data x = (t_data)i;
00257           t_data r = volume[vcoord(i,j,k)];
00258           t_data c = volume[vcoord(i,j,k) + 1];
00259           tv[vcoord(i,j,k)] = rfilter(r, c, x);
00260           tv[vcoord(i,j,k) + 1] = cfilter(r, c, x);
00261         }
00262 
00263   doStuff(tv, g_x3d, g_y3d, g_z3d);
00264 
00265   //cout << "Done with FFT" << flush;
00266   for (dword i = 0; i < g_x3d; i++)
00267     for(dword j = 0; j < g_y3d; j++)
00268       for(dword k = 0; k < g_z3d; k++)
00269         v[vcoord(i,j,k)/2][0] = roundThings(tv[vcoord(i,j,k)]);
00270 
00271   cout << "Done.\n" << flush;
00272 
00273 
00274   cout << "Computing Y Component of frequency moment..." << flush;
00275   //create Wy; Y component of frequency moment
00276   //memcpy(&tv[0], &volume[0], g_x3d * g_y3d * g_z3d * 2 * sizeof(t_data));
00277   for (dword i = 0; i < g_x3d; i++)
00278     for(dword j = 0; j < g_y3d; j++)
00279       for(dword k = 0; k < g_z3d; k++)
00280         {
00281           t_data x = (t_data)j;
00282           t_data r = volume[vcoord(i,j,k)];
00283           t_data c = volume[vcoord(i,j,k) + 1];
00284           tv[vcoord(i,j,k)] = rfilter(r, c, x);
00285           tv[vcoord(i,j,k) + 1] = cfilter(r, c, x);
00286         }
00287   doStuff(tv, g_x3d, g_y3d, g_z3d);
00288 
00289 
00290   for (dword i = 0; i < g_x3d; i++)
00291     for(dword j = 0; j < g_y3d; j++)
00292       for(dword k = 0; k < g_z3d; k++)
00293         v[vcoord(i,j,k)/2][1] = roundThings(tv[vcoord(i,j,k)]);
00294 
00295   cout << "Done.\n" << flush;
00296 
00297   cout << "Computing Z Component of frequency moment..." << flush;
00298   //create Wz; Z component of frequency moment
00299   //memcpy(&tv[0], &volume[0], g_x3d * g_y3d * g_z3d * 2 * sizeof(t_data));
00300   for (dword i = 0; i < g_x3d; i++)
00301     for(dword j = 0; j < g_y3d; j++)
00302       for(dword k = 0; k < g_z3d; k++)
00303         {
00304           t_data x = (t_data)k;
00305           t_data r = volume[vcoord(i,j,k)];
00306           t_data c = volume[vcoord(i,j,k) + 1];
00307           tv[vcoord(i,j,k)] = rfilter(r, c, x);
00308           tv[vcoord(i,j,k) + 1] = cfilter(r, c, x);
00309         }
00310 
00311   doStuff(tv, g_x3d, g_y3d, g_z3d);
00312 
00313   for (dword i = 0; i < g_x3d; i++)
00314     for(dword j = 0; j < g_y3d; j++)
00315       for(dword k = 0; k < g_z3d; k++)
00316         v[vcoord(i,j,k)/2][2] = roundThings(tv[vcoord(i,j,k)]);
00317 
00318   cout << "Done.\n" << flush;
00319 
00320 
00321   delete tv;
00322   return v;
00323 }
00324 
00325 } // NS_END

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