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

FVR.cpp

Go to the documentation of this file.
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 // #define NO_THREADS
00026 
00027 float diff_intensity = .8;
00028 float cueCoeff = 15.0f;//.9;
00029 // the harmonic transform coefficients for
00030 // diffuse BRDF with max function
00031 float my_dc[]={
00032   .1, // ambient coefficient
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 //------------------------- The default constructor --------------------------
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   // slice
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 //------------------------- The copy constructor -----------------------------
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 //------------------------- The destructor -----------------------------------
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 //------------------------- public setViewVectors() --------------------------
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 //------------------------- The assignment operator --------------------------
00193 //----------------------------------------------------------------------------
00194 
00195 vu1112117& vu1112117::operator=(const vu1112117& rhs)
00196 {
00197   if (this != &rhs) {
00198   }
00199   return *this;
00200 }
00201 
00202 //----------------------------------------------------------------------------
00203 //------------------------- public read() ------------------------------------
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) { // myFile must contain at least the file extension
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     { // read preprocessed files
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 //------------------------- public readRaw() ---------------------------------
00290 //----------------------------------------------------------------------------
00291 
00292 bool vu1112117::readRaw(void)
00293 {
00294   return false;
00295 }
00296 
00297 //----------------------------------------------------------------------------
00298 //------------------------- public render() ----------------------------------
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 //------------------------- public initOpenGL() ------------------------------
00328 //----------------------------------------------------------------------------
00329 
00330 void vu1112117::initOpenGL(void)
00331 {
00332   glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
00333 }
00334 
00335 //----------------------------------------------------------------------------
00336 //---------------------------- public methods --------------------------------
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   //realLight(&m_shLight[1], lv);
00359 }
00360 
00361 #if 0
00362 vuVector &vu1112117::getLightPosition(void) {
00363   
00364 }
00365 #endif
00366 
00367 void vu1112117::write_fvr(char* fileName)// const
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   // copy m_Data to myVol and make dimensions equal
00411   read_raw(m_Data, myVol, m_Dim1Size, m_Dim2Size, m_Dim3Size,
00412            sizeX, sizeY, sizeZ, sizeof(byte));
00413 
00414   // applying transfer function
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) { //0 
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   // don't allow undersampling
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 // fill in zeroed wrap-arounds
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 // and may the gods of good-style forgive me ...
00709 // inner convolution loop, should be fast!
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         /* vuVector *vFGrad = (fgrad++);  */                                \
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         /* vuVector *vFGrad = (fgrad++);*/                                  \
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); // unlocks m_Mutex[i] if ready
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); // unlocks m_Mutex[i] if ready
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   // compute ambient light
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     // compute diffuse light
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  ***                            not used                                  ***
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) // remove neg values due to approx in image
01081 {
01082   t_data phaser; // undo the shift
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

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