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

vuFourierVolume.cpp

Go to the documentation of this file.
00001 #include "vuFourierVolume.h"
00002 
00003 #include <stdio.h>
00004 #include "vuMatrix.h"
00005 #include "math.h"
00006 
00007 #define WISDOM "fvr_fftw.wis"
00008 
00009 #define ODD(x) ((x)&1)
00010 
00011 template <int S>
00012 vuFourierVolume<S>::vuFourierVolume() : m_XStep(1.0f, 0.0f, 0.0f),
00013                                         m_YStep(0.0f, 1.0f, 0.0f),
00014                                         m_XAxis(1.0f, 0.0f, 0.0f),
00015                                         m_YAxis(0.0f, 1.0f, 0.0f),
00016                                         m_ZAxis(0.0f, 0.0f, 1.0f)
00017 {
00018   m_Volume        = NULL;
00019   m_Slice         = NULL;
00020   m_Filter        = NULL;
00021 
00022   m_Wrap          = 3; // largest filter.width()/2
00023 
00024   m_ImageStep     = 0;
00025   m_SlicePtr      = NULL;
00026 
00027   loadWisdom();
00028 
00029   m_SliceXSize    = 0;
00030   m_SliceYSize    = 0;
00031 
00032   m_RenderMethod  = 0; // spatial image
00033   m_Scale         = 1.0f;
00034   m_Bias          = 0.0f;
00035   m_UnscaledImage = NULL;
00036   m_IsReRendering = true;
00037   for (int i=0; i<S; i++) m_Channels[i] = true;
00038 }
00039 
00040 template <int S>
00041 vuFourierVolume<S>::vuFourierVolume(vuFourierVolume<S>& inst)
00042 {
00043   cerr << "vuFourierVolume::copyConstructor is not implemented!" << endl;
00044   throw "vuFourierVolume::copyConstructor is not implemented!";
00045 }
00046 
00047 template <int S>
00048 vuFourierVolume<S>::~vuFourierVolume()
00049 {
00050   CHECKNDELETE(m_Volume);
00051   CHECKNDELETE(m_Slice);
00052   CHECKNDELETE(m_UnscaledImage);
00053   destroyTransform2D();
00054   saveWisdom();
00055 }
00056 
00057 template <int S>
00058 vuFourierVolume<S>& vuFourierVolume<S>::operator=(vuFourierVolume<S>& rhs)
00059 {
00060   if (this != &rhs)
00061   {
00062   }
00063   return *this;
00064 }
00065 
00066 template <int S>
00067 dword vuFourierVolume<S>::computeDimensions(dword XSize, dword YSize,
00068                                             dword ZSize, float m_pad,
00069                                             dword a_pad)
00070 {
00071   dword dim = (XSize > YSize) ? XSize : YSize;
00072   dim = (dim > ZSize) ? dim : ZSize;
00073   dim = (dword) ceil((float)dim * m_pad) + a_pad;
00074   dim = (ODD(dim)) ? dim + 1 : dim;
00075 
00076   return dim;
00077 }
00078 
00079 template <int S>
00080 void vuFourierVolume<S>::setOversampling(float s) {
00081   // don't allow undersampling
00082   if (s < 1.0) s = 1.0;
00083 
00084   dword sliceXSize = (dword) ceil(m_ImageXSize * s);
00085   dword sliceYSize = (dword) ceil(m_ImageYSize * s);
00086   sliceXSize = (ODD(sliceXSize)) ? sliceXSize + 1 : sliceXSize;
00087   sliceYSize = (ODD(sliceYSize)) ? sliceYSize + 1 : sliceYSize;
00088 
00089   // initialize m_Slice and fftw
00090   if (sliceXSize != m_SliceXSize || sliceYSize != m_SliceYSize || !m_Slice) {
00091     destroyTransform2D();
00092 
00093     m_SliceXSize = sliceXSize;
00094     m_SliceYSize = sliceYSize;
00095 
00096     if (m_Slice) delete [] m_Slice;
00097     m_Slice = new float[m_SliceXSize * m_SliceYSize * 2 * S];
00098   
00099     if ((sliceXSize > 0) && (sliceYSize > 0)) {
00100       cerr << "init plan ";
00101 
00102       m_Plan2D = fftw2d_create_plan_specific(sliceXSize, sliceYSize,
00103                                              FFTW_BACKWARD,
00104                                              FFTW_MEASURE | FFTW_IN_PLACE |
00105                                              FFTW_USE_WISDOM,
00106                                              (fftw_complex*)m_Slice, S,
00107                                              NULL, 0);
00108       cerr << " done." << endl;
00109     }
00110   }
00111 
00112   m_XStep = m_XAxis * m_ImageXSize / m_SliceXSize;
00113   m_YStep = m_YAxis * m_ImageYSize / m_SliceYSize;
00114 
00115   m_ImageStep = (m_SliceXSize - m_ImageXSize) * 2 * S; // ?_?
00116 
00117   cerr << "Slice size: " << m_SliceXSize << ' ' << m_SliceYSize << endl;
00118 
00119   m_SlicePtr = &m_Slice[scoord( (m_SliceXSize - m_ImageXSize) / 2,
00120                                 (m_SliceYSize - m_ImageYSize) / 2)];
00121 
00122   m_InnerXSize = (dword) floor((float) m_SliceXSize / M_SQRT2);
00123   m_InnerYSize = (dword) floor((float) m_SliceYSize / M_SQRT2);
00124   m_InnerXSize = (ODD(m_InnerXSize)) ? m_InnerXSize + 1 : m_InnerXSize;
00125   m_InnerYSize = (ODD(m_InnerYSize)) ? m_InnerYSize + 1 : m_InnerYSize;
00126 }
00127 
00128 template <int S>
00129 void vuFourierVolume<S>::setWrap(dword wrap)
00130 {
00131   if (m_Filter) {
00132     if ((m_Filter->getWidth() / 2) > wrap)
00133       m_Wrap = m_Filter->getWidth()/2;
00134     else
00135       m_Wrap = wrap;
00136   }
00137   else
00138     m_Wrap = wrap;
00139 }
00140 
00141 template <int S>
00142 dword vuFourierVolume<S>::getWrap()
00143 {
00144   return m_Wrap;
00145 }
00146 
00147 template <int S>
00148 void vuFourierVolume<S>::setFilter(vuFourierFilter* filter)
00149 {
00150   if (m_Filter != filter) m_IsReRendering = true;
00151   m_Filter = filter;
00152 }
00153 
00154 template <int S>
00155 vuFourierFilter *vuFourierVolume<S>::getFilter()
00156 {
00157   return m_Filter;
00158 }
00159 
00160 
00161 template <int S>
00162 void vuFourierVolume<S>::setCamera(const vuCamera *camera)
00163 {
00164   if (!isInitialized()) return;
00165 
00166   setViewVectors(camera->getLookAtVector(), 
00167                  camera->getUpVector(),
00168                  camera->getRightVector());
00169 }
00170 
00171 template <int S>
00172 void vuFourierVolume<S>::setViewVectors(const vuVector& view,
00173                                         const vuVector& up,
00174                                         const vuVector& right)
00175 {
00176   if (!isInitialized()) return;
00177 
00178   m_XAxis  = right;
00179   m_YAxis  = up;
00180   m_ZAxis  = view;
00181 
00182   m_XStep  = right;
00183   m_YStep  = up;
00184   
00185   m_XStep.makeUnit();
00186   m_YStep.makeUnit();
00187 
00188   m_XStep *= (float)(m_ImageXSize / m_SliceXSize);
00189   m_YStep *= (float)(m_ImageYSize / m_SliceYSize);
00190 
00191   m_Origin  = m_XAxis*(m_SliceXSize*-0.5f);
00192   m_Origin += m_YAxis*(m_SliceYSize*-0.5f);
00193   m_IsReRendering = true;
00194 }
00195 
00196 template <int S>
00197 void vuFourierVolume<S>::calcViewVectors(vuVector& lookAt,
00198                                          vuVector& up,
00199                                          vuVector& right)
00200 {
00201   lookAt.makeUnit();
00202 
00203   if (fabs(lookAt.dot(up)) > 0.999) up = vuVector(0,0,1).makeUnit();
00204 
00205   up = (up - (lookAt.dot(up) * lookAt)).makeUnit();
00206   right  = up.cross(lookAt).makeUnit() * -1;
00207 }
00208 
00209 template <int S>
00210 dword vuFourierVolume<S>::getImageWidth() const
00211 {
00212   return m_ImageXSize;
00213 }
00214 
00215 template <int S>
00216 dword vuFourierVolume<S>::getImageHeight() const
00217 {
00218   return m_ImageYSize;
00219 }
00220 
00221 template <int S>
00222 dword vuFourierVolume<S>::getSliceWidth() const
00223 {
00224   return m_SliceXSize;
00225 }
00226 
00227 template <int S>
00228 dword vuFourierVolume<S>::getSliceHeight() const
00229 {
00230   return m_SliceYSize;
00231 }
00232 
00233 template <int S>
00234 float *vuFourierVolume<S>::getSliceData() const
00235 {
00236   return m_SlicePtr;
00237 }
00238 
00239 template <int S>
00240 dword vuFourierVolume<S>::getXSize()
00241 {
00242   return m_XSize;
00243 }
00244 
00245 template <int S>
00246 dword vuFourierVolume<S>::getYSize()
00247 {
00248   return m_YSize;
00249 }
00250 
00251 template <int S>
00252 dword vuFourierVolume<S>::getZSize()
00253 {
00254   return m_ZSize;
00255 }
00256 
00257 template <int S>
00258 void vuFourierVolume<S>::clearSlice()
00259 {
00260   float* s_ptr = m_Slice;
00261   dword  count = m_SliceXSize * m_SliceYSize * 2 * S;
00262   
00263   for (dword i = 0; i < count; i++) *(s_ptr++) = 0.0f;
00264 }
00265 
00266 #define CONVOLVE_NO_MOD                                                   \
00267 {                                                                         \
00268   int x = (int) p[0];                                                     \
00269   int y = (int) p[1];                                                     \
00270   int z = (int) p[2];                                                     \
00271                                                                           \
00272   t[0] = p[0] - x;                                                        \
00273   t[1] = p[1] - y;                                                        \
00274   t[2] = p[2] - z;                                                        \
00275                                                                           \
00276   x = x - XSize + m_Wrap;                                                 \
00277   y = y - YSize + m_Wrap;                                                 \
00278   z = z - ZSize + m_Wrap;                                                 \
00279                                                                           \
00280   register float* fweight = m_Filter->getWeights(t);                      \
00281   register float* v_ptr   = &m_Volume[vcoord(x + f2, y + f2, z + f2)];    \
00282                                                                           \
00283   for (dword l = 0; l < f_len; l++) {                                     \
00284     for (dword m = 0; m < f_len; m++) {                                   \
00285       for (dword n = 0; n < f_len; n++) {                                 \
00286         register float tmp     = *(fweight++);                            \
00287         register float *ss_ptr = s_ptr;                                   \
00288         register float *vv_ptr = v_ptr;                                   \
00289         for (int k=0; k<S; k++) {                                         \
00290           *(ss_ptr++) += *(vv_ptr++) * tmp;                               \
00291           *(ss_ptr++) += *(vv_ptr++) * tmp;                               \
00292         }                                                                 \
00293         v_ptr -= (S+S);                                                   \
00294       }                                                                   \
00295       v_ptr += fXStep;                                                    \
00296     }                                                                     \
00297     v_ptr += fYStep;                                                      \
00298   }                                                                       \
00299                                                                           \
00300   s_ptr += (S+S);                                                         \
00301   p     += m_XStep;                                                       \
00302 }
00303 
00304 #define CONVOLVE                                                          \
00305 {                                                                         \
00306   int x = (int) p[0];                                                     \
00307   int y = (int) p[1];                                                     \
00308   int z = (int) p[2];                                                     \
00309                                                                           \
00310   t[0] = p[0] - x;                                                        \
00311   t[1] = p[1] - y;                                                        \
00312   t[2] = p[2] - z;                                                        \
00313                                                                           \
00314   x = x % XSize + m_Wrap;                                                 \
00315   y = y % YSize + m_Wrap;                                                 \
00316   z = z % ZSize + m_Wrap;                                                 \
00317                                                                           \
00318   register float* fweight = m_Filter->getWeights(t);                      \
00319   register float* v_ptr   = &m_Volume[vcoord(x + f2, y + f2, z + f2)];    \
00320                                                                           \
00321   for (dword l = 0; l < f_len; l++) {                                     \
00322     for (dword m = 0; m < f_len; m++) {                                   \
00323       for (dword n = 0; n < f_len; n++) {                                 \
00324         register float tmp     = *(fweight++);                            \
00325         register float *ss_ptr = s_ptr;                                   \
00326         register float *vv_ptr = v_ptr;                                   \
00327         for (int k=0; k<S; k++) {                                         \
00328           *(ss_ptr++) += *(vv_ptr++) * tmp;                               \
00329           *(ss_ptr++) += *(vv_ptr++) * tmp;                               \
00330         }                                                                 \
00331         v_ptr -= (S+S);                                                   \
00332       }                                                                   \
00333       v_ptr += fXStep;                                                    \
00334     }                                                                     \
00335     v_ptr += fYStep;                                                      \
00336   }                                                                       \
00337                                                                           \
00338   s_ptr += (S+S);                                                         \
00339   p     += m_XStep;                                                       \
00340 }
00341 
00342 template <int S>
00343 void vuFourierVolume<S>::interpolateSlice(dword x_stop, dword y_stop,
00344                                           dword x_pass, dword y_pass)
00345 {
00346   if (m_Filter == NULL) {
00347     cerr << "Warning: interpolateSlice(refine): no filter set!" << endl;
00348     return; 
00349   }
00350 
00351   vuVector t;
00352   dword f_len = m_Filter->getWidth();
00353   int f2 = f_len / 2;
00354 
00355   x_stop *= 2;
00356   y_stop *= 2;
00357   dword stop_l = (m_SliceXSize - x_stop - 2 * x_pass) / 2;
00358   dword stop_r =  m_SliceXSize - x_stop - 2 * x_pass - stop_l;
00359   dword stop_b = (m_SliceYSize - y_stop - 2 * y_pass) / 2;
00360   dword x_hi = x_stop + 2 * x_pass;
00361   dword s_step = stop_l + stop_r;
00362 
00363   dword XSize = m_XSize - m_Wrap * 2;
00364   dword YSize = m_YSize - m_Wrap * 2;
00365   dword ZSize = m_ZSize - m_Wrap * 2;
00366 
00367   dword fXStep = (f_len - m_XSize) * 2 * S;
00368   dword fYStep = (f_len - m_YSize) * m_XSize * 2 * S;
00369 
00370   vuVector q = m_Origin;
00371   q[0] += XSize / 2 + XSize + stop_l * m_XStep[0] + stop_b * m_YStep[0];
00372   q[1] += YSize / 2 + YSize + stop_l * m_XStep[1] + stop_b * m_YStep[1];
00373   q[2] += ZSize / 2 + ZSize + stop_l * m_XStep[2] + stop_b * m_YStep[2];
00374   
00375   clearSlice();
00376   
00377   float* s_ptr = &m_Slice[scoord(stop_l, stop_b)];
00378 
00379   for (dword j = 0; j < y_pass; ++j) {
00380     vuVector p = q;
00381     for (dword i = 0; i < x_hi; ++i) CONVOLVE;
00382     s_ptr += s_step * 2 * S;
00383     q += m_YStep;
00384   }
00385 
00386   for (dword j = 0; j < y_stop; ++j) {
00387     vuVector p = q;
00388     for (dword i = 0; i < x_pass; ++i) CONVOLVE;
00389 
00390     s_ptr += x_stop * 2 * S;
00391     p[0] += x_stop * m_XStep[0];
00392     p[1] += x_stop * m_XStep[1];
00393     p[2] += x_stop * m_XStep[2];
00394 
00395     for (dword i = 0; i < x_pass; ++i) CONVOLVE;
00396     s_ptr += s_step * 2 * S;
00397     q += m_YStep;
00398   }
00399 
00400   for (dword j = 0; j < y_pass; ++j) {
00401     vuVector p = q;
00402     for (dword i = 0; i < x_hi; ++i) CONVOLVE;
00403     s_ptr += s_step * 2 * S;
00404     q += m_YStep;
00405   }
00406 }
00407 
00408 template <int S>
00409 void vuFourierVolume<S>::interpolateSlice()
00410 {
00411   if (m_Filter == NULL) {
00412     cerr << "Warning: interpolateSlice(): no filter set!" << endl;
00413     return; 
00414   }
00415 
00416   dword i, j;
00417   vuVector t;
00418   dword f_len = m_Filter->getWidth();
00419   int f2 = f_len / 2;
00420 
00421   dword x_pass = (m_SliceXSize - m_InnerXSize) / 2;
00422   dword y_pass = (m_SliceYSize - m_InnerYSize) / 2;
00423 
00424   dword XSize = m_XSize - m_Wrap * 2;
00425   dword YSize = m_YSize - m_Wrap * 2;
00426   dword ZSize = m_ZSize - m_Wrap * 2;
00427 
00428   dword fXStep = (f_len - m_XSize) * 2 * S;
00429   dword fYStep = (f_len - m_YSize) * m_XSize * 2 * S;
00430 
00431   vuVector q = m_Origin;
00432   q[0] += XSize / 2 + XSize;
00433   q[1] += YSize / 2 + YSize;
00434   q[2] += ZSize / 2 + ZSize;
00435 
00436   clearSlice();
00437 
00438   float* s_ptr = m_Slice;
00439 
00440   for (j = 0; j < y_pass; ++j) {
00441     vuVector p = q;
00442     for (i = 0; i < m_SliceXSize; ++i) CONVOLVE;
00443     q += m_YStep;
00444   }
00445 
00446   for (j = 0; j < m_InnerYSize; ++j) {
00447     vuVector p = q;
00448     for (i = 0; i < x_pass; ++i) CONVOLVE;
00449     for (i = 0; i < m_InnerXSize; ++i) CONVOLVE_NO_MOD;
00450     for (i = 0; i < x_pass; ++i) CONVOLVE;
00451     q += m_YStep;
00452   }
00453 
00454   for (j = 0; j < y_pass; ++j) {
00455     vuVector p = q;
00456     for (i = 0; i < m_SliceXSize; ++i) CONVOLVE;
00457     q += m_YStep;
00458   }
00459 }
00460 
00461 template <int S>
00462 void vuFourierVolume<S>::computeSlice()
00463 {
00464   interpolateSlice();
00465   transformSlice();
00466 }
00467 
00468 template <int S>
00469 void vuFourierVolume<S>::transformSlice()
00470 {
00471   shift2D(m_Slice, m_SliceXSize, m_SliceYSize);
00472   transform2D(m_Slice);
00473   shift2D(m_Slice, m_SliceXSize, m_SliceYSize);  
00474 }
00475 
00476 
00477 template <int S>
00478 void vuFourierVolume<S>::computeSlice(dword x_stop, dword y_stop,
00479                                       dword x_pass, dword y_pass)
00480 {
00481   interpolateSlice(x_stop, y_stop, x_pass, y_pass);
00482   
00483   shift2D(m_Slice, m_SliceXSize, m_SliceYSize);
00484   transform2D(m_Slice);
00485   shift2D(m_Slice, m_SliceXSize, m_SliceYSize);
00486 }
00487 
00488 template <int S>
00489 bool vuFourierVolume<S>::isInitialized()
00490 {
00491   return (m_SliceXSize > 0 || m_SliceYSize > 0);
00492 }
00493 
00494 /* ************************************************************************* */
00495 /* *** private functions *************************************************** */
00496 /* ************************************************************************* */
00497 
00498 template <int S>
00499 void vuFourierVolume<S>::wrapAndInitialize(float overSampling)
00500 {
00501   wrapVolume(m_Volume);
00502 
00503   dword XSize = m_XSize - m_Wrap * 2;
00504   dword YSize = m_YSize - m_Wrap * 2;
00505 
00506   m_Origin[0] = (float)(XSize / 2) * -1.0f;
00507   m_Origin[1] = (float)(YSize / 2) * -1.0f;
00508   m_Origin[2] = 0.0f;
00509 
00510   m_ImageXSize = XSize;
00511   m_ImageYSize = YSize;
00512   CHECKNDELETE(m_UnscaledImage);
00513   m_UnscaledImage = new vuFixelMap<S,float>(m_ImageXSize, m_ImageYSize);
00514 
00515   setOversampling(overSampling);
00516 }
00517 
00518 // copy corresponding data to the zeroed wrap-arounds
00519 template <int S>
00520 void vuFourierVolume<S>::wrapVolume(float *volume)
00521 {
00522   // ?_?: to be done: depend wrapping on SIZE (number of channels)
00523   if (S==1) {
00524     dword index;
00525     dword ii, jj ,kk;
00526 
00527     dword hiX = m_XSize - m_Wrap - 1;
00528     dword hiY = m_YSize - m_Wrap - 1;
00529     dword hiZ = m_ZSize - m_Wrap - 1;
00530 
00531     dword XSize = m_XSize - m_Wrap * 2;
00532     dword YSize = m_YSize - m_Wrap * 2;
00533     dword ZSize = m_ZSize - m_Wrap * 2;
00534 
00535     index = 0;
00536     for(dword k = 0; k < m_ZSize; ++k) {
00537       for(dword j = 0; j < m_YSize; ++j) {
00538         for(dword i = 0; i < m_XSize; ++i) {
00539           if (i < m_Wrap)
00540             ii = i + XSize;
00541           else if (i > hiX)
00542             ii = i - XSize;
00543           else {
00544             index += 2;
00545             continue;
00546           }
00547 
00548           if (j < m_Wrap)
00549             jj = j + YSize;
00550           else if (j > hiY)
00551             jj = j - YSize;
00552           else {
00553             index += 2;
00554             continue;
00555           }
00556 
00557           if (k < m_Wrap)
00558             kk = k + ZSize;
00559           else if (k > hiZ)
00560             kk = k - ZSize;
00561           else {
00562             index += 2;
00563             continue;
00564           }
00565 
00566           ii = vcoord(ii, jj, kk);
00567 
00568           volume[index++] = volume[ii];
00569           volume[index++] = volume[ii+1];
00570         }
00571       }
00572     }
00573   }
00574 }
00575 
00576 /* ************************************************************************* */
00577 /* *** private fourier transformation functions **************************** */
00578 /* ************************************************************************* */
00579 
00580 template <int S>
00581 void vuFourierVolume<S>::loadWisdom()
00582 {
00583   cerr << "loading fftw wisdom..." << endl;
00584   FILE* wisdom_file = fopen(WISDOM, "a+");
00585   rewind(wisdom_file);
00586   fftw_import_wisdom_from_file(wisdom_file);
00587   fclose(wisdom_file);
00588 }
00589 
00590 template <int S>
00591 void vuFourierVolume<S>::saveWisdom()
00592 {
00593   cerr << "saving fftw wisdom..." << endl;
00594   FILE* wisdom_file = fopen(WISDOM, "w");
00595   fftw_export_wisdom_to_file(wisdom_file);
00596   fclose(wisdom_file);
00597 }
00598 
00599 template <int S>
00600 void vuFourierVolume<S>::transform2D(float *data)
00601 {
00602   fftwnd(m_Plan2D, S, (fftw_complex*)data, S, 1, NULL, 0, 0);
00603 }
00604 
00605 template <int S>
00606 void vuFourierVolume<S>::destroyTransform2D()
00607 {
00608   if (m_SliceXSize > 0 || m_SliceYSize > 0)
00609     fftwnd_destroy_plan(m_Plan2D);
00610   m_SliceXSize = 0;
00611   m_SliceYSize = 0;
00612 }
00613 
00614 template <int S>
00615 void vuFourierVolume<S>::shift2D(float* slice, dword XSize, dword YSize)
00616 {
00617   dword xstep = XSize * 2 * S;
00618   dword index;
00619   int start = 0;
00620   int val = -1;
00621 
00622   index = 0;
00623   for (dword j=0; j<YSize; j++) {
00624     for (dword i=start; i<xstep; i+=4*S) {
00625       float *ptr = &slice[index+i];
00626       for (int k=0; k<S; k++) {
00627         *(ptr++) *= -1.0f;
00628         *(ptr++) *= -1.0f;
00629       }
00630     }
00631     val   *= -1;
00632     start += val * 2 * S;
00633     index += xstep;
00634   }
00635 }
00636 
00637 // ?_?
00638 template <int S>
00639 void vuFourierVolume<S>::shiftCopy2D(float* dest, float* src, dword XSize,
00640                                      dword YSize)
00641 {
00642   dword i, j;
00643   dword xlast = XSize / 2;
00644   dword ylast = YSize / 2;
00645 
00646   float* src_ptr = src;
00647   float* dest_ptr = dest;
00648 
00649   for (j = 0; j < ylast; ++j) {
00650     for (i = 0; i < xlast; ++i) {
00651       *(dest_ptr++) = *(src_ptr++) * -1.0f;
00652       *(dest_ptr++) = *(src_ptr++) * -1.0f;
00653       *(dest_ptr++) = *(src_ptr++);
00654       *(dest_ptr++) = *(src_ptr++);
00655     }
00656     for (i = 0; i < xlast; ++i) {
00657       *(dest_ptr++) = *(src_ptr++);
00658       *(dest_ptr++) = *(src_ptr++);
00659       *(dest_ptr++) = *(src_ptr++) * -1.0f;
00660       *(dest_ptr++) = *(src_ptr++) * -1.0f;
00661     }
00662   }
00663 }
00664 
00665 /* ************************************************************************* */
00666 /* *** rendering *********************************************************** */
00667 /* ************************************************************************* */
00668 
00669 template <int S>
00670 void vuFourierVolume<S>::setRenderMethod(dword renderMethod)
00671 {
00672   if (renderMethod > 4) renderMethod = 0;
00673   m_RenderMethod = renderMethod;
00674 }
00675 template <int S>
00676 dword vuFourierVolume<S>::getRenderMethod()
00677 {
00678   return m_RenderMethod;
00679 }
00680 
00681 template <int S>
00682 void vuFourierVolume<S>::setScale(float scale)
00683 {
00684   m_Scale = scale;
00685 }
00686 template <int S>
00687 float vuFourierVolume<S>::getScale()
00688 {
00689   return m_Scale;
00690 }
00691 template <int S>
00692 void vuFourierVolume<S>::fitScale()
00693 {
00694   if (m_UnscaledImage == NULL) return;
00695 
00696   float minValue, maxValue;
00697 
00698   m_UnscaledImage->getMinAndMaxValue(minValue, maxValue);
00699 
00700   cerr << "fitScale=" << maxValue << endl;
00701 
00702   m_Scale = 1 / maxValue;
00703 }
00704 
00705 
00706 template <int S>
00707 void vuFourierVolume<S>::setBias(float bias)
00708 {
00709   m_Bias = bias;
00710 }
00711 template <int S>
00712 float vuFourierVolume<S>::getBias()
00713 {
00714   return m_Bias;
00715 }
00716 template <int S>
00717 void vuFourierVolume<S>::fitBias()
00718 {
00719   if (m_UnscaledImage == NULL) return;
00720 
00721   float minValue, maxValue;
00722 
00723   m_UnscaledImage->getMinAndMaxValue(minValue, maxValue);
00724 
00725   m_Bias = - minValue;
00726 }
00727 
00728 template <int S>
00729 bool vuFourierVolume<S>::getIsChannelActive(dword idx)
00730 {
00731   if (idx < S) 
00732     return m_Channels[idx];
00733   else
00734     return false;
00735 }
00736 
00737 template <int S>
00738 void vuFourierVolume<S>::setIsChannelActive(dword idx, bool flag)
00739 {
00740   if (idx < S)  m_Channels[idx] = flag;
00741 }
00742 
00743 template <int S>
00744 void vuFourierVolume<S>::setIsReRendering(bool flag)
00745 {
00746   m_IsReRendering = flag;
00747 }
00748 template <int S>
00749 bool vuFourierVolume<S>::getIsReRendering()
00750 {
00751   return m_IsReRendering;
00752 }
00753 
00754 template <int S>
00755 void vuFourierVolume<S>::glResize(dword width, dword height)
00756 {
00757   if (m_UnscaledImage == NULL) return;
00758   m_UnscaledImage->glResize(width, height);
00759 }
00760 
00761 template <int S>
00762 vuFixelMap<S,float> *vuFourierVolume<S>::getUnscaledImage()
00763 {
00764   return m_UnscaledImage;
00765 }
00766 
00767 
00768 template <int S>
00769 void vuFourierVolume<S>::computeUnscaledImage(bool doSlicing)
00770 {
00771   if (m_IsReRendering) {
00772     computeUnscaledImage(m_UnscaledImage, m_RenderMethod, doSlicing);
00773     m_IsReRendering = false;
00774   }
00775 }
00776 
00777 // chose the channels from m_Image, apply scale and bias  and display
00778 template <int S>
00779 void vuFourierVolume<S>::render()
00780 {
00781   computeUnscaledImage();
00782   // compute unscaled image (m_UnscaledImage)
00783   if (m_UnscaledImage == NULL) {
00784     cerr << "vuFourierVolume::displayImage: m_UnscaledImage is NULL" << endl;
00785     return;
00786   }
00787   vuFixelMap<S,float> *img = NULL;
00788 
00789   img = new vuFixelMap<S,float>(m_UnscaledImage->getWidth(),
00790                                 m_UnscaledImage->getHeight());
00791 
00792   bool status = true;
00793 
00794   for (int i=0; i<S; i++) status=status && m_Channels[i];
00795 
00796   if (status)
00797     *img = *m_UnscaledImage;
00798   else {
00799     vuFixelMap<1,float> *imag = NULL;
00800 
00801     imag = new vuFixelMap<1,float>(m_UnscaledImage->getWidth(),
00802                                    m_UnscaledImage->getHeight());
00803 
00804     for (int i=0; i<S; i++) {
00805       if (m_Channels[i]) {
00806         m_UnscaledImage->getChannel(imag,i);
00807         img->copyMapToChannel(imag,i);
00808       }
00809     }
00810     CHECKNDELETE(imag);
00811   }
00812 
00813   *img *= m_Scale;
00814   *img += m_Bias;
00815   img->glRender();
00816   CHECKNDELETE(img);
00817 }
00818 
00819 /* ************************************************************************* */
00820 /* *** advanced rendering ************************************************** */
00821 /* ************************************************************************* */
00822 
00823 /* method: 0 ... spatial image
00824            1 ... frequency amplitude
00825            2 ... frequency phase
00826            3 ... frequency real part
00827            4 ... frequency imaginary part
00828 */
00829 
00830 template <int S>
00831 void vuFourierVolume<S>::computeUnscaledImage(vuFixelMap<S,float>* &image,
00832                                               word method, bool doSlicing)
00833 {
00834   float minValue, maxValue;
00835   computeUnscaledImage(image, minValue, maxValue, method, doSlicing);
00836 }
00837 
00838 template <int S>
00839 void vuFourierVolume<S>::computeUnscaledImage(vuFixelMap<S,float>* &image,
00840                                               float &minVal,
00841                                               float &maxVal,
00842                                               word method, bool doSlicing)
00843 {
00844   if (!isInitialized()) {
00845     if (image == NULL) image = new vuFixelMap<S,float>(16,16);
00846     return;
00847   }
00848 
00849   if (method > 4) {
00850     cerr << "vuFourierVolume<S>::computeUnscaledImage():method not supported!";
00851     cerr << " Chose method=0 instead" << endl;
00852     method = 0;
00853   }
00854 
00855   dword width; 
00856   dword height;
00857   float *s_ptr;
00858 
00859   if (method == 0) { // spatial method
00860     if (doSlicing) computeSlice();
00861     width  = getImageWidth();
00862     height = getImageHeight();
00863     s_ptr  = m_SlicePtr;
00864   }
00865   else {             // frequency methods...
00866     if (doSlicing) interpolateSlice();
00867     width  = getSliceWidth();
00868     height = getSliceHeight();
00869     s_ptr  = m_Slice;
00870   }
00871 
00872   if (image == NULL)
00873     image = new vuFixelMap<S,float>(width, height);
00874   else if (width != image->getWidth() || height != image->getHeight()) {
00875     CHECKNDELETE(image);
00876     image = new vuFixelMap<S,float>(width, height);
00877   }
00878 
00879   _copyImageAndFindMinMax(image, s_ptr, minVal, maxVal, method);
00880 }
00881 
00882 template <int S>
00883 void vuFourierVolume<S>::_copyImageAndFindMinMax(vuFixelMap<S,float> *image,
00884                                                  float *slice,
00885                                                  float &minVal,
00886                                                  float &maxVal,
00887                                                  word  method)
00888 {
00889   dword width  = image->getWidth();
00890   dword height = image->getHeight();
00891 
00892   float *s_ptr = slice;
00893   float *i_ptr = image->getBuffer();
00894 
00895   minVal = + 10000000000.0f;
00896   maxVal = - 10000000000.0f;
00897 
00898   for (dword j=0; j<height; j++) {
00899     for (dword i=0; i<width; i++) {
00900       for (int k=0; k<S; k++) {
00901         register float value = *(s_ptr++);
00902 
00903         switch (method) {
00904         case 0: // image
00905           s_ptr++;
00906           break;
00907         case 1: // freq. amplitude
00908           value = sqrt(value*value + (*s_ptr) * (*s_ptr)); 
00909           s_ptr++;
00910           break;
00911         case 2: // freq. phase
00912           if ((*s_ptr) == 0)
00913             value = (value >= 0) ? 0.0f : M_PI;
00914           else
00915             value = atan(value/(*s_ptr));
00916           s_ptr++;
00917           break;
00918         case 3: // freq. real part
00919           s_ptr++;
00920           break;
00921         case 4: // freq. imaginary part
00922           value = *(s_ptr++);
00923           break;
00924         }
00925 
00926         *(i_ptr++) = value;
00927         if (minVal > value) minVal = value;
00928         if (maxVal < value) maxVal = value;
00929       }
00930       // ?_? ...
00931       if (method == 0) s_ptr += m_ImageStep;
00932     }
00933   }
00934 }
00935 
00936 
00937 
00938 
00939 
00940 
00941 // --- ??? kann weg ----------------------------------------------------------
00942 #if 1
00943 #include "vuFourier_IO.h"
00944 // ?_? todo
00945 // write the fourier volume to a file (WITHOUT m_Wrap)!!!
00946 template <int S>
00947 bool vuFourierVolume<S>::writeFourierToFile(const char *fileName,
00948                                             vuProgressHandler *handler) const
00949 {
00950   ofstream fout;
00951   fout.open(fileName, ios::out|ios::binary);
00952 
00953   dword XSize = m_XSize; // - (m_Wrap * 2);
00954   dword YSize = m_YSize; // - (m_Wrap * 2);
00955   dword ZSize = m_ZSize; // - (m_Wrap * 2);
00956 
00957   //??? ToDO: wrapping does not work!!!
00958   vuFourier_IO::write_fvr_head(fout, XSize, YSize, ZSize, sizeof(float));
00959   dword count = XSize*YSize*ZSize*2*sizeof(float);
00960 
00961   if (handler == NULL)
00962     fout.write((const char*)m_Volume, count);
00963   else {
00964 #define MEGABYTE 1000000
00965     int   cnt = count / MEGABYTE; // each 1MB one progress update
00966     char *ptr = (char *)m_Volume;
00967 
00968     int   oldProg  = handler->getCurrentProgress();
00969     float progress = (float)oldProg;
00970     float step     = (float)handler->getRange()/cnt;
00971 
00972     vuString str = "Write file: ";
00973     str += vuString(fileName).getLastPathComponent();
00974     handler->update(oldProg, str.c_str());
00975 
00976     for (int i=0; i<cnt; i++) {
00977       fout.write((const char*)ptr, MEGABYTE);
00978       ptr += MEGABYTE;
00979       progress += step;
00980 
00981       if (!handler->update((int)progress)) { // abort if user presses CANCEL
00982         fout.close();
00983         return false;
00984       }
00985     }
00986     int len = count % MEGABYTE;
00987     if (len) fout.write((const char*)ptr, len);
00988     handler->update(oldProg+handler->getRange());
00989   }
00990   fout.close();
00991 
00992   return true;
00993 }
00994 
00995 // ?_? todo
00996 // read the fourier volume from a file (data are stored without m_Wrap)!!!
00997 template <int S>
00998 bool vuFourierVolume<S>::readFourierFromFile(const char *fileName, dword wrap,
00999                                              float s)
01000 {
01001   dword XSize, YSize, ZSize, d_size;
01002 
01003   setWrap(wrap);
01004 
01005   ifstream fin(fileName, ios::in|ios::binary
01006 #ifdef IS_NOCREATE_NEEDED
01007 |ios::nocreate
01008 #endif
01009 );
01010   int ret = vuFourier_IO::read_head(fin, XSize, YSize, ZSize, d_size);
01011 
01012   if (ret != 2) {
01013     cerr << "readFourierFromFile(): wrong header type" << endl;
01014     return false;
01015   }
01016 
01017   if (d_size != sizeof(float)) {
01018     cerr << "readFourierFromFile(): wrong data type" << endl;
01019     return false;
01020   }
01021 
01022   m_XSize = XSize + m_Wrap * 2;
01023   m_YSize = YSize + m_Wrap * 2;
01024   m_ZSize = ZSize + m_Wrap * 2;
01025 
01026   m_Volume = new float[m_XSize * m_YSize * m_ZSize * 2];
01027 
01028   vuFourier_IO::read_raw_r_fast(fin, (byte*)m_Volume,
01029                                 m_XSize, m_YSize, m_ZSize,
01030                                 XSize, YSize, ZSize, d_size * 2);
01031   fin.close();
01032 
01033   wrapAndInitialize(s);
01034 
01035   return true;
01036 }
01037 
01038 // ?_? todo
01039 template <int S>
01040 bool vuFourierVolume<S>::writeSpatialVolume(const char *fileName,
01041                                             vuProgressHandler *handler) const
01042 {
01043   return true;
01044 }
01045 #endif

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