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
00190 fftwnd_one(myplan, (fftw_complex*)x, NULL);
00191
00192
00193
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
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;
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;
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;
00240 }
00241
00242 vuVector* computeGradient(t_data* volume)
00243 {
00244
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
00250
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
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
00276
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
00299
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 }