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

SHarmonics.cpp

Go to the documentation of this file.
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 // the harmonic transform coefficients for
00014 // diffuse BRDF with max function
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   /* Real version of the spherical harmonic transform
00031      (aka tesseract harmonic transform)
00032      The main differences between this an complex are:
00033      
00034      for Ylm m!= 0, multiply by sqrt(2)
00035      all coefficients are positive
00036 
00037      this implementation follows Ravi Rammamorhi's in
00038      "Efficient Representation of Radiance Maps" from
00039      SIGGRAPH 2001
00040 
00041      x,y,z = unit normal
00042 
00043      Note that since we use real values. If I had used the real
00044      version of the fftw software I could gain a factor of 2 in time
00045      and storage, but this works and was easiest to get started.
00046 
00047      ***
00048      put Yarray[0] the volume itself
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             //int sourceIdx = indexOf(subSample * i, subSample * N,
00076             //subSample * j, subSample * M,
00077             //subSample * k, subSample * L);
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   //shift3D(volume, N, M, L);
00110   destroyTransform3D();
00111 
00112   initTransform3D(N, M, L);
00113   vuVector* grads = computeGradient(volume);
00114   destroyTransform3D();
00115   // spitVec3d(grads, N*M*L, "/tmp/grad.exact");
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   //return;
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           x = (rinp[l+1][m][n]-rinp[l-1][m][n]);
00161           y = (rinp[l][m+1][n]-rinp[l][m-1][n]);
00162           z = (rinp[l][m][n+1]-rinp[l][m][n-1]);
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         //normal = normal.MakeUnit();
00176         vuVector cd(x,y,z);
00177         //cd = cd.MakeUnit();
00178         
00179         //if(inRange(n,N/2.0, 2) && inRange(m,M/2.0, 2) && inRange(l,L/2.0, 2))
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         //if(cd.norm() > 0.005 || normal.norm() > 0.005)
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         //use the central differencing method.
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           // shift 0 freq to center of array with phshift
00232           phshift= tst % 2 == 0  ? 1: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           // the real spherical harmonics for l=0..2
00245           // in terms of surface normal (x,y,z)
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   /* Real version of the spherical harmonic transform This is the same
00270      as all_realY except it does only one evaluation. I use it to
00271      compute the coefficients for a directional light, which is the
00272      harmonic evaluated for the light direction. The calculation is
00273      exactly the same as for all_realY, but the coefficients are
00274      calculated explicitly to show how they are obtained.
00275  
00276      x,y,z = unit normal
00277      sphy = array with light coefficients */
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 #define clip(a) (a<0 ? 0 : a) // remove neg values due to approx in image
00297 void shade(t_data lv[3], t_data*& pic)
00298 {
00299   t_data lc[9];
00300   realLight(lv, lc);
00301   if(pic == NULL)
00302     pic = new t_data[Max*Max];
00303   ofstream out;
00304   int l;
00305   int m;
00306   int n;
00307   int L = Max; int M = Max; int N = Max;
00308   //  t_data sf=1.0/sqrt((t_data)(L*M*N)); // zero working array
00309   for(l=0;l<L;l++)
00310     for(m=0;m<M;m++) {
00311       Image[l][m].re =0.0;
00312       Image[l][m].im =0.0;
00313     }     
00314                     
00315   //  Color col;
00316   n=(int)ceil((t_data)Max/2.0);
00317   //  for(int yl=0;yl<9;yl++) {
00318   //    t_data c=dc[yl]*lc[yl]; // combined BRDF and light coef.
00319     for(l=0;l<L;l++) // loop over fft slice at origin 
00320       for(m=0;m<M;m++) // plane perpendicular to xy axis
00321         {
00322           //      Image[l][m].re +=c*Yarray[yl][indexOf(n,N,m,M,l,L)].re;
00323           //      Image[l][m].im +=c*Yarray[yl][indexOf(n,N,m,M,l,L)].im;
00324           //t_data* currSlice = slices[yl];
00325           //Image[l][m].re += c * currSlice[l * M + m];
00326           //Image[l][m].im += c * currSlice[l * M + m + 1];
00327         }
00328     //  }
00329     //fftwnd_one(rp2, &Image[0][0], NULL);
00330     fftwnd_one(rp2, (fftw_complex*)g_fSlice, NULL);
00331     for(l=0; l<L;l++)
00332       for(m=0;m<M;m++)
00333         {
00334           t_data val = g_fSlice[2*(l*M+m)];
00335           if(val > 1.0f) val = 1.0f;
00336           if(val < 0.0f) val = 0.0f;
00337           pic[l*M+m] = val;
00338         }
00339     return;
00340   
00341   t_data phaser; // undo the shift
00342   for(l=0;l<L;l++)
00343     for(m=0;m<M;m++) {
00344       phaser=(l+m) % 2 ==0 ? 1:-1;
00345       Image[l][m].re*=phaser;
00346       int val=(int)(clip(Image[l][m].re));
00347       //      col.set(val,val,val);
00348       //      pic->SetPixel(l,m,col);
00349       pic[l * M + m] = (t_data)val;
00350     }
00351   
00352   
00353   t_data maxval=-HUGE; 
00354   // should not need if I worked out the 
00355   // fft scaling correcly, I think.
00356   for(l=0;l<L;l++)  // scale dynamic range to max
00357     for(m=0;m<M;m++) {
00358       if(clip(Image[l][m].re)>maxval)
00359         maxval=clip(Image[l][m].re);
00360     }
00361   if(maxval>0)maxval=1.0/maxval*255;
00362   for(l=0;l<L;l++)
00363     for(m=0;m<M;m++) {
00364       int val=(int)(clip(Image[l][m].re)*maxval);
00365       //      col.set(val,val,val);
00366       //      pic->SetPixel(l,m,col);
00367       pic[l * M + m] = (t_data)val;///(t_data)maxval;
00368     }
00369 }  
00370 
00371 */

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