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

vuFourierCluster.cpp

Go to the documentation of this file.
00001 #include "vuFourierCluster.h"
00002 #include <math.h>
00003 #include <stddef.h>
00004 
00005 #define MAX_FILTER_SIZE 6 // max filter size
00006 
00007 #define ODD(x) ((x)&1)
00008 
00009 template <int S, class T>
00010 vuFourierCluster<S,T>::vuFourierCluster()
00011 {
00012   m_PlanWidth                = 0;
00013   m_PlanHeight               = 0;
00014   m_IsPreprocessed           = false;
00015   m_IsPreparedForInteractive = false;
00016   m_WeightVolume             = NULL;
00017   m_SliceFilter              = NULL;
00018   m_CacheVolume              = NULL;
00019 }
00020 
00021 template <int S, class T>
00022 vuFourierCluster<S,T>::~vuFourierCluster()
00023 {
00024   destroyPlan();
00025   CHECKNDELETE(m_WeightVolume);
00026   CHECKNDELETE(m_CacheVolume);
00027   // m_SliceFilter is an outside living object
00028 }
00029 
00030 template <int S, class T>
00031 bool vuFourierCluster<S,T>::isPreprocessed()
00032 {
00033   return m_IsPreprocessed;
00034 }
00035 
00036 template <int S, class T>
00037 bool vuFourierCluster<S,T>::isPreparedForInteractive()
00038 {
00039   return m_IsPreparedForInteractive;
00040 }
00041 
00042 template <int S, class T>
00043 void vuFourierCluster<S,T>::setNoInteractiveMode()
00044 {
00045   m_IsPreparedForInteractive = false;
00046   CHECKNDELETE(m_CacheVolume);
00047   CHECKNDELETE(m_WeightVolume);
00048   // ToDo: delete in case of true
00049 }
00050 
00051 template <int S, class T>
00052 void vuFourierCluster<S,T>::prepareForInteractive(dword width, dword height)
00053 {
00054   initializeVolume(width, height);
00055 
00056   CHECKNDELETE(m_CacheVolume);
00057   
00058   // initialize interactive cache
00059   {
00060     dword cnt     = m_XSize * m_YSize * m_ZSize * 2 *S;
00061     m_CacheVolume = new float[cnt];
00062   
00063     float *ptr = m_CacheVolume;
00064 
00065     for (dword i=0; i<cnt; i++) *(ptr++) = 0.0f;
00066   }
00067   
00068   wrapAndInitialize(1.0);
00069   
00070   m_IsPreparedForInteractive = true;
00071   m_IsPreprocessed           = true;
00072 }
00073 
00074 template <int S, class T>
00075 void vuFourierCluster<S,T>::addView(vuSphericView<S,T> *view)
00076 {
00077   if (m_SliceFilter == NULL) {
00078     cerr << "vuFourierCluster::addView() no slice filter set" << endl;
00079     return;
00080   }
00081 
00082   if (!m_IsPreparedForInteractive) {
00083     cerr << "vuFourierCluster::addView() is not prepared for interactive";
00084     cerr << " reconstruction!" << endl;
00085     return;
00086   }
00087  
00088   // add slice to m_CacheVolume
00089   {  
00090     float *tmpVolume = m_Volume;
00091 
00092     // use cache volume for slicing...
00093     m_Volume = m_CacheVolume;
00094   
00095     addViewToVolume(view);
00096 
00097     m_Volume = tmpVolume;
00098   }
00099   
00100   //weight m_CacheVolume with m_WeightVolume...
00101   {    
00102     dword sliceWidth, sliceHeight;
00103     calcSliceDimensions(view, sliceWidth, sliceHeight);  
00104     handleSlice(NULL, sliceWidth, sliceHeight, view);
00105   }
00106 }
00107 
00108 template <int S, class T>
00109 void vuFourierCluster<S,T>::preprocess(dword numOfViews,
00110                                        vuSphericView<S,T> **views,
00111                                        vuProgressHandler *handler)
00112 {
00113   if (numOfViews == 0) return;  
00114 
00115   dword width  = views[0]->getWidth();
00116   dword height = views[0]->getHeight();
00117 
00118   // progress handler variables...
00119   dword constFact    = 10;
00120   dword range        = 0;
00121   float progress     = 0.0f;
00122   float delta        = 0.0f;
00123 
00124   if (handler) {
00125     range = handler->getRange();
00126     delta = (float)range / (float)(numOfViews+2*constFact);    
00127   }
00128 
00129   if (handler) handler->update((int)progress, "Initialize Volume");
00130   initializeVolume(width, height);
00131   progress += delta*constFact;
00132 
00133   if (m_SliceFilter != NULL) {
00134     for (dword i=0; i<numOfViews; i++) {
00135 
00136       if (handler) {
00137         vuString msg("add View: ");
00138         msg += i;
00139         progress += delta;
00140         if (!handler->update((int)progress, msg)) break;
00141       }
00142       else
00143         cerr << "add slice " << i << "...";
00144 
00145       addViewToVolume(views[i]);
00146 
00147       if (!handler) cerr << "done." << endl;
00148     }
00149     if (handler) handler->update((int)progress, "Normalize Volume");
00150 
00151     normalizeVolume();
00152 
00153     if (handler) {
00154       progress += delta * constFact;
00155       handler->update((int)progress, "Initialize FFT");
00156     }
00157   }
00158   else {    
00159     cerr << "vuFourierCluster::preprocess(): m_SliceFilter is not set!"<<endl;
00160   }
00161 
00162   wrapAndInitialize(1.0);
00163   
00164   m_IsPreprocessed           = true;
00165   m_IsPreparedForInteractive = false;// not prepared for interactive reconstr.
00166 }
00167 
00168 template <int S, class T>
00169 void vuFourierCluster<S,T>::normalizeVolume()
00170 {
00171   dword size   = m_XSize * m_YSize * m_ZSize;
00172   float *v_ptr = m_Volume;
00173   float *w_ptr = m_WeightVolume;
00174   int   step   = S * 2;
00175   
00176   for (dword i=0; i<size; i++) {
00177     for (int j=0; j<step; j++) {
00178       if (*w_ptr == 0)
00179         *(v_ptr++) = 0.0f;
00180       else 
00181         *(v_ptr++) /= *w_ptr;
00182     }
00183     w_ptr++;
00184   }
00185   
00186   CHECKNDELETE(m_WeightVolume);
00187 }
00188 
00189 template <int S, class T>
00190 void vuFourierCluster<S,T>::setSliceFilter(vuSliceFilter *sliceFilter)
00191 {
00192   m_SliceFilter = sliceFilter;
00193 }
00194 
00195 template <int S, class T>
00196 vuSliceFilter *vuFourierCluster<S,T>::getSliceFilter()
00197 {
00198   return m_SliceFilter;
00199 }
00200 
00201 
00202 /* ************************************************************************* */
00203 /* *** private functions *************************************************** */
00204 /* ************************************************************************* */
00205 
00206 template <int S, class T>
00207 void vuFourierCluster<S,T>::initializeVolume(dword width, dword height)
00208 {
00209   // don't change the next line, it would have strong effect on the algorithm! 
00210   dword dim = computeDimensions(width,height,height, M_SQRT2,MAX_FILTER_SIZE);
00211   dword cnt = m_XSize * m_YSize * m_ZSize;
00212 
00213   if (m_Volume == NULL || m_XSize != dim || m_YSize != dim || m_ZSize != dim) {
00214     m_XSize = m_YSize = m_ZSize = dim;
00215     cnt     = m_XSize * m_YSize * m_ZSize;
00216 
00217     CHECKNDELETE(m_Volume);    
00218     m_Volume       = new float[2*cnt*S];
00219   }
00220 
00221   CHECKNDELETE(m_WeightVolume);
00222   CHECKNDELETE(m_CacheVolume); // just in case ;-)
00223 
00224   m_WeightVolume = new float[cnt];
00225 
00226   float *v_ptr = m_Volume;
00227   float *w_ptr = m_WeightVolume;
00228 
00229   for (dword i=0; i<cnt; i++) *(w_ptr++) = 0.0f;
00230   cnt = cnt * S * 2;
00231   for (dword i=0; i<cnt; i++) *(v_ptr++) = 0.0f;
00232 }
00233 
00234 template <int S, class T>
00235 void vuFourierCluster<S,T>::transformSlice(float *slice, 
00236                                            dword width, dword height)
00237 {
00238   fftwnd(m_Plan, S, (fftw_complex*)slice, S, 1, NULL, 0, 0);
00239   dword count  = width * height * 2 * S;
00240   float size_f = (float) width * height;
00241   float *ptr   = slice;
00242   for (dword i = 0; i < count; ++i) *(ptr++) /= size_f;
00243 }
00244 
00245 template <int S, class T>
00246 void vuFourierCluster<S,T>::ensurePlan(dword width, dword height)
00247 {
00248   if ((m_PlanWidth != width) || (m_PlanHeight != height)) {
00249     destroyPlan();
00250       
00251     m_Plan  = fftw2d_create_plan(width, height,
00252                                  FFTW_FORWARD, FFTW_MEASURE |
00253                                  FFTW_IN_PLACE | FFTW_USE_WISDOM);
00254     m_PlanWidth  = width;
00255     m_PlanHeight = height;
00256   }
00257 }
00258 
00259 template <int S, class T>
00260 void vuFourierCluster<S,T>::destroyPlan()
00261 {
00262   if ((m_PlanWidth) != 0 || (m_PlanHeight != 0)) {
00263     fftwnd_destroy_plan(m_Plan);
00264     m_PlanWidth = m_PlanHeight = 0;
00265   }
00266 }
00267 
00268 template <int S, class T>
00269 void vuFourierCluster<S,T>::calcSliceDimensions(vuSphericView<S,T> *view,
00270                                                 dword &sliceWidth,
00271                                                 dword &sliceHeight)
00272 {
00273   sliceWidth  = (dword)((m_XSize - MAX_FILTER_SIZE) / M_SQRT2);
00274   sliceHeight = (dword)((m_YSize - MAX_FILTER_SIZE) / M_SQRT2);
00275   int deltaWidth  = view->getWidth()  - sliceWidth;
00276   int deltaHeight = view->getHeight() - sliceHeight;
00277   
00278   if (deltaWidth  == 0) deltaWidth  = 1;
00279   if (deltaHeight == 0) deltaHeight = 1;
00280 
00281   if (abs(deltaWidth)  > 1) {
00282     cerr << "vuFourierCluster::addViewToVolume: sliceWidth does not fit"<<endl;
00283     return;
00284   }
00285   
00286   if (abs(deltaHeight) > 1) {
00287     cerr <<"vuFourierCluster::addViewToVolume: sliceHeight does not fit"<<endl;
00288     return;
00289   }
00290   
00291   if (ODD(sliceWidth))  sliceWidth  += deltaWidth;
00292   if (ODD(sliceHeight)) sliceHeight += deltaHeight;
00293 
00294   if ((sliceWidth != view->getWidth()) || (sliceHeight != view->getHeight())) {
00295     cerr << "vuFourierCluster::calcSliceDimensions:slice and view have not the"
00296             " same size" << endl;
00297     throw("vuFourierCluster::calcSliceDimensions: slice and view have not the"
00298           " same size");
00299   }
00300 }
00301 
00302 template <int S, class T>
00303 void vuFourierCluster<S,T>::addViewToVolume(vuSphericView<S,T> *view)
00304 {
00305   dword sliceWidth, sliceHeight;
00306 
00307   // ToDo: simplify calcSliceDimension, e.g. sliceWidth = view->getWidth()
00308   calcSliceDimensions(view, sliceWidth, sliceHeight);
00309   
00310   float *slice = new float [sliceWidth * sliceHeight * 2 * S];
00311 
00312   ensurePlan(sliceWidth, sliceHeight);
00313 
00314   // copy view to slice
00315   T     *v_ptr = view->getMap()->getBuffer();
00316   float *s_ptr = slice;
00317   
00318   for (dword j=0; j< sliceHeight; j++) {
00319     for (dword i = 0; i < sliceWidth; ++i) {
00320       for (int k=0; k<S; k++) {
00321         *(s_ptr++) = (float)*(v_ptr++);
00322         *(s_ptr++) = 0.0f;
00323       }
00324     }
00325   }
00326 
00327   shift2D(slice, sliceWidth, sliceHeight);
00328   transformSlice(slice, sliceWidth, sliceHeight);
00329   shift2D(slice, sliceWidth, sliceHeight);
00330 
00331   //premultiplySlice(slice);
00332   handleSlice(slice, sliceWidth, sliceHeight, view); // add slice to volume
00333 
00334   CHECKNDELETE(slice);
00335 }
00336 
00337 template <int S, class T>
00338 void vuFourierCluster<S,T>::handleSlice(float *slice,
00339                                         dword sliceWidth,
00340                                         dword sliceHeight,
00341                                         vuSphericView<S,T> *view)
00342 {
00343   vuVector lookAt = view->getLookFrom() * -1;
00344   vuVector yStep  = view->getUp();
00345   vuVector xStep  = vuVector(1,0,0);
00346 
00347   calcViewVectors(lookAt, yStep, xStep);
00348 
00349   xStep.makeUnit();
00350   yStep.makeUnit();
00351 
00352   dword innerXSize = sliceWidth;
00353   dword innerYSize = sliceHeight;
00354 
00355   vuVector origin(m_XSize, m_YSize, m_ZSize);
00356   origin *= 0.5;
00357 
00358   origin += xStep*(innerXSize*-0.5f);
00359   origin += yStep*(innerYSize*-0.5f);
00360 
00361   dword XStart = 0;
00362   dword YStart = 0;
00363   dword XStop  = innerXSize;
00364   dword YStop  = innerYSize;
00365 
00366   vuVector pos = origin;
00367 
00368   if (m_SliceFilter->getLowPassFactor() < 1.0) {
00369     float lowPass = m_SliceFilter->getLowPassFactor();
00370     lowPass = 1.0 - sqrt(lowPass);
00371 
00372     dword cutOffX = innerXSize / 2;
00373     dword cutOffY = innerYSize / 2;
00374 
00375     cutOffX = dword(cutOffX * lowPass);
00376     cutOffY = dword(cutOffY * lowPass);
00377 
00378     XStart += cutOffX;
00379     YStart += cutOffY;
00380     XStop  -= cutOffX;
00381     YStop  -= cutOffY;
00382 
00383     pos += xStep * cutOffX;
00384     pos += yStep * cutOffY;
00385   }
00386 
00387   dword width   = m_SliceFilter->getWidth();
00388   mc_HalfWidth  = (float)width / 2;          // float
00389   mc_High       = (width/2);                 // int
00390   mc_Low        = -mc_High + (1- width % 2); // int
00391   mc_IsWidthOdd = (ODD(width));              // bool
00392 
00393   //mc_HalfWidth *= 1.5;
00394 
00395   KernelCallback kernelCB;
00396 
00397   if (slice == NULL)
00398     kernelCB = &vuFourierCluster<S,T>::doWeighting;
00399   else if (m_SliceFilter->getKind() == 0)
00400     kernelCB = &vuFourierCluster<S,T>::doFilteringSeparable;
00401   else
00402     kernelCB = &vuFourierCluster<S,T>::doFilteringSpheric;
00403 
00404   for (dword j = YStart; j < YStop; ++j) {
00405     vuVector p = pos;
00406     for (dword i = XStart; i < XStop; ++i) {
00407       (this->*kernelCB)(p, &slice[(j*innerXSize+i) * 2 * S]);
00408       p += xStep;
00409     }
00410     pos += yStep;
00411   }
00412 }
00413 
00414 template <int S, class T>
00415 inline void vuFourierCluster<S,T>::doFilteringSeparable(vuVector& pos, float *value)
00416 {
00417   dword center[3]; // center in Volume (grid coordinate)
00418   float POS[3];    // 
00419 
00420   if (mc_IsWidthOdd) {
00421     center[0] = (dword)(pos[0]+0.5f);
00422     center[1] = (dword)(pos[1]+0.5f);
00423     center[2] = (dword)(pos[2]+0.5f);
00424   }
00425   else {
00426     center[0] = (dword)pos[0];
00427     center[1] = (dword)pos[1];
00428     center[2] = (dword)pos[2];
00429   }
00430 
00431   POS[0] = pos[0] - center[0];
00432   POS[1] = pos[1] - center[1];
00433   POS[2] = pos[2] - center[2];
00434 
00435   for (int k=mc_Low; k<=mc_High; k++) {
00436     register float facZ = (POS[2]-k) / mc_HalfWidth;
00437 
00438     facZ = m_SliceFilter->getWeight(facZ);
00439 
00440     for (int l=mc_Low; l<=mc_High; l++) {
00441       register float facYZ = (POS[1]-l) / mc_HalfWidth;
00442 
00443       facYZ = facZ *m_SliceFilter->getWeight(facYZ);
00444 
00445       for (int m=mc_Low; m<=mc_High; m++) {
00446         register dword idx   = wcoord(center[0]+m,center[1]+l,center[2]+k);
00447         register float fac   = (POS[0]-m)/mc_HalfWidth;
00448 
00449         fac = facYZ * m_SliceFilter->getWeight(fac);
00450 
00451         if (fac < 0.0) cerr << fac << " "; // this should never happen
00452 
00453         m_WeightVolume[idx] += fac;
00454         idx *= (S+S);
00455 
00456         float *v_ptr = &m_Volume[idx];
00457         float *s_ptr = value;
00458         for (int i=0; i<S; i++) {
00459           *(v_ptr++) += fac * *(s_ptr++);
00460           *(v_ptr++) += fac * *(s_ptr++);
00461         }
00462       }
00463     }
00464   }
00465 }
00466 
00467 template <int S, class T>
00468 inline void vuFourierCluster<S,T>::doFilteringSpheric(vuVector& pos, float *value)
00469 {
00470   dword center[3]; // center in Volume (grid coordinate)
00471   float POS[3];    // 
00472 
00473   if (mc_IsWidthOdd) {
00474     center[0] = (dword)(pos[0]+0.5f);
00475     center[1] = (dword)(pos[1]+0.5f);
00476     center[2] = (dword)(pos[2]+0.5f);
00477   }
00478   else {
00479     center[0] = (dword)pos[0];
00480     center[1] = (dword)pos[1];
00481     center[2] = (dword)pos[2];
00482   }
00483 
00484   POS[0] = pos[0] - center[0];
00485   POS[1] = pos[1] - center[1];
00486   POS[2] = pos[2] - center[2];
00487 
00488   for (int k=mc_Low; k<=mc_High; k++) {
00489     for (int l=mc_Low; l<=mc_High; l++) {
00490       for (int m=mc_Low; m<=mc_High; m++) {
00491         register float xx    = POS[0] - m;
00492         register float yy    = POS[1] - l;
00493         register float zz    = POS[2] - k;
00494         register float fac   = sqrt(xx*xx+yy*yy+zz*zz) / mc_HalfWidth;
00495         register dword idx   = wcoord(center[0]+m,center[1]+l,center[2]+k);
00496 
00497         fac = m_SliceFilter->getWeight(fac);
00498 
00499         if (fac < 0.0) cerr << fac << " "; // this should never happen
00500 
00501         m_WeightVolume[idx] += fac;
00502         idx *= (S+S);
00503 
00504         float *v_ptr = &m_Volume[idx];
00505         float *s_ptr = value;
00506         for (int i=0; i<S; i++) {
00507           *(v_ptr++) += fac * *(s_ptr++);
00508           *(v_ptr++) += fac * *(s_ptr++);
00509         }
00510       }
00511     }
00512   }
00513 }
00514 
00515 // --------------------------------------------------------------------------
00516 // --- private - interactive reconstruction - weighting - functions ---------
00517 // --------------------------------------------------------------------------
00518 
00519 // ?_? to be done:  (channels)
00520 template <int S, class T>
00521 void vuFourierCluster<S,T>::weightView(vuSphericView<S,T> *view)
00522 {
00523   vuVector lookAt = view->getLookFrom();
00524   vuVector up     = view->getUp();
00525   vuVector right  = vuVector(1,0,0);
00526 
00527   calcViewVectors(lookAt, up, right);
00528 
00529   // ----------------------------------------------------------------------
00530 
00531   vuVector xStep  = right.makeUnit(); //  * (float)(m_XSize/m_XSize);
00532   vuVector yStep  = up.makeUnit();    //  * (float)(m_YSize/m_YSize);
00533 
00534   vuVector origin(m_XSize, m_YSize, m_ZSize);
00535   origin *= 0.5;
00536   
00537   dword innerXSize = (dword)((m_XSize - MAX_FILTER_SIZE) / M_SQRT2);
00538   dword innerYSize = (dword)((m_YSize - MAX_FILTER_SIZE) / M_SQRT2);
00539 
00540   if (ODD(innerXSize)) innerXSize += 1;
00541   if (ODD(innerYSize)) innerYSize += 1;
00542 
00543   origin += xStep*(innerXSize*-0.5f);
00544   origin += yStep*(innerYSize*-0.5f);
00545 
00546   dword XStart = (m_XSize - innerXSize) / 2;
00547   dword YStart = (m_YSize - innerYSize) / 2;
00548   dword XStop  = XStart + innerXSize;
00549   dword YStop  = YStart + innerYSize;
00550 
00551   vuVector pos = origin;
00552 
00553   dword width       = m_SliceFilter->getWidth();
00554   mc_High           = (width/2);                 // int
00555   mc_Low            = -mc_High + (1- width % 2); // int
00556   mc_IsWidthOdd     = (ODD(width));              // bool
00557 
00558   for (dword j = YStart; j < YStop; ++j) {
00559     vuVector p = pos;
00560     for (dword i = XStart; i < XStop; ++i) {
00561       doWeighting(p, NULL);
00562       p += xStep;
00563     }
00564     pos += yStep;
00565   }
00566 }
00567 
00568 // ?_?: to be done (channels)
00569 template <int S, class T>
00570 inline void vuFourierCluster<S,T>::doWeighting(vuVector& pos, float *value)
00571 {
00572   dword center[3]; // center in Volume (grid coordinate)
00573 
00574   if (mc_IsWidthOdd) {
00575     center[0] = (dword)(pos[0]+0.5f);
00576     center[1] = (dword)(pos[1]+0.5f);
00577     center[2] = (dword)(pos[2]+0.5f);
00578   }
00579   else {
00580     center[0] = (dword)pos[0];
00581     center[1] = (dword)pos[1];
00582     center[2] = (dword)pos[2];
00583   }
00584 
00585   for (int k=mc_Low; k<=mc_High; k++) {
00586     for (int l=mc_Low; l<=mc_High; l++) {
00587       for (int m=mc_Low; m<=mc_High; m++) {
00588         dword idx = vcoord(center[0]+m,center[1]+l,center[2]+k);
00589 
00590         m_Volume[idx]   = m_CacheVolume[idx]   / m_WeightVolume[idx/2];
00591         m_Volume[idx+1] = m_CacheVolume[idx+1] / m_WeightVolume[idx/2];
00592       }
00593     }
00594   }
00595 }
00596 
00597 // --------------------------------------------------------------------------
00598 // --- private - premultiplication ------------------------------------------
00599 // --------------------------------------------------------------------------
00600 
00601 // ?_? to be done: (channels)
00602 template <int S, class T>
00603 void vuFourierCluster<S,T>::premultiplySlice(vuFixelMap2F *slice)
00604 {
00605   dword width      = slice->getWidth();
00606   dword height     = slice->getHeight();
00607   float halfWidth  = (float)width  / 2;
00608   float halfHeight = (float)height / 2;
00609   float *ptr       = slice->getBuffer();
00610 
00611   for (dword j=0; j<height; j++) {
00612     float zFac = m_SliceFilter->getTransformedWeight(j-halfHeight);
00613 
00614     zFac *= m_SliceFilter->getWidth();
00615     
00616     for (dword i=0; i<width; i++) {
00617       float fac = zFac * m_SliceFilter->getTransformedWeight(i-halfWidth);
00618         
00619       fac *= m_SliceFilter->getWidth();
00620       
00621       // cerr << "fac=" << fac << endl;
00622       
00623       ptr[(j*width+i)*2] /= fac;
00624     }
00625   }
00626 }

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