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
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
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
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
00089 {
00090 float *tmpVolume = m_Volume;
00091
00092
00093 m_Volume = m_CacheVolume;
00094
00095 addViewToVolume(view);
00096
00097 m_Volume = tmpVolume;
00098 }
00099
00100
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
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;
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
00204
00205
00206 template <int S, class T>
00207 void vuFourierCluster<S,T>::initializeVolume(dword width, dword height)
00208 {
00209
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);
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
00308 calcSliceDimensions(view, sliceWidth, sliceHeight);
00309
00310 float *slice = new float [sliceWidth * sliceHeight * 2 * S];
00311
00312 ensurePlan(sliceWidth, sliceHeight);
00313
00314
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
00332 handleSlice(slice, sliceWidth, sliceHeight, view);
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;
00389 mc_High = (width/2);
00390 mc_Low = -mc_High + (1- width % 2);
00391 mc_IsWidthOdd = (ODD(width));
00392
00393
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];
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 << " ";
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];
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 << " ";
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
00517
00518
00519
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();
00532 vuVector yStep = up.makeUnit();
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);
00555 mc_Low = -mc_High + (1- width % 2);
00556 mc_IsWidthOdd = (ODD(width));
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
00569 template <int S, class T>
00570 inline void vuFourierCluster<S,T>::doWeighting(vuVector& pos, float *value)
00571 {
00572 dword center[3];
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
00599
00600
00601
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
00622
00623 ptr[(j*width+i)*2] /= fac;
00624 }
00625 }
00626 }