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;
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;
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
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
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
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
00519 template <int S>
00520 void vuFourierVolume<S>::wrapVolume(float *volume)
00521 {
00522
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
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
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
00778 template <int S>
00779 void vuFourierVolume<S>::render()
00780 {
00781 computeUnscaledImage();
00782
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
00821
00822
00823
00824
00825
00826
00827
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) {
00860 if (doSlicing) computeSlice();
00861 width = getImageWidth();
00862 height = getImageHeight();
00863 s_ptr = m_SlicePtr;
00864 }
00865 else {
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:
00905 s_ptr++;
00906 break;
00907 case 1:
00908 value = sqrt(value*value + (*s_ptr) * (*s_ptr));
00909 s_ptr++;
00910 break;
00911 case 2:
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:
00919 s_ptr++;
00920 break;
00921 case 4:
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
00942 #if 1
00943 #include "vuFourier_IO.h"
00944
00945
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;
00954 dword YSize = m_YSize;
00955 dword ZSize = m_ZSize;
00956
00957
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;
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)) {
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
00996
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
01039 template <int S>
01040 bool vuFourierVolume<S>::writeSpatialVolume(const char *fileName,
01041 vuProgressHandler *handler) const
01042 {
01043 return true;
01044 }
01045 #endif