00001 #include <dirent.h>
00002 #include <string.h>
00003 #include <sys/stat.h>
00004 #include <stdio.h>
00005 #include <iostream.h>
00006 #include "SimpleDefs.h"
00007 #include "FVR.h"
00008 #include "vuMatrix.h"
00009 #include "Image_io.h"
00010 #include "Transform.h"
00011 #include "SHarmonics.h"
00012
00013 #include <assert.h>
00014 #include <fftw.h>
00015 #include <unistd.h>
00016
00017 #include "vuTFunc/vuTFDesign.h"
00018
00019 using namespace FVR_NS;
00020
00021 #define ODD(x) ((x)&1)
00022
00023 #define TEMP_PATH "/tmp/vuVolume/"
00024
00025
00026
00027 float diff_intensity = .8;
00028 float cueCoeff = 15.0f;
00029
00030
00031 float my_dc[]={
00032 .1,
00033 3.14198,
00034 2.09440,2.09440,2.09440,
00035 0.78520,0.78520,0.78520,0.78520,0.78520};
00036
00037 float* vu1112117::getScale()
00038 {
00039 return &m_Scale;
00040 }
00041
00042 struct vuFVRarg {
00043 dword x_stop;
00044 dword y_stop;
00045 dword x_pass;
00046 dword y_pass;
00047 };
00048
00049
00050
00051
00052
00053 vu1112117::vu1112117() : m_currLight(1.0f, 0.0f, 0.0f),
00054 m_AmbientColour(1.0f, 1.0f, 1.0f),
00055 m_DiffuseColour(0.0f, 0.0f, 1.0f),
00056 m_XStep(1.0f, 0.0f, 0.0f),
00057 m_YStep(0.0f, 1.0f, 0.0f),
00058 m_XAxis(1.0f, 0.0f, 0.0f),
00059 m_YAxis(0.0f, 1.0f, 0.0f),
00060 m_ZAxis(0.0f, 0.0f, 1.0f)
00061 {
00062 m_volumeNumber = VOLUME_MAX_NUM;
00063 m_Scale = 1.5f;
00064 m_Bias = 0.0f;
00065 m_Wrap = 0;
00066 m_depthCue = true;
00067 m_shLight[0] = 1.0f;
00068
00069 for(int i = 0; i < VOLUME_MAX_NUM; i++) m_Yarray[i] = NULL;
00070 setLightPosition(vuVector(-1.0, 0.0, 0.0));
00071 initTransforms();
00072
00073
00074 for (int i=0; i<VOLUME_MAX_NUM; i++) m__SliceArray[i] = NULL;
00075 m_Slice = NULL;
00076
00077 vuTFDesign myTrans;
00078
00079 myTrans.addOpacity(1,0.01f);
00080 myTrans.addOpacity(255,1.0f);
00081 myTrans.setOpacitySmoothing(0);
00082 myTrans.setColourSmoothing(0);
00083
00084 myTrans.generateFunction();
00085 setTransferFunc(myTrans);
00086 }
00087
00088
00089
00090
00091
00092 #if 0
00093 vu1112117::vu1112117(const vu1112117& inst) : vu111211(inst)
00094 {
00095 assert("This has to be implemented: "
00096 "vu1112117::vu1112117(const vu1112117& inst)");
00097 }
00098 #endif
00099
00100
00101
00102
00103
00104 vu1112117::~vu1112117()
00105 {
00106 cerr << "destroying\n";
00107 destroyTransform2D();
00108 destroyTransforms();
00109
00110 deleteVolumes();
00111 destroySlices();
00112 if (m_Data) delete m_Data;
00113 }
00114
00115 void vu1112117::setIsDepthCueing(bool flag)
00116 {
00117 if (m_depthCue != flag) {
00118 m_depthCue = flag;
00119 }
00120 }
00121
00122 bool vu1112117::IsDepthCueing(void)
00123 {
00124 return m_depthCue;
00125 }
00126
00127 void vu1112117::deleteVolumes()
00128 {
00129 for(int i = 0; i < VOLUME_MAX_NUM; i++)
00130 if(m_Yarray[i] != NULL) {
00131 delete m_Yarray[i];
00132 m_Yarray[i] = NULL;
00133 }
00134 }
00135
00136 void vu1112117::setIsDiffuseShading(bool flag) {
00137 if (IsDiffuseShading() == flag)
00138 return;
00139 else if (flag)
00140 m_volumeNumber = VOLUME_MAX_NUM;
00141 else
00142 m_volumeNumber = 1;
00143 destroySlices();
00144 }
00145
00146 bool vu1112117::IsDiffuseShading() {
00147 if (this->m_volumeNumber == VOLUME_MAX_NUM)
00148 return true;
00149 else
00150 return false;
00151 }
00152
00153 void vu1112117::setAmbientColour(const vuColourRGBa _colour) {
00154 m_AmbientColour = _colour;
00155 }
00156 vuColourRGBa vu1112117::getAmbientColour(void) {
00157 return m_AmbientColour;
00158 }
00159
00160 void vu1112117::setDiffuseColour(const vuColourRGBa _colour) {
00161 m_DiffuseColour = _colour;
00162 }
00163 vuColourRGBa vu1112117::getDiffuseColour(void) {
00164 return m_DiffuseColour;
00165 }
00166
00167
00168
00169
00170
00171 void vu1112117::setViewVectors(const vuVector& view,const vuVector& up,const vuVector& right)
00172 {
00173 m_XAxis = right;
00174 m_YAxis = up;
00175 m_ZAxis = view;
00176
00177 m_XStep = right;
00178 m_YStep = up;
00179
00180 m_XStep.makeUnit();
00181 m_YStep.makeUnit();
00182
00183 m_XStep *= (float)(m_Image.getWidth() / m_SliceXSize);
00184 m_YStep *= (float)(m_Image.getHeight() / m_SliceYSize);
00185
00186 updateShLight();
00187 m_Origin = m_XAxis*(m_SliceXSize*-0.5f);
00188 m_Origin += m_YAxis*(m_SliceYSize*-0.5f);
00189 }
00190
00191
00192
00193
00194
00195 vu1112117& vu1112117::operator=(const vu1112117& rhs)
00196 {
00197 if (this != &rhs) {
00198 }
00199 return *this;
00200 }
00201
00202
00203
00204
00205
00206 bool vu1112117::doTempFilesExist(string fileName) {
00207 FILE *fileStream;
00208
00209 fileStream = fopen(fileName.c_str(), "r");
00210 if (fileStream != NULL) {
00211 fclose(fileStream);
00212 return true;
00213 }
00214 else
00215 return false;
00216 }
00217
00218 void vu1112117::ensureTempDirectory() {
00219 DIR *tmpDir;
00220
00221 tmpDir = opendir(TEMP_PATH);
00222 if (tmpDir == NULL) {
00223 #ifdef WIN32
00224 mkdir(TEMP_PATH);
00225 #else
00226 mkdir(TEMP_PATH,
00227 S_IRUSR | S_IWUSR | S_IXUSR |
00228 S_IRGRP | S_IWGRP | S_IXGRP |
00229 S_IROTH | S_IWOTH | S_IXOTH);
00230 #endif
00231 }
00232 else
00233 closedir(tmpDir);
00234 }
00235
00236 bool vu1112117::read()
00237 {
00238 string myFile(this->getFileName());
00239
00240 unsigned int len = myFile.length();
00241 unsigned int pos = myFile.rfind("/", len);
00242
00243 if ((pos != string::npos) && (pos < len)) {
00244 myFile = myFile.substr(pos+1, len);
00245 }
00246 len = myFile.length();
00247
00248 if (len >4) {
00249 myFile.replace(len-4, len, "0.fvr");
00250 myFile = string(TEMP_PATH) + myFile;
00251
00252 if (!this->doTempFilesExist(myFile)) {
00253 bool success = vu111211::read();
00254 if (!success) return false;
00255
00256 this->ensureTempDirectory();
00257
00258 preprocess();
00259
00260 write_fvr((char *)myFile.c_str());
00261 }
00262 {
00263 setWrap(m_Filter->getWidth()/ 2);
00264 read_shfvr((char *)myFile.c_str());
00265 _init();
00266 return true;
00267 }
00268 }
00269 return false;
00270 }
00271
00272 void vu1112117::_init(void) {
00273 for(int yl = 0; yl < VOLUME_MAX_NUM; yl++) {
00274 wrapVolume(m_Yarray[yl]);
00275 }
00276
00277 m_Image.init(m_Dim1Size, m_Dim2Size);
00278
00279 m_Origin[0] = (t_data)(m_Dim1Size / 2) * -1.0f;
00280 m_Origin[1] = (t_data)(m_Dim2Size / 2) * -1.0f;
00281 m_Origin[2] = 0.0f;
00282
00283 setOversampling(1.0);
00284
00285 initTransform2D(m_SliceXSize, m_SliceYSize);
00286 }
00287
00288
00289
00290
00291
00292 bool vu1112117::readRaw(void)
00293 {
00294 return false;
00295 }
00296
00297
00298
00299
00300
00301 void vu1112117::setIsPostProcessing(bool flag) {
00302 m_isPostProcessing = flag;
00303 }
00304 bool vu1112117::IsPostProcessing(void) {
00305 return m_isPostProcessing;
00306 }
00307
00308 void vu1112117::render(void) {
00309 getBuffer()->blit();
00310 }
00311
00312 vuImage* vu1112117::getBuffer ()
00313
00314 {
00315 if (IsReRendering()) {
00316 computeSlice();
00317 setIsReRendering(false);
00318 }
00319 if (IsPostProcessing()) {
00320 drawImageFromSlices();
00321 setIsPostProcessing(false);
00322 }
00323 return &m_Image;
00324 }
00325
00326
00327
00328
00329
00330 void vu1112117::initOpenGL(void)
00331 {
00332 glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
00333 }
00334
00335
00336
00337
00338
00339 void vu1112117::setLightPosition(const vuVector &_light)
00340 {
00341 t_data lv[3];
00342 lv[0] = _light[0]; lv[1] = _light[1]; lv[2] = _light[2];
00343
00344 t_data llv = sqrt(lv[0] * lv[0] + lv[1] * lv[1] + lv[2] * lv[2]);
00345 if(llv > 0) {
00346 lv[0] = lv[0] / llv;
00347 lv[1] = lv[1] / llv;
00348 lv[2] = lv[2] / llv;
00349 }
00350 else
00351 printf("Woooo %f\n", llv);
00352
00353 m_currLight[0] = lv[0];
00354 m_currLight[1] = lv[1];
00355 m_currLight[2] = lv[2];
00356
00357 updateShLight();
00358
00359 }
00360
00361 #if 0
00362 vuVector &vu1112117::getLightPosition(void) {
00363
00364 }
00365 #endif
00366
00367 void vu1112117::write_fvr(char* fileName)
00368 {
00369 cerr << "Writing transformed spherical harmonics basis to files...\n";
00370
00371 int len = strlen(fileName) + 1;
00372 char *out = new char[len];
00373 memset(out, 0, len);
00374 strcpy(out, fileName);
00375
00376 for(int yl = 0; yl < VOLUME_MAX_NUM; yl++) {
00377 ofstream fout;
00378 out[len - 6] = char('0' + yl);
00379 cerr << out << endl;
00380 fout.open(out, ios::out|ios::binary);
00381 write_fvr_head(fout, m_Dim1Size, m_Dim2Size, m_Dim3Size, sizeof(t_data));
00382 fout.write((char*)&m_Yarray[yl][0], m_Dim1Size * m_Dim2Size * m_Dim3Size * 2 * sizeof(t_data));
00383 fout.flush();
00384 fout.close();
00385 }
00386 delete out;
00387 }
00388
00389 int vu1112117::computeDimension(void)
00390 {
00391 dword dim = (m_Dim1Size > m_Dim2Size) ? m_Dim1Size : m_Dim2Size;
00392 dim = (dim > m_Dim3Size) ? dim : m_Dim3Size;
00393 dim = (ODD(dim)) ? dim + 1 : dim;
00394
00395 cerr << "dimension = " << dim << endl;
00396
00397 return dim;
00398 }
00399
00400 void vu1112117::preprocess(void)
00401 {
00402 int dim = computeDimension();
00403 int sizeX = m_Dim1Size;
00404 int sizeY = m_Dim2Size;
00405 int sizeZ = m_Dim3Size;
00406 t_data *myVol = new t_data[dim * dim * dim * 2];
00407
00408 m_Dim1Size = m_Dim2Size = m_Dim3Size = dim;
00409
00410
00411 read_raw(m_Data, myVol, m_Dim1Size, m_Dim2Size, m_Dim3Size,
00412 sizeX, sizeY, sizeZ, sizeof(byte));
00413
00414
00415 t_data *v_ptr = myVol;
00416 for (dword i = dim*dim*dim; i > 0; i--) {
00417 *v_ptr = (t_data)(m_TFunc.getOpacityAtPos((dword)(*v_ptr))*255);
00418 v_ptr+=2;
00419 }
00420 delete m_Data; m_Data = NULL;
00421
00422
00423 cerr << "Transforming to Real Spherical Harmonics...\n" << flush;
00424 realTrans(myVol, (int&)m_Dim1Size, (int&)m_Dim2Size, (int&)m_Dim3Size, m_Yarray);
00425 delete myVol;
00426 cerr << "Done.\n" << flush;
00427
00428 cerr << "Transforming the Spherical Harmonics basis to Frequency domain...";
00429 cerr << endl;
00430
00431 for(int yl = 0; yl < VOLUME_MAX_NUM; yl++) {
00432 if(yl == yl) {
00433 cerr << "shifting spherical basis " << yl << " ... " << flush;
00434 shift3D(m_Yarray[yl], m_Dim1Size, m_Dim2Size, m_Dim3Size);
00435 cerr << "done\ntransforming spherical basis " << yl << " ... " << flush;
00436 }
00437 initTransform3D(m_Dim1Size, m_Dim2Size, m_Dim3Size);
00438 transform3D(m_Yarray[yl]);
00439 destroyTransform3D();
00440
00441 if(yl == yl) {
00442 cerr << "done\nshifting spherical basis " << yl << " ... " << flush;
00443 shift3D(m_Yarray[yl], m_Dim1Size, m_Dim2Size, m_Dim3Size);
00444 cerr << "done\n" << flush;
00445 }
00446 }
00447 cerr << "Done.\n" << flush;
00448 }
00449
00450 bool vu1112117::read_fvr(t_data* &_volume,
00451 ifstream& fin,
00452 dword XSize,
00453 dword YSize,
00454 dword ZSize,
00455 dword d_size)
00456 {
00457 if (d_size != sizeof(t_data)) return false;
00458
00459 m_Dim1Size = XSize + m_Wrap * 2;
00460 m_Dim2Size = YSize + m_Wrap * 2;
00461 m_Dim3Size = ZSize + m_Wrap * 2;
00462
00463 _volume = new t_data[m_Dim1Size * m_Dim2Size * m_Dim3Size * 2];
00464
00465 read_raw_fast(fin, (byte*)_volume, m_Dim1Size, m_Dim2Size, m_Dim3Size, XSize, YSize, ZSize, d_size * 2);
00466
00467 return true;
00468 }
00469
00470 bool vu1112117::read_shfvr(char* fileName)
00471 {
00472 dword XSize, YSize, ZSize, d_size;
00473 int len = strlen(fileName);
00474 char *in = new char[len + 1];
00475 strcpy(in, fileName);
00476
00477 for(int yl = 0; yl < VOLUME_MAX_NUM; yl++) {
00478 in[len - 5] = '0' + yl;
00479 ifstream fin(in, ios::in|ios::binary
00480 #if IS_NOCREATE_NEEDED
00481 |ios::nocreate
00482 #endif
00483 );
00484
00485 int ret = read_head(fin, XSize, YSize, ZSize, d_size);
00486
00487 if(ret != 2)
00488 cerr << "Error reading basis " << yl << " from file " << in << endl;
00489 cerr << "reading spherical harmonics basis " << yl << " fvr file " << in << " ... " << flush;
00490
00491 if (!read_fvr(m_Yarray[yl], fin, XSize, YSize, ZSize, d_size))
00492 cerr << "Error reading basis " << yl << " from file " << in << endl;
00493 cerr << "done" << endl;
00494 fin.close();
00495 }
00496 delete in;
00497
00498 return true;
00499 }
00500
00501 void vu1112117::setOversampling(t_data s) {
00502
00503 if (s < 1.0) s = 1.0;
00504
00505 m_SliceXSize = (dword) ceil(m_Image.getWidth() * s);
00506 m_SliceYSize = (dword) ceil(m_Image.getHeight() * s);
00507 m_SliceXSize = (ODD(m_SliceXSize)) ? m_SliceXSize + 1 : m_SliceXSize;
00508 m_SliceYSize = (ODD(m_SliceYSize)) ? m_SliceYSize + 1 : m_SliceYSize;
00509
00510 m_XStep = m_XAxis * (float)(m_Image.getWidth() / m_SliceXSize);
00511 m_YStep = m_YAxis * (float)(m_Image.getHeight() / m_SliceYSize);
00512
00513 cerr << "Slice size: " << m_SliceXSize << ' ' << m_SliceYSize << endl;
00514
00515 m_InnerXSize = (dword) floor((t_data) m_SliceXSize / ROOT_TWO);
00516 m_InnerYSize = (dword) floor((t_data) m_SliceYSize / ROOT_TWO);
00517 m_InnerXSize = (ODD(m_InnerXSize)) ? m_InnerXSize + 1 : m_InnerXSize;
00518 m_InnerYSize = (ODD(m_InnerYSize)) ? m_InnerYSize + 1 : m_InnerYSize;
00519
00520 destroySlices();
00521 }
00522
00523 void vu1112117::setWrap(dword wrap)
00524 {
00525 m_Wrap = wrap;
00526 }
00527
00528 void vu1112117::setFilter(Filter* filter)
00529 {
00530 m_Filter = filter;
00531 }
00532
00533 void vu1112117::setSliceScale(t_data scale)
00534 {
00535 m_Scale = scale;
00536 }
00537
00538 void vu1112117::setSliceBias(t_data bias)
00539 {
00540 m_Bias = bias;
00541 }
00542
00543 void vu1112117::rotateSliceX(t_data alpha)
00544 {
00545 vuMatrix rot;
00546
00547 rot.makeRotate(m_XAxis, alpha);
00548
00549 m_YAxis = rot * m_YAxis;
00550 m_ZAxis = rot * m_ZAxis;
00551 m_XStep = rot * m_XStep;
00552 m_YStep = rot * m_YStep;
00553 m_Origin = rot * m_Origin;
00554
00555 updateShLight();
00556 }
00557
00558 void vu1112117::rotateSliceY(t_data alpha)
00559 {
00560 vuMatrix rot;
00561
00562 rot.makeRotate(m_YAxis, alpha);
00563
00564 m_XAxis = rot * m_XAxis;
00565 m_ZAxis = rot * m_ZAxis;
00566 m_XStep = rot * m_XStep;
00567 m_YStep = rot * m_YStep;
00568 m_Origin = rot * m_Origin;
00569
00570 updateShLight();
00571 }
00572
00573 void vu1112117::rotateSliceZ(t_data alpha)
00574 {
00575 vuMatrix rot;
00576
00577 rot.makeRotate(m_ZAxis, alpha);
00578
00579 m_XAxis = rot * m_XAxis;
00580 m_YAxis = rot * m_YAxis;
00581 m_XStep = rot * m_XStep;
00582 m_YStep = rot * m_YStep;
00583 m_Origin = rot * m_Origin;
00584 updateShLight();
00585 }
00586
00587 dword vu1112117::getImageWidth(void) const
00588 {
00589 return m_Image.getWidth();
00590 }
00591
00592 dword vu1112117::getImageHeight(void) const
00593 {
00594 return m_Image.getHeight();
00595 }
00596
00597 dword vu1112117::getSliceWidth(void) const
00598 {
00599 return m_SliceXSize;
00600 }
00601
00602 dword vu1112117::getSliceHeight(void) const
00603 {
00604 return m_SliceYSize;
00605 }
00606
00607
00608 void vu1112117::wrapVolume(t_data* &_volume)
00609 {
00610 dword index;
00611 dword ii, jj ,kk;
00612
00613 dword hiX = m_Dim1Size - m_Wrap - 1;
00614 dword hiY = m_Dim2Size - m_Wrap - 1;
00615 dword hiZ = m_Dim3Size - m_Wrap - 1;
00616
00617 dword XSize = m_Dim1Size - m_Wrap * 2;
00618 dword YSize = m_Dim2Size - m_Wrap * 2;
00619 dword ZSize = m_Dim3Size - m_Wrap * 2;
00620
00621 index = 0;
00622 for(dword k = 0; k < m_Dim3Size; ++k) {
00623 for(dword j = 0; j < m_Dim2Size; ++j) {
00624 for(dword i = 0; i < m_Dim1Size; ++i) {
00625 if (i < m_Wrap)
00626 ii = i + XSize;
00627 else if (i > hiX)
00628 ii = i - XSize;
00629 else {
00630 index += 2;
00631 continue;
00632 }
00633
00634 if (j < m_Wrap)
00635 jj = j + YSize;
00636 else if (j > hiY)
00637 jj = j - YSize;
00638 else {
00639 index += 2;
00640 continue;
00641 }
00642
00643 if (k < m_Wrap)
00644 kk = k + ZSize;
00645 else if (k > hiZ)
00646 kk = k - ZSize;
00647 else {
00648 index += 2;
00649 continue;
00650 }
00651
00652 ii = vcoord(ii, jj, kk);
00653
00654 _volume[index++] = _volume[ii];
00655 _volume[index++] = _volume[ii+1];
00656 }
00657 }
00658 }
00659 }
00660
00661 void vu1112117::destroySlices(void)
00662 {
00663 for (int i=0; i< VOLUME_MAX_NUM; i++) {
00664 if (m__SliceArray[i]) {
00665 delete [] m__SliceArray[i];
00666 m__SliceArray[i] = NULL;
00667 }
00668 }
00669 if (m_Slice != NULL) {
00670 delete [] m_Slice;
00671 m_Slice = NULL;
00672 }
00673 }
00674
00675 void vu1112117::ensureSlices(void)
00676 {
00677 if (m__SliceArray[0] == NULL) {
00678 dword length = m_SliceXSize * m_SliceYSize * 2;
00679 for (int i=0; i<m_volumeNumber; i++) m__SliceArray[i] = new t_data[length];
00680 _clearSlices();
00681 }
00682 }
00683
00684 void vu1112117::clearSlices(void)
00685 {
00686 ensureSlices();
00687 _clearSlices();
00688 }
00689
00690 void vu1112117::_clearSlices(void)
00691 {
00692 dword cnt = m_SliceXSize * m_SliceYSize * 2;
00693 for (int num = 0; num < m_volumeNumber; num++) {
00694 t_data* s_ptr = m__SliceArray[num];
00695 for (dword i = 0; i < cnt; i++) *(s_ptr++) = 0.0f;
00696 }
00697 }
00698
00699 void vu1112117::_ensureSliceDummy(void) {
00700 dword length = m_SliceXSize * m_SliceYSize * 2;
00701 if (m_Slice == NULL)
00702 m_Slice = new t_data[length];
00703
00704 t_data* s_ptr = m_Slice;
00705 for (dword i = 0; i < length; i++) *(s_ptr++) = 0.0f;
00706 }
00707
00708
00709
00710 #define CONVOLVE_NO_MOD(_volume_) \
00711 { \
00712 int x = (int) p[0]; \
00713 int y = (int) p[1]; \
00714 int z = (int) p[2]; \
00715 \
00716 t[0] = p[0] - x; \
00717 t[1] = p[1] - y; \
00718 t[2] = p[2] - z; \
00719 \
00720 x = x - XSize + m_Wrap; \
00721 y = y - YSize + m_Wrap; \
00722 z = z - ZSize + m_Wrap; \
00723 \
00724 register t_data* fweight = m_Filter->getWeights(t); \
00725 register vuVector* fgrad = m_Filter->getGrad(t); \
00726 \
00727 register t_data real = 0.0f; \
00728 register t_data imag = 0.0f; \
00729 \
00730 register t_data* v_ptr = &_volume_[vcoord(x + f2, y + f2, z + f2)]; \
00731 \
00732 for (dword l = 0; l < f_len; l++) { \
00733 for (dword m = 0; m < f_len; m++) { \
00734 for (dword n = 0; n < f_len; n++) { \
00735 t_data tmp = *(fweight++); \
00736 t_data tmpi = 0.0; \
00737 \
00738 if(m_depthCue)tmpi = (fgrad++)->dot(m_ZAxis) * cueCoeff/(2.0*M_PI); \
00739 real += (*(v_ptr) * tmp - *(v_ptr + 1) * tmpi); \
00740 imag += (*(v_ptr + 1) * tmp + *(v_ptr) * tmpi); \
00741 v_ptr -= 2; \
00742 } \
00743 v_ptr += fXStep; \
00744 } \
00745 v_ptr += fYStep; \
00746 } \
00747 *(s_ptr++) += light_brdf * real; \
00748 *(s_ptr++) += light_brdf * imag; \
00749 p += m_XStep; \
00750 }
00751
00752 #define CONVOLVE(_volume_) \
00753 { \
00754 int x = (int) p[0]; \
00755 int y = (int) p[1]; \
00756 int z = (int) p[2]; \
00757 \
00758 t[0] = p[0] - x; \
00759 t[1] = p[1] - y; \
00760 t[2] = p[2] - z; \
00761 \
00762 x = x % XSize + m_Wrap; \
00763 y = y % YSize + m_Wrap; \
00764 z = z % ZSize + m_Wrap; \
00765 \
00766 register t_data* fweight = m_Filter->getWeights(t); \
00767 register vuVector* fgrad = m_Filter->getGrad(t); \
00768 \
00769 register t_data real = 0.0; \
00770 register t_data imag = 0.0; \
00771 \
00772 register t_data* v_ptr = &_volume_[vcoord(x + f2, y + f2, z + f2)]; \
00773 \
00774 for (dword l = 0; l < f_len; l++) { \
00775 for (dword m = 0; m < f_len; m++) { \
00776 for (dword n = 0; n < f_len; n++) { \
00777 t_data tmp = *(fweight++); \
00778 t_data tmpi = 0.0; \
00779 \
00780 if(m_depthCue)tmpi = (fgrad++)->dot(m_ZAxis) * cueCoeff/(2.0*M_PI); \
00781 real += (*(v_ptr) * tmp - *(v_ptr + 1) * tmpi); \
00782 imag += (*(v_ptr + 1) * tmp + *(v_ptr) * tmpi); \
00783 v_ptr -= 2; \
00784 } \
00785 v_ptr += fXStep; \
00786 } \
00787 v_ptr += fYStep; \
00788 } \
00789 *(s_ptr++) += light_brdf * real; \
00790 *(s_ptr++) += light_brdf * imag; \
00791 p += m_XStep; \
00792 }
00793
00794 void vu1112117::interpolateSlice(t_data* &_volume, t_data* &_slice,
00795 t_data &_shLight, t_data &_currDC,
00796 dword x_stop, dword y_stop,
00797 dword x_pass, dword y_pass)
00798 {
00799 vuVector t;
00800 dword f_len = m_Filter->getWidth();
00801 int f2 = f_len / 2;
00802 x_stop *= 2;
00803 y_stop *= 2;
00804 dword stop_l = (m_SliceXSize - x_stop - 2 * x_pass) / 2;
00805 dword stop_r = m_SliceXSize - x_stop - 2 * x_pass - stop_l;
00806 dword stop_b = (m_SliceYSize - y_stop - 2 * y_pass) / 2;
00807 dword x_hi = x_stop + 2 * x_pass;
00808 dword s_step = stop_l + stop_r;
00809
00810 dword XSize = m_Dim1Size - m_Wrap * 2;
00811 dword YSize = m_Dim2Size - m_Wrap * 2;
00812 dword ZSize = m_Dim3Size - m_Wrap * 2;
00813
00814 dword fXStep = (f_len - m_Dim1Size) * 2;
00815 dword fYStep = (f_len - m_Dim2Size) * m_Dim1Size * 2;
00816
00817 vuVector q = m_Origin;
00818 q[0] += XSize / 2 + XSize + stop_l * m_XStep[0] + stop_b * m_YStep[0];
00819 q[1] += YSize / 2 + YSize + stop_l * m_XStep[1] + stop_b * m_YStep[1];
00820 q[2] += ZSize / 2 + ZSize + stop_l * m_XStep[2] + stop_b * m_YStep[2];
00821
00822 t_data* s_ptr = &_slice[scoord(stop_l, stop_b)];
00823
00824 t_data light_brdf = _shLight * _currDC;
00825
00826 for (dword j = 0; j < y_pass; ++j) {
00827 vuVector p = q;
00828 for (dword i = 0; i < x_hi; ++i) CONVOLVE(_volume);
00829 s_ptr += s_step * 2;
00830 q += m_YStep;
00831 }
00832
00833 for (dword j = 0; j < y_stop; ++j) {
00834 vuVector p = q;
00835 for (dword i = 0; i < x_pass; ++i) CONVOLVE(_volume);
00836
00837 s_ptr += x_stop * 2;
00838 p[0] += x_stop * m_XStep[0];
00839 p[1] += x_stop * m_XStep[1];
00840 p[2] += x_stop * m_XStep[2];
00841
00842 for (dword i = 0; i < x_pass; ++i) CONVOLVE(_volume);
00843 s_ptr += s_step * 2;
00844 q += m_YStep;
00845 }
00846
00847 for (dword j = 0; j < y_pass; ++j) {
00848 vuVector p = q;
00849 for (dword i = 0; i < x_hi; ++i) CONVOLVE(_volume);
00850 s_ptr += s_step * 2;
00851 q += m_YStep;
00852 }
00853 }
00854
00855 void vu1112117::interpolateSlice(t_data* &_volume, t_data* &_slice,
00856 t_data &_shLight, t_data &_currDC)
00857 {
00858 dword i, j;
00859 vuVector t;
00860 dword f_len = m_Filter->getWidth();
00861 int f2 = f_len / 2;
00862
00863 dword x_pass = (m_SliceXSize - m_InnerXSize) / 2;
00864 dword y_pass = (m_SliceYSize - m_InnerYSize) / 2;
00865
00866 dword XSize = m_Dim1Size - m_Wrap * 2;
00867 dword YSize = m_Dim2Size - m_Wrap * 2;
00868 dword ZSize = m_Dim3Size - m_Wrap * 2;
00869
00870 dword fXStep = (f_len - m_Dim1Size) * 2;
00871 dword fYStep = (f_len - m_Dim2Size) * m_Dim1Size * 2;
00872
00873 vuVector q = m_Origin;
00874 q[0] += XSize / 2 + XSize;
00875 q[1] += YSize / 2 + YSize;
00876 q[2] += ZSize / 2 + ZSize;
00877
00878 t_data* s_ptr = _slice;
00879
00880 t_data light_brdf = _shLight * _currDC;
00881
00882 for (j = 0; j < y_pass; ++j) {
00883 vuVector p = q;
00884 for (i = 0; i < m_SliceXSize; ++i) CONVOLVE(_volume);
00885 q += m_YStep;
00886 }
00887
00888 for (j = 0; j < m_InnerYSize; ++j) {
00889 vuVector p = q;
00890 for (i = 0; i < x_pass; ++i) CONVOLVE(_volume);
00891 for (i = 0; i < m_InnerXSize; ++i) CONVOLVE_NO_MOD(_volume);
00892 for (i = 0; i < x_pass; ++i)CONVOLVE(_volume);
00893 q += m_YStep;
00894 }
00895
00896 for (j = 0; j < y_pass; ++j) {
00897 vuVector p = q;
00898 for (i = 0; i < m_SliceXSize; ++i) CONVOLVE(_volume);
00899 q += m_YStep;
00900 }
00901 }
00902
00903 void vu1112117::accumulateSlices(t_data* _slice, int _from, int _to) {
00904 dword cnt = m_SliceXSize * m_SliceYSize * 2;
00905 for (int num = _from; num < _to; num++) {
00906 t_data* src_ptr = m__SliceArray[num];
00907 t_data* dest_ptr = _slice;
00908 for (dword i = 0; i < cnt; i++) *(dest_ptr++) += *(src_ptr++);
00909 }
00910 }
00911
00912 void vu1112117::run(int whatsup, void* data) {
00913 vuFVRarg *arg = (vuFVRarg*)data;
00914 if (data == NULL)
00915 computeSlice(whatsup);
00916 else
00917 refineSlice(whatsup, arg->x_stop, arg->y_stop, arg->x_pass, arg->y_pass);
00918 m_Mutex[whatsup].unlock();
00919 }
00920
00921 void vu1112117::computeSlice(int num)
00922 {
00923 t_data currDC = (num==0?1.0:diff_intensity) * my_dc[num];
00924 interpolateSlice(m_Yarray[num], m__SliceArray[num], m_shLight[num], currDC);
00925 }
00926
00927 void vu1112117::refineSlice(int num,
00928 dword x_stop, dword y_stop,
00929 dword x_pass, dword y_pass)
00930 {
00931 t_data currDC = (num==0?1.0:diff_intensity) * my_dc[num];
00932 interpolateSlice(m_Yarray[num], m__SliceArray[num],
00933 m_shLight[num], currDC,
00934 x_stop, y_stop, x_pass, y_pass);
00935 }
00936
00937 bool vu1112117::computeSlicesUsingThreads(void) {
00938 #ifdef NO_THREADS
00939 return false;
00940 #endif
00941 for (int i=0; i < m_volumeNumber; i++) m_Mutex[i].lock();
00942 if (startThread(0)) {
00943 for (int i=1; i < m_volumeNumber; i++) {
00944 startThread(i);
00945 }
00946 }
00947 else {
00948 return false;
00949 }
00950 for (int i=0; i<m_volumeNumber; i++) m_Mutex[i].lock();
00951 for (int i=0; i<m_volumeNumber; i++) m_Mutex[i].unlock();
00952 return true;
00953 }
00954
00955 bool vu1112117::refineSlicesUsingThreads(dword x_stop, dword y_stop,
00956 dword x_pass, dword y_pass)
00957 {
00958 vuFVRarg arg;
00959 arg.x_stop = x_stop; arg.y_stop = y_stop;
00960 arg.x_pass = x_pass; arg.y_pass = y_pass;
00961
00962 #ifdef NO_THREADS
00963 return false;
00964 #endif
00965 for (int i=0; i < m_volumeNumber; i++) m_Mutex[i].lock();
00966 if (startThread(0, &arg)) {
00967 for (int i=1; i < m_volumeNumber; i++) {
00968 startThread(i, &arg);
00969 }
00970 }
00971 else {
00972 return false;
00973 }
00974 for (int i=0; i<m_volumeNumber; i++) m_Mutex[i].lock();
00975 for (int i=0; i<m_volumeNumber; i++) m_Mutex[i].unlock();
00976 return true;
00977 }
00978
00979 void vu1112117::drawImageFromSlices(void) {
00980 t_data *slice = NULL;
00981
00982 ensureSlices();
00983
00984 m_Image.clear();
00985 #if 1
00986
00987 _ensureSliceDummy();
00988 slice = m_Slice;
00989 accumulateSlices(slice, 0, 1);
00990 shift2D(slice, m_SliceXSize, m_SliceYSize);
00991 transform2D(slice);
00992 shift2D(slice, m_SliceXSize, m_SliceYSize);
00993 scaleAndBias(slice, m_AmbientColour);
00994 #endif
00995
00996 _ensureSliceDummy();
00997 if (IsDiffuseShading()) {
00998
00999 slice = m_Slice;
01000 accumulateSlices(slice, 1, m_volumeNumber);
01001 shift2D(slice, m_SliceXSize, m_SliceYSize);
01002 transform2D(slice);
01003 shift2D(slice, m_SliceXSize, m_SliceYSize);
01004 scaleAndBias(slice, m_DiffuseColour);
01005 }
01006
01007 }
01008
01009 void vu1112117::computeSlice(void)
01010 {
01011 clearSlices();
01012
01013 if (!computeSlicesUsingThreads()) {
01014 for(int yl = 0; yl < this->m_volumeNumber; yl++) computeSlice(yl);
01015 }
01016 setIsPostProcessing(true);
01017 }
01018
01019 void vu1112117::refineSlice(dword x_stop, dword y_stop,
01020 dword x_pass, dword y_pass)
01021 {
01022 ensureSlices();
01023
01024 if (!refineSlicesUsingThreads(x_stop, y_stop, x_pass, y_pass)) {
01025 for(int yl = 0; yl < this->m_volumeNumber; yl++) {
01026 refineSlice(yl, x_stop, y_stop, x_pass, y_pass);
01027 }
01028 }
01029 drawImageFromSlices();
01030 }
01031
01032 void vu1112117::scaleAndBias(t_data* _slice, const vuColourRGBa &color) {
01033 #define clip(a) (a<0?0:a>1?1:a)
01034 #define clip255(a) (a>255?255:a)
01035 dword imgWidth = m_Image.getWidth();
01036 dword imgHeight = m_Image.getHeight();
01037 dword imgStep = (m_SliceXSize - imgWidth) * 2;
01038 byte* i_ptr = m_Image.get_buffer();
01039 t_data* s_ptr = _slice;
01040
01041 for (dword j = 0; j < imgHeight; j++) {
01042 for (dword i = 0; i < imgWidth; i++) {
01043 register dword tmp;
01044
01045 tmp = (dword)(clip(*s_ptr * m_Scale + m_Bias) * color[0] * 255);
01046 *(i_ptr++) = (byte)clip255(*i_ptr+tmp);
01047 tmp = (dword)(clip(*s_ptr * m_Scale + m_Bias) * color[1] * 255);
01048 *(i_ptr++) = (byte)clip255(*i_ptr+tmp);
01049 tmp = (dword)(clip(*s_ptr * m_Scale + m_Bias) * color[2] * 255);
01050 *(i_ptr++) = (byte)clip255(*i_ptr+tmp);
01051 s_ptr += 2;
01052 }
01053 s_ptr += imgStep;
01054 }
01055 }
01056
01057 void vu1112117::updateShLight()
01058 {
01059 vuVector newLight(0,0,0);
01060 vuVector vx(0,0,0), vy(0,0,0), vz(0,0,0);
01061 vx = m_currLight[0] * m_XAxis;
01062 vy = m_currLight[1] * m_YAxis;
01063 vz = m_currLight[2] * m_ZAxis;
01064 newLight = vx + vy + vz;
01065 realLight(&m_shLight[1], newLight.getData());
01066 }
01067
01068
01069
01070
01071
01072
01073 #if 0
01074 int indexOf(int x, int xMax, int y, int yMax)
01075 {
01076 return 2 * (y * xMax + x);
01077 }
01078
01079 void refine(t_data *fSlice, int M, int L)
01080 #define clip(a) (a<0 ? 0 : a)
01081 {
01082 t_data phaser;
01083 t_data maxval = -HUGE;
01084 int l, m;
01085 for(l=0;l<L;l++) {
01086 for(m=0;m<M;m++) {
01087 phaser=(l+m) % 2 ==0 ? 1:-1;
01088 if(clip(fSlice[indexOf(m, M, l, L)]) > maxval)
01089 maxval = clip(fSlice[indexOf(m, M, l, L)]);
01090 }
01091 }
01092 cerr << "maxval is " << maxval << endl;
01093 if(maxval >0) maxval = 1.0/maxval;
01094 for(l=0;l<L;l++) {
01095 for(m=0;m<M;m++) {
01096 fSlice[indexOf(m, M, l, L)] = (clip(fSlice[indexOf(m, M, l, L)])*maxval);
01097 }
01098 }
01099 }
01100 #endif