00001 #include "SHarmonics.h"
00002 #include "TorstensFilters.h"
00003 #include "vuVector.h"
00004 #include "Transform.h"
00005 #include "Image_io.h"
00006 #include <math.h>
00007 #include <string.h>
00008 #include <stdio.h>
00009
00010 using namespace FVR_NS;
00011
00012 #define ANGLE_TOL 20.0/180.0 * M_PI
00013
00014
00015 t_data dc[]={3.14198,
00016 2.09440,2.09440,2.09440,
00017 0.78520,0.78520,0.78520,0.78520,0.78520};
00018 bool downsample = false;
00019 bool inRange(float v, float value, float error)
00020 {
00021 return (v < value + error) && (v > value - error);
00022 }
00023 int indexOf(int x, int xMax, int y, int yMax, int z, int zMax)
00024 {
00025 return 2 * (z * yMax * xMax + y * xMax + x);
00026 }
00027
00028 void realTrans(t_data* volume, int& xMax, int& yMax, int& zMax, t_data* Yarray[10])
00029 {
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 int subSample = 2;
00053 if ((xMax > 128 || yMax > 128 || zMax > 128) && false )
00054 {
00055 printf("We subsample!\n");
00056 xMax /= subSample;
00057 yMax /= subSample;
00058 zMax /= subSample;
00059 downsample = true;
00060 }
00061
00062
00063 int L = zMax; int M = yMax; int N = xMax;
00064 t_data x,y,z,delf,phshift;
00065 int volumeSize = L * M * N * 2;
00066 Yarray[0] = new t_data[volumeSize];
00067
00068 if(downsample)
00069 for(int k = 0; k < L; k++)
00070 for(int j = 0; j < M; j++)
00071 for(int i = 0; i < N; i++)
00072 {
00073 int destIdx = indexOf(i,N,j,M,k,L);
00074
00075
00076
00077
00078
00079 t_data accum = 0;
00080 for(int avgx = 0; avgx < subSample; avgx++)
00081 for(int avgy = 0; avgy < subSample; avgy++)
00082 for(int avgz = 0; avgz < subSample; avgz++)
00083 accum += volume[indexOf(subSample * i + avgx, subSample * N,
00084 subSample * j + avgy, subSample * M,
00085 subSample * k + avgz, subSample * L)];
00086 Yarray[0][destIdx] = accum/(t_data)(subSample * subSample * subSample);
00087 accum = 0;
00088 for(int avgx = 0; avgx < subSample; avgx++)
00089 for(int avgy = 0; avgy < subSample; avgy++)
00090 for(int avgz = 0; avgz < subSample; avgz++)
00091 accum += volume[indexOf(subSample * i + avgx + 1, subSample * N,
00092 subSample * j + avgy + 1, subSample * M,
00093 subSample * k + avgz + 1, subSample * L)];
00094
00095 Yarray[0][destIdx+1] = accum/(t_data)(subSample * subSample * subSample);
00096 }
00097
00098 else
00099 memcpy(Yarray[0], volume, volumeSize * sizeof(t_data));
00100
00106 initTransform3D(N, M, L);
00107 shift3D(volume, N, M, L);
00108 transform3D(volume);
00109
00110 destroyTransform3D();
00111
00112 initTransform3D(N, M, L);
00113 vuVector* grads = computeGradient(volume);
00114 destroyTransform3D();
00115
00116 memcpy(volume, Yarray[0], volumeSize * sizeof(t_data));
00117
00119
00120 Yarray[1] = new t_data[volumeSize];
00121 memset(Yarray[1], 0, volumeSize * sizeof(t_data));
00122 Yarray[2] = new t_data[volumeSize];
00123 memset(Yarray[2], 0, volumeSize * sizeof(t_data));
00124 Yarray[3] = new t_data[volumeSize];
00125 memset(Yarray[3], 0, volumeSize * sizeof(t_data));
00126 Yarray[4] = new t_data[volumeSize];
00127 memset(Yarray[4], 0, volumeSize * sizeof(t_data));
00128 Yarray[5] = new t_data[volumeSize];
00129 memset(Yarray[5], 0, volumeSize * sizeof(t_data));
00130 Yarray[6] = new t_data[volumeSize];
00131 memset(Yarray[6], 0, volumeSize * sizeof(t_data));
00132 Yarray[7] = new t_data[volumeSize];
00133 memset(Yarray[7], 0, volumeSize * sizeof(t_data));
00134 Yarray[8] = new t_data[volumeSize];
00135 memset(Yarray[8], 0, volumeSize * sizeof(t_data));
00136 Yarray[9] = new t_data[volumeSize];
00137 memset(Yarray[9], 0, volumeSize * sizeof(t_data));
00138
00139 float minAll = 100.0;
00140 float maxAll = 100.0;
00141 float sumLen = 0.0;
00142 float minLen = 100.0;
00143 float maxLen = -100.0;
00144 dword numOfOvers = 0;
00145 float largestTheta = 0.0;
00146 for(int l=1;l<L-1;l++)
00147 for(int m=1;m<M-1;m++)
00148 for(int n=1;n<N-1;n++) {
00149 Yarray[1][indexOf(n,N,m,M,l,L)+1]=0;
00150 Yarray[2][indexOf(n,N,m,M,l,L)+1]=0;
00151 Yarray[3][indexOf(n,N,m,M,l,L)+1]=0;
00152 Yarray[4][indexOf(n,N,m,M,l,L)+1]=0;
00153 Yarray[5][indexOf(n,N,m,M,l,L)+1]=0;
00154 Yarray[6][indexOf(n,N,m,M,l,L)+1]=0;
00155 Yarray[7][indexOf(n,N,m,M,l,L)+1]=0;
00156 Yarray[8][indexOf(n,N,m,M,l,L)+1]=0;
00157 Yarray[9][indexOf(n,N,m,M,l,L)+1]=0;
00158
00159
00160
00161
00162
00163
00164
00165
00166 x = (Yarray[0][indexOf(n+1,N,m,M,l,L)] -
00167 Yarray[0][indexOf(n-1,N,m,M,l,L)]);
00168 y = (Yarray[0][indexOf(n,N,m+1,M,l,L)] -
00169 Yarray[0][indexOf(n,N,m-1,M,l,L)]);
00170 z = (Yarray[0][indexOf(n,N,m,M,l+1,L)] -
00171 Yarray[0][indexOf(n,N,m,M,l-1,L)]);
00172
00173 vuVector normal = grads[indexOf(n,N,m,M,l,L)/2];
00174
00175
00176 vuVector cd(x,y,z);
00177
00178
00179
00180 float nl = normal.norm();
00181 float cdl = cd.norm();
00182 float theta = acos(normal.dot(cd) / (nl * cdl));
00183 if(minAll > nl)
00184 minAll = nl;
00185 if(maxAll < nl)
00186 maxAll = nl;
00187
00188 if(theta > ANGLE_TOL || cdl == 0)
00189 {
00190 numOfOvers++;
00191 float len = normal.norm();
00192 if(maxLen < len)
00193 maxLen = len;
00194 if(minLen > len)
00195 minLen = len;
00196 sumLen += len;
00197 normal = vuVector(0,0,0);
00198 }
00199
00200
00201 if(largestTheta < theta)
00202 largestTheta = theta;
00203
00204 if(false && theta > ANGLE_TOL)
00205 {
00206 cout << "(" << n << ", " << m << ", " << l << ")";
00207 cout << "Grad "; normal.print();
00208 cout << " CD "; cd.print();
00209 cout << endl << flush;
00210 }
00211
00212 normal = cd;
00213
00214 x = normal[0];
00215 y = normal[1];
00216 z = normal[2];
00217
00218 if(x==0 && y==0 && z==0) {
00219 Yarray[1][indexOf(n,N,m,M,l,L)]=0;
00220 Yarray[2][indexOf(n,N,m,M,l,L)]=0;
00221 Yarray[3][indexOf(n,N,m,M,l,L)]=0;
00222 Yarray[4][indexOf(n,N,m,M,l,L)]=0;
00223 Yarray[5][indexOf(n,N,m,M,l,L)]=0;
00224 Yarray[6][indexOf(n,N,m,M,l,L)]=0;
00225 Yarray[7][indexOf(n,N,m,M,l,L)]=0;
00226 Yarray[8][indexOf(n,N,m,M,l,L)]=0;
00227 Yarray[9][indexOf(n,N,m,M,l,L)]=0;
00228 }
00229 else {
00230 int tst=l+m+n;
00231
00232 phshift= tst % 2 == 0 ? 1:1;
00233 phshift*=Yarray[0][indexOf(n,N,m,M,l,L)];
00234
00235 delf = x*x;
00236 delf += y*y;
00237 delf += z*z;
00238 delf = sqrt(delf);
00239 x/=delf;
00240 y/=delf;
00241 z/=delf;
00242
00243
00244
00245
00246 Yarray[1][indexOf(n,N,m,M,l,L)]=.282096*phshift;
00247 Yarray[2][indexOf(n,N,m,M,l,L)]=.488603*y*phshift;
00248 Yarray[3][indexOf(n,N,m,M,l,L)]=.488603*z*phshift;
00249 Yarray[4][indexOf(n,N,m,M,l,L)]=.488603*x*phshift;
00250 Yarray[5][indexOf(n,N,m,M,l,L)]=1.092548*x*y*phshift;
00251 Yarray[6][indexOf(n,N,m,M,l,L)]=1.092548*z*y*phshift;
00252 Yarray[7][indexOf(n,N,m,M,l,L)]=.315392*(3*z*z-1)*phshift;
00253 Yarray[8][indexOf(n,N,m,M,l,L)]=1.092548*z*x*phshift;
00254 Yarray[9][indexOf(n,N,m,M,l,L)]=.546274*(x*x-y*y)*phshift;
00255 }
00256 }
00257
00258 cout << "Largest Difference between CD and Exact was " << largestTheta * 180.0 /M_PI << endl << flush;
00259
00260 cout << numOfOvers << " vectors thrown away. Avg: " << sumLen/numOfOvers
00261 << " min: " << minLen << " max: " << maxLen;
00262 cout << "Min Normal: " << minAll << "max Normal: " << maxAll << endl << flush;
00263 return;
00264 }
00265
00266 void realLight(t_data sphy[], t_data lv[3])
00267 {
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279 t_data sqrt2=sqrt(2.0);
00280
00281 t_data x=lv[0], y=lv[1], z=lv[2];
00282
00283 sphy[0]=sqrt(1.0/4.0/M_PI);
00284 sphy[1]=y*sqrt(3.0/8.0/M_PI)*sqrt2;
00285 sphy[2]=sqrt(3.0/4.0/M_PI)*z;
00286 sphy[3]=x*sqrt(3.0/8.0/M_PI)*sqrt2;
00287 sphy[4]=2*x*y*sqrt(15.0/2.0/M_PI)/4.0*sqrt2;
00288 sphy[5]=z*y*sqrt(15.0/8.0/M_PI)*sqrt2;
00289 sphy[6]=sqrt(5.0/4.0/M_PI)*(pow(1.5*z,2)-.5);
00290 sphy[7]=z*x*sqrt(15.0/8.0/M_PI)*sqrt2;
00291 sphy[8]=(x*x-y*y)*sqrt(15.0/2.0/M_PI)/4.0*sqrt2;
00292 return;
00293 }
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371