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

Regular/Unimodal/3d/1B/Intensity/ShearWarp/shearWarp.cpp

Go to the documentation of this file.
00001 
00002 //                                                                    //
00003 // Informatikpraktikum I "Implementation of the ShearWarp Algorithm"  //
00004 // Februar - October 2002                                             //
00005 //                                                                    //
00006 // author: Sebastian Zambal                                           //
00007 //         e9826978@student.tuwien.ac.at                              //
00008 //                                                                    //
00009 // file:   shearWarp.cpp (regular)                                    //
00010 //                                                                    //
00012 
00013 #include "shearWarp.h"
00014 #include "vuGui/vuDrawTools.h"
00015 #include "../simpleFuncs.h"
00016 #include <fstream.h>
00017 #include <math.h>
00018 
00019 //-----------------------------------------------------------------------
00020 //-- Constructor                                                      ---
00021 //-----------------------------------------------------------------------
00022 vu111211A::vu111211A()
00023 {
00024     //initially no zoom
00025     zoomValue = 1;
00026 
00027     //initially don't use OpenGL for warping
00028     orthoWarpOpenGL = 0;
00029 
00030     //initially no normals
00031     m_Normals = 0;
00032 
00033     //initially there is no runlength encoded data
00034     dataX = 0;
00035     dataY = 0;
00036     dataZ = 0;
00037 
00038     //initially no current runlength encoding
00039     curr = 0;
00040     //initially no intermediate image
00041     intermediate = 0;
00042     //initially orthogonal projection
00043     viewingMode = VIEWING_MODE_ORTHO;
00044     //initially no specular shading
00045     specular = 0;
00046     //initially no fast classification
00047     fastClassification = 0;
00048 
00049     m_GLShearWarp = 0;
00050 }
00051 
00052 //-----------------------------------------------------------------------
00053 //-- Destructor                                                       ---
00054 //-----------------------------------------------------------------------
00055 vu111211A::~vu111211A() 
00056 {
00057     if (glImage) {
00058         delete [] glImage;
00059         glImage = 0;
00060     }
00061     if (m_Normals) delete [] m_Normals;
00062     if (dataX) {
00063         for (dword i = 0; i < m_Dim1Size; ++i) { 
00064             if (dataX[i].size>0) {
00065                 delete [] dataX[i].voxel;
00066                 delete [] dataX[i].runlength;
00067             }
00068         }
00069         delete [] dataX;
00070     }
00071     if (dataY) {
00072         for (dword i = 0; i < m_Dim2Size; ++i) { 
00073             if (dataY[i].size>0) {
00074                 delete [] dataY[i].voxel;
00075                 delete [] dataY[i].runlength;
00076             }
00077         }
00078         delete [] dataY;
00079     }
00080     if (dataZ) {
00081         for (dword i = 0; i < m_Dim3Size; ++i) { 
00082             if (dataZ[i].size>0) {
00083                 delete [] dataZ[i].voxel;
00084                 delete [] dataZ[i].runlength;
00085             }
00086         }
00087         delete [] dataZ;
00088     }
00089     if (intermediate) {
00090         delete [] intermediate;
00091     }
00092     if (fastClassificationObj) {
00093         delete fastClassificationObj;
00094     }
00095 }
00096 
00097 //-----------------------------------------------------------------------
00098 //-- returns the size of dimensions of the volume data                ---
00099 //-----------------------------------------------------------------------
00100 void vu111211A::getDimensions(int &x, int &y, int &z) {
00101     x = m_Dim1Size;
00102     y = m_Dim2Size;
00103     z = m_Dim3Size;
00104 }
00105 
00106 //-----------------------------------------------------------------------
00107 //-- Construct OpenGL image                                           ---
00108 //-----------------------------------------------------------------------
00109 void vu111211A::createGLImage(void) {
00110 
00111     word glSize = 64;
00112     while (glSize < maxSize*2) {
00113         glSize *= 2;
00114     }
00115 //    glSize *= 2;
00116 
00117     glImageWidth = glSize;
00118     glImageHeight = glSize;
00119     glImage = new GLubyte[glImageWidth*glImageHeight*4];
00120 }
00121 
00122 //-----------------------------------------------------------------------
00123 //-- Compute maximum stretch-out of volume                            ---
00124 //-----------------------------------------------------------------------
00125 void vu111211A::computeMaxSize(void) 
00126 {
00127     maxSize = m_Dim1Size;
00128     if (m_Dim2Size > maxSize) maxSize = m_Dim2Size;
00129     if (m_Dim3Size > maxSize) maxSize = m_Dim3Size;
00130     maxSize += 2;                         
00131 }
00132 
00133 //-----------------------------------------------------------------------
00134 //-- Set view vectors                                                 ---
00135 //-----------------------------------------------------------------------
00136 void vu111211A::setViewVectors(const vuVector& view, const vuVector& up, const vuVector& right)
00137 {
00138   vuVector view1;
00139   vuVector up1;
00140   vuVector right1;
00141 
00142   m_View = view;
00143 
00144   m_Up = up;
00145   m_Right = right;
00146 
00147 }
00148 
00149 //-----------------------------------------------------------------------
00150 //-- Set viewing Mode                                                 ---
00151 //-----------------------------------------------------------------------
00152 void vu111211A::setViewing(int vm)
00153 {
00154     viewingMode = vm;
00155 }
00156 
00157 //-----------------------------------------------------------------------
00158 //-- Turn on/off specular light                                       ---
00159 //-----------------------------------------------------------------------
00160 void vu111211A::setSpecular(int spec)
00161 {
00162     specular = spec;
00163 }
00164 
00165 //-----------------------------------------------------------------------
00166 //-- Turn on/off fast classification                                  ---
00167 //-----------------------------------------------------------------------
00168 void vu111211A::setFastClassification(int fastClass)
00169 {
00170     fastClassification = fastClass;
00171     if (fastClassification == 1) {
00172         initFastClassification();
00173     } else {
00174         runlengthEncode();
00175     }
00176 }
00177 
00178 //-----------------------------------------------------------------------
00179 //-- Get the classification mode                                      ---
00180 //-----------------------------------------------------------------------
00181 int vu111211A::getFastClassification()
00182 {
00183     return fastClassification;
00184 }
00185 
00186 //-----------------------------------------------------------------------
00187 //-- Classifies the volume represented by the octree                  ---
00188 //-----------------------------------------------------------------------
00189 void vu111211A::classify() 
00190 {
00191     fastClassificationObj->classify();
00192 }
00193 
00194 //-----------------------------------------------------------------------
00195 //-- Zoom                                                             ---
00196 //-----------------------------------------------------------------------
00197 void vu111211A::zoom(float value) 
00198 {
00199     if (value > 0) zoomValue = value;
00200 }
00201 
00202 void vu111211A::setOrthogonalWarpOpenGL(int useOpenGL) 
00203 {
00204     orthoWarpOpenGL = useOpenGL;
00205 }
00206 
00207 //-----------------------------------------------------------------------
00208 //-- Set size of canvas                                               ---
00209 //-----------------------------------------------------------------------
00210 void vu111211A::setCanvasSize(int width, int height) 
00211 {
00212     canvasWidth = width;
00213     canvasHeight = height;
00214 }
00215 
00216 //-----------------------------------------------------------------------
00217 //-- Reimplements the read() method to do some extra volume data      ---
00218 //-- processing.                                                      ---
00219 //-----------------------------------------------------------------------
00221 bool vu111211A::read()
00222 {
00223     bool success = vu111211::read();
00224     if (!success) return false;
00225 
00226     //Allocate memory for the normals.
00227     if (m_Normals != 0) delete [] m_Normals;
00228 
00229     m_Normals = new float[m_DataSize*3];
00230 
00231     computeMaxSize();
00232 
00233     //initially value for distance between eye and projection plane
00234     eye_distance = getMinEyeDistance();
00235 
00236     createGLImage();
00237     computeNormals();
00238 
00239     //Allocate memory for intermediate image
00240     intermediate = new intermediatePixel[maxSize*maxSize*4];
00241 
00242     if (fastClassification == 1) {
00243         initFastClassification();
00244     } else {
00245         runlengthEncode();
00246     }
00247 
00248     return true;
00249 }
00250 
00251 //-----------------------------------------------------------------------
00252 //-- Implements readRaw                                               ---
00253 //-----------------------------------------------------------------------
00254 bool vu111211A::readRaw(void)
00255 {
00256     dword len;
00257     ifstream in;
00258 
00259     in.open(m_FileName, ios::in|ios::binary
00260 #ifdef IS_NOCREATE_NEEDED
00261             |ios::nocreate
00262 #endif
00263         );
00264     if (!in.is_open()) return false;
00265 
00266     in.seekg(0, ios::end);
00267     len = in.tellg();
00268     in.seekg(0, ios::beg);
00269 
00270     in >> m_Dim1Size >> m_Dim2Size >> m_Dim3Size;
00271     if (in.fail()) return false;
00272     m_DataSize = m_Dim1Size*m_Dim2Size*m_Dim3Size;
00273 
00274     m_Data = new byte[m_DataSize];
00275     in.read((char*)m_Data, m_DataSize);
00276     if (in.fail()) return false;
00277 
00278     in.close();
00279 
00280     //Allocate memory for the normals.
00281     if (m_Normals != 0) delete [] m_Normals;
00282 
00283     m_Normals = new float[m_DataSize];
00284 
00285     computeNormals();
00286 
00287     return true;
00288 }
00289 
00290 //-----------------------------------------------------------------------
00291 //-- Removes the runlength encoding                                   ---
00292 //-----------------------------------------------------------------------
00293 void vu111211A::removeRunlengthEncoding() 
00294 {
00295     // if neccessary delete existing runlength-encoding
00296     if (dataX) {
00297         for (dword i = 0; i < m_Dim1Size; ++i) { 
00298             if (dataX[i].size>0) {
00299                 delete [] dataX[i].voxel;
00300                 delete [] dataX[i].runlength;
00301             }
00302         }
00303         delete [] dataX;
00304     }
00305     if (dataY) {
00306         for (dword i = 0; i < m_Dim2Size; ++i) { 
00307             if (dataY[i].size>0) {
00308                 delete [] dataY[i].voxel;
00309                 delete [] dataY[i].runlength;
00310             }
00311         }
00312         delete [] dataY;
00313     }
00314     if (dataZ) {
00315         for (dword i = 0; i < m_Dim3Size; ++i) { 
00316             if (dataZ[i].size>0) {
00317                 delete [] dataZ[i].voxel;
00318                 delete [] dataZ[i].runlength;
00319             }
00320         }
00321         delete [] dataZ;
00322     }
00323 }
00324 
00325 //-----------------------------------------------------------------------
00326 //-- Do the runlength encoding                                        ---
00327 //-----------------------------------------------------------------------
00328 void vu111211A::runlengthEncode() 
00329 {
00330     //If there exists already a runlength encoding, remove it
00331     removeRunlengthEncoding();
00332 
00333     dword voxels = 0;           // the number of voxels in the current plane
00334     dword runs = 0;             // the number of runs in the current plane
00335     byte lengthr = 0;           // length of the run
00336     byte state = 0;             // are we in a run of transparent (0) 
00337                                 // or non-transparent (1) voxels?
00338     dword voxelIndex = 0;       // index of voxel in array voxel
00339     dword runIndex = 0;         // index of run in array runlength
00340     dword indexRunlength = 0;   // index in the runlength encoded data
00341     dword indexVolume = 0;      // index in the original volume data
00342     dword maxPlaneSize = m_Dim2Size * m_Dim1Size;  // uncompressed plane size
00343 
00344     dataZ = new RLEplane[m_Dim3Size];
00345 
00346     for (word z = 0; z < m_Dim3Size; ++z) {
00347 
00348         voxels = 0;
00349         runs = 0;
00350         state = 0;
00351         lengthr = 0;
00352 
00353         // find out the number of voxels and runs in the plane
00354         for (dword i = 0; i < maxPlaneSize; ++i) {
00355             if (m_TFunc[m_Data[maxPlaneSize*z + i]][3] > THRESHOLD_RUNLENGTH) {
00356                 ++voxels;
00357                 if (state == 0) {
00358                     ++runs;
00359                     lengthr = 0;
00360                     state = 1;
00361                 }
00362             } else {
00363                 if (state == 1) {
00364                     ++runs;
00365                     lengthr = 0;
00366                     state = 0;
00367                 }
00368             }
00369             if ((lengthr == 255) || ((i+1) % m_Dim1Size == 0)) {
00370                 runs += 3;
00371                 lengthr = 0;
00372                 if (state == 0) state = 1; else state = 0;
00373             }
00374         }
00375 
00376         // set data of the plane
00377         dataZ[z].size = voxels;
00378         dataZ[z].dim1_pos = 0;
00379         dataZ[z].dim2_pos = 0;
00380         dataZ[z].scale = 1;
00381 
00382         // if plane is not empty...
00383         if (voxels > 0) {
00384 
00385             dataZ[z].voxel = new RLEvoxel[voxels+1];
00386             dataZ[z].runlength = new byte[runs+1];
00387 
00388             state = 0;
00389             indexRunlength = 0;
00390             voxelIndex = 0;
00391             runIndex = 0;
00392 
00393             dataZ[z].runlength[0] = 0; // init first runlength
00394 
00395             // construct data for current plane
00396             for (dword i = 0; i < maxPlaneSize; ++i) {
00397 
00398                 if (m_TFunc[m_Data[maxPlaneSize*z + i]][3] > THRESHOLD_RUNLENGTH) {
00399 
00400                     byte val = m_Data[maxPlaneSize*z + i];
00401                     dataZ[z].voxel[voxelIndex].value = val;
00402 
00403                     dataZ[z].voxel[voxelIndex].normal[XDIR] = 
00404                         m_Normals[(maxPlaneSize*z + i)*3 + XDIR];
00405                     dataZ[z].voxel[voxelIndex].normal[YDIR] = 
00406                         m_Normals[(maxPlaneSize*z + i)*3 + YDIR];
00407                     dataZ[z].voxel[voxelIndex].normal[ZDIR] = 
00408                         m_Normals[(maxPlaneSize*z + i)*3 + ZDIR];
00409 
00410                     dataZ[z].voxel[voxelIndex].red = m_TFunc[val][0];
00411                     dataZ[z].voxel[voxelIndex].green = m_TFunc[val][1];
00412                     dataZ[z].voxel[voxelIndex].blue = m_TFunc[val][2];
00413                     dataZ[z].voxel[voxelIndex].opacity = m_TFunc[val][3];
00414 
00415                     ++voxelIndex;
00416                     if (state == 0) {
00417                         ++runIndex;
00418                         dataZ[z].runlength[runIndex] = 0;
00419                         state = 1;
00420                     }
00421                     ++dataZ[z].runlength[runIndex];
00422                 } else {
00423                     if (state == 1) {
00424                         ++runIndex;
00425                         dataZ[z].runlength[runIndex] = 0;
00426                         state = 0;
00427                     }
00428                     ++dataZ[z].runlength[runIndex];
00429                 }
00430                 if ((dataZ[z].runlength[runIndex] == 255) || ((i+1) % m_Dim1Size == 0)) {
00431                     ++runIndex;
00432                     dataZ[z].runlength[runIndex] = 0;
00433                     if (state == 0) state = 1; else state = 0;
00434                 }
00435             }
00436 
00437         } else {
00438             dataZ[z].voxel = 0;
00439             dataZ[z].runlength = 0;
00440         }
00441     }
00442 
00443     //--------------------------------------------------------------
00444     // Now we do the same for the main viewing direction in y
00445 
00446     dword start = 0;
00447     indexVolume = 0;
00448     maxPlaneSize = m_Dim1Size * m_Dim3Size;  // uncompressed plane size
00449 
00450     dataY = new RLEplane[m_Dim2Size];
00451 
00452     for (word y = 0; y < m_Dim2Size; ++y) {
00453 
00454         voxels = 0;
00455         runs = 0;
00456         state = 0;
00457         lengthr = 0;
00458 
00459         dword step = m_Dim1Size * m_Dim2Size;  // size of the steps in a row
00460         dword row_length = (m_Dim3Size-1)*step;  // length of a row
00461 
00462 
00463         indexVolume = y * m_Dim1Size;
00464 
00465         start = indexVolume;
00466 
00467         // find out the number of voxels and runs in the plane
00468         for (dword i = 0; i < maxPlaneSize; ++i) {
00469             if (m_TFunc[m_Data[indexVolume]][3] > THRESHOLD_RUNLENGTH) {
00470                 ++voxels;
00471                 if (state == 0) {
00472                     ++runs;
00473                     lengthr = 0;
00474                     state = 1;
00475                 }
00476             } else {
00477                 if (state == 1) {
00478                     ++runs;
00479                     lengthr = 0;
00480                     state = 0;
00481                 }
00482             }
00483             if ((lengthr == 255) || ((i+1) % m_Dim3Size == 0)) {
00484                 runs += 3;
00485                 lengthr = 0;
00486                 if (state == 0) state = 1; else state = 0;
00487             }
00488 
00489 
00490             if (i > 0) {
00491                 //if one row is finished -> increment start position for next row and
00492                 //continue at this new start position
00493                 if (indexVolume >= y*m_Dim1Size+row_length) {
00494                 ++start;
00495                     indexVolume = start;
00496                 //...otherwise just go on tracing through the row
00497                 } else if (i > 0) {    
00498                     indexVolume += step;
00499                 }
00500             }
00501         }
00502 
00503         // set data of the plane
00504         dataY[y].size = voxels;
00505         dataY[y].dim1_pos = 0;
00506         dataY[y].dim2_pos = 0;
00507         dataY[y].scale = 1;
00508 
00509         // if plane is not empty...
00510         if (voxels > 0) {
00511 
00512             dataY[y].voxel = new RLEvoxel[voxels+1];
00513             dataY[y].runlength = new byte[runs+1];
00514 
00515             state = 0;
00516             indexRunlength = 0;
00517             voxelIndex = 0;
00518             runIndex = 0;
00519 
00520             dataY[y].runlength[0] = 0; // init first runlength
00521 
00522             indexVolume = y * m_Dim1Size;
00523             start = indexVolume;
00524 
00525             // construct data for current plane
00526             for (dword i = 0; i < maxPlaneSize; ++i) {
00527 
00528                 if (i > 0) {
00529                     //if one row is finished -> increment start position for next row and
00530                     //continue at this new start position
00531                     if (indexVolume >= y*m_Dim1Size+row_length) {
00532                         ++start;
00533                         indexVolume = start;
00534                     //...otherwise just go on tracing through the row
00535                     } else if (i > 0) {    
00536                         indexVolume += step;
00537                     }
00538                 }
00539 
00540                 if (m_TFunc[m_Data[indexVolume]][3] > THRESHOLD_RUNLENGTH) {
00541 
00542                     byte val = m_Data[indexVolume];
00543                     dataY[y].voxel[voxelIndex].value = val;
00544                     dataY[y].voxel[voxelIndex].normal[XDIR] = m_Normals[indexVolume*3 + XDIR];
00545                     dataY[y].voxel[voxelIndex].normal[YDIR] = m_Normals[indexVolume*3 + YDIR];
00546                     dataY[y].voxel[voxelIndex].normal[ZDIR] = m_Normals[indexVolume*3 + ZDIR];
00547                     dataY[y].voxel[voxelIndex].red = m_TFunc[val][0];
00548                     dataY[y].voxel[voxelIndex].green = m_TFunc[val][1];
00549                     dataY[y].voxel[voxelIndex].blue = m_TFunc[val][2];
00550                     dataY[y].voxel[voxelIndex].opacity = m_TFunc[val][3];
00551 
00552                     ++voxelIndex;
00553                     if (state == 0) {
00554                         ++runIndex;
00555                         dataY[y].runlength[runIndex] = 0;
00556                         state = 1;
00557                     }
00558                     ++dataY[y].runlength[runIndex];
00559                 } else {
00560                     if (state == 1) {
00561                         ++runIndex;
00562                         dataY[y].runlength[runIndex] = 0;
00563                         state = 0;
00564                     }
00565                     ++dataY[y].runlength[runIndex];
00566                 }
00567 
00568                 //if the capacity of the current runlength is exhausted or if there
00569                 //is a line-break: increase runIndex (new run with length 0) and change 
00570                 //current state 
00571                 if ((dataY[y].runlength[runIndex] == 255) || ((i+1) % m_Dim3Size == 0)) {
00572                     ++runIndex;
00573                     dataY[y].runlength[runIndex] = 0;
00574                     if (state == 0) state = 1; else state = 0;
00575                 }
00576             }
00577 
00578         } else {
00579             dataY[y].voxel = 0;
00580             dataY[y].runlength = 0;
00581         }
00582     }
00583 
00584     //--------------------------------------------------------------
00585     // Now we do the same for the main viewing direction in x
00586 
00587     indexVolume = 0;
00588     maxPlaneSize = m_Dim2Size * m_Dim3Size;  // uncompressed plane size
00589 
00590     dataX = new RLEplane[m_Dim1Size];
00591 
00592     for (word x = 0; x < m_Dim1Size; ++x) {
00593 
00594         voxels = 0;
00595         runs = 0;
00596         state = 0;
00597         lengthr = 0;
00598 
00599         indexVolume = x;
00600 
00601         // find out the number of voxels and runs in the plane
00602         for (dword i = 0; i < maxPlaneSize; ++i) {
00603             if (m_TFunc[m_Data[indexVolume]][3] > THRESHOLD_RUNLENGTH) {
00604                 ++voxels;
00605                 if (state == 0) {
00606                     ++runs;
00607                     lengthr = 0;
00608                     state = 1;
00609                 }
00610             } else {
00611                 if (state == 1) {
00612                     ++runs;
00613                     lengthr = 0;
00614                     state = 0;
00615                 }
00616             }
00617             if ((lengthr == 255) || ((i+1) % m_Dim2Size == 0)) {
00618                 runs += 3;
00619                 lengthr = 0;
00620                 if (state == 0) state = 1; else state = 0;
00621             }
00622 
00623             indexVolume += m_Dim1Size;
00624         }
00625 
00626         // set data of the plane
00627         dataX[x].size = voxels;
00628         dataX[x].dim1_pos = 0;
00629         dataX[x].dim2_pos = 0;
00630         dataX[x].scale = 1;
00631 
00632         // if plane is not empty...
00633         if (voxels > 0) {
00634 
00635             dataX[x].voxel = new RLEvoxel[voxels+1];
00636             dataX[x].runlength = new byte[runs+1];
00637 
00638             state = 0;
00639             indexRunlength = 0;
00640             voxelIndex = 0;
00641             runIndex = 0;
00642 
00643             dataX[x].runlength[0] = 0; // init first runlength
00644 
00645             indexVolume = x;
00646 
00647             // construct data for current plane
00648             for (dword i = 0; i < maxPlaneSize; ++i) {
00649 
00650                 if (m_TFunc[m_Data[indexVolume]][3] > THRESHOLD_RUNLENGTH) {
00651 
00652                     byte val = m_Data[indexVolume];
00653                     dataX[x].voxel[voxelIndex].value = val;
00654                     dataX[x].voxel[voxelIndex].normal[XDIR] = m_Normals[indexVolume*3 + XDIR];
00655                     dataX[x].voxel[voxelIndex].normal[YDIR] = m_Normals[indexVolume*3 + YDIR];
00656                     dataX[x].voxel[voxelIndex].normal[ZDIR] = m_Normals[indexVolume*3 + ZDIR];
00657 
00658                     dataX[x].voxel[voxelIndex].red = m_TFunc[val][0];
00659                     dataX[x].voxel[voxelIndex].green = m_TFunc[val][1];
00660                     dataX[x].voxel[voxelIndex].blue = m_TFunc[val][2];
00661                     dataX[x].voxel[voxelIndex].opacity = m_TFunc[val][3];
00662 
00663                     ++voxelIndex;
00664                     if (state == 0) {
00665                         ++runIndex;
00666                         dataX[x].runlength[runIndex] = 0;
00667                         state = 1;
00668                     }
00669                     ++dataX[x].runlength[runIndex];
00670                 } else {
00671                     if (state == 1) {
00672                         ++runIndex;
00673                         dataX[x].runlength[runIndex] = 0;
00674                         state = 0;
00675                     }
00676                     ++dataX[x].runlength[runIndex];
00677                 }
00678                 if ((dataX[x].runlength[runIndex] == 255) || ((i+1) % m_Dim2Size == 0)) {
00679                     ++runIndex;
00680                     dataX[x].runlength[runIndex] = 0;
00681                     if (state == 0) state = 1; else state = 0;
00682                 }
00683 
00684                 indexVolume += m_Dim1Size;
00685             }
00686 
00687         } else {
00688             dataX[x].voxel = 0;
00689             dataX[x].runlength = 0;
00690         }
00691     }
00692 }
00693 
00694 //-----------------------------------------------------------------------
00695 //-- Preprocesses volume data for rendering once it's been read       ---
00696 //-----------------------------------------------------------------------
00697 void vu111211A::computeNormals(void)
00698 {
00699     float norm[3];
00700     float len;
00701     dword w, wh;
00702     dword i, j, k, index;
00703 
00704     w = m_Dim1Size;
00705     wh = m_Dim1Size*m_Dim2Size;
00706 
00707     //
00708     // Compute all of the normals
00709     //
00710     index = 0;
00711 
00712     for(k=0;k<m_Dim3Size;++k)
00713     {
00714         for(j=0;j<m_Dim2Size;++j)
00715         {
00716             for(i=0;i<m_Dim1Size;++i)
00717             {
00718                 if (i == 0) {
00719                     norm[0] = m_Data[index+1]-m_Data[index];
00720                 }
00721                 else if (i == m_Dim1Size-1) {
00722                     norm[0] = m_Data[index]-m_Data[index-1];
00723                 }
00724                 else {
00725                     norm[0] = (m_Data[index+1]-m_Data[index-1])*1;
00726 
00727                     if ((j > 0) && (j < m_Dim2Size) && (k > 0) && (k < m_Dim3Size)) {
00728                     norm[0] += (m_Data[index  +w +1]-m_Data[index  +w -1])*0.2;
00729                     norm[0] += (m_Data[index  -w +1]-m_Data[index  -w -1])*0.2;
00730                     norm[0] += (m_Data[index +wh +1]-m_Data[index +wh -1])*0.2;
00731                     norm[0] += (m_Data[index -wh +1]-m_Data[index -wh -1])*0.2;
00732 
00733                     norm[0] += (m_Data[index +wh +w +1]-m_Data[index +wh +w -1])*0.1;
00734                     norm[0] += (m_Data[index +wh -w +1]-m_Data[index +wh -w -1])*0.1;
00735                     norm[0] += (m_Data[index -wh +w +1]-m_Data[index -wh +w -1])*0.1;
00736                     norm[0] += (m_Data[index -wh -w +1]-m_Data[index -wh -w -1])*0.1;
00737                     }
00738                 }
00739 
00740                 if (j == 0) {
00741                     norm[1] = m_Data[index+w]-m_Data[index];
00742                 }
00743                 else if (j == m_Dim2Size-1) {
00744                     norm[1] = m_Data[index]-m_Data[index-w];
00745                 }
00746                 else {
00747                     norm[1] = (m_Data[index+w]-m_Data[index-w])*1;
00748 
00749                     if ((i > 0) && (i < m_Dim1Size) && (k > 0) && (k < m_Dim3Size)) {
00750                     norm[1] += (m_Data[index  +1 +w]-m_Data[index +1  -w])*0.2;
00751                     norm[1] += (m_Data[index  -1 +w]-m_Data[index -1  -w])*0.2;
00752                     norm[1] += (m_Data[index +wh +w]-m_Data[index +wh -w])*0.2;
00753                     norm[1] += (m_Data[index -wh +w]-m_Data[index -wh -w])*0.2; 
00754 
00755                     norm[1] += (m_Data[index +wh +1 +w]-m_Data[index +wh +1 -w])*0.1;
00756                     norm[1] += (m_Data[index +wh -1 +w]-m_Data[index +wh -1 -w])*0.1;
00757                     norm[1] += (m_Data[index -wh +1 +w]-m_Data[index -wh +1 -w])*0.1;
00758                     norm[1] += (m_Data[index -wh -1 +w]-m_Data[index -wh -1 -w])*0.1;
00759                     }
00760                 }
00761 
00762                 if (k == 0) {
00763                     norm[2] = m_Data[index+wh]-m_Data[index];
00764                 }
00765                 else if (k == m_Dim3Size-1) {
00766                     norm[2] = m_Data[index]-m_Data[index-wh];
00767                 }
00768                 else {
00769                     norm[2] = (m_Data[index+wh]-m_Data[index-wh])*1;
00770 
00771                     if ((i > 0) && (i < m_Dim1Size) && (j > 0) && (j < m_Dim2Size)) {
00772                     norm[2] += (m_Data[index +1 +wh]-m_Data[index +1 -wh])*0.2;
00773                     norm[2] += (m_Data[index -1 +wh]-m_Data[index -1 -wh])*0.2;
00774                     norm[2] += (m_Data[index +w +wh]-m_Data[index +w -wh])*0.2;
00775                     norm[2] += (m_Data[index -w +wh]-m_Data[index -w -wh])*0.2;
00776 
00777                     norm[2] += (m_Data[index +w +1 +wh]-m_Data[index +w +1 -wh])*0.1;
00778                     norm[2] += (m_Data[index +w -1 +wh]-m_Data[index +w -1 -wh])*0.1;
00779                     norm[2] += (m_Data[index -w +1 +wh]-m_Data[index -w +1 -wh])*0.1;
00780                     norm[2] += (m_Data[index -w -1 +wh]-m_Data[index -w -1 -wh])*0.1;
00781                     }
00782                 }
00783 
00784                 len = (float)sqrt((double)((norm[0]*norm[0])+
00785                                            (norm[1]*norm[1])+
00786                                            (norm[2]*norm[2])));
00787 
00788                 if (len > 0.0f)
00789                 {
00790                     norm[0] /= len;
00791                     norm[1] /= len;
00792                     norm[2] /= len;
00793                 }
00794 
00795                 m_Normals[index*3+XDIR] = norm[0];
00796                 m_Normals[index*3+YDIR] = norm[1];
00797                 m_Normals[index*3+ZDIR] = norm[2];
00798 
00799                 ++index;
00800             }
00801         }
00802     }
00803 }
00804 
00805 //-----------------------------------------------------------------------
00806 //-- Initializes the Fast Classification                              ---
00807 //-----------------------------------------------------------------------
00808 void vu111211A::initFastClassification() 
00809 {
00810     removeRunlengthEncoding();
00811 
00812     //Set up the Arrays for saving the sheared positions of the planes
00813     dataX = new RLEplane[m_Dim1Size];
00814     for (unsigned int i = 0; i < m_Dim1Size; ++i) dataX[i].size = 0;
00815     dataY = new RLEplane[m_Dim2Size];
00816     for (unsigned int i = 0; i < m_Dim2Size; ++i) dataY[i].size = 0;
00817     dataZ = new RLEplane[m_Dim3Size];
00818     for (unsigned int i = 0; i < m_Dim3Size; ++i) dataZ[i].size = 0;
00819 
00820     //init the FastClassification object
00821     fastClassificationObj = new FastClassification(m_Data,m_Dim1Size,m_Dim2Size,m_Dim3Size);
00822 }
00823 
00824 //-----------------------------------------------------------------------
00825 //-- Determines the main viewing direction from the viewing vector    ---
00826 //-----------------------------------------------------------------------
00827 void vu111211A::computeMainViewingDirection() {
00828     mainViewingDir = XDIR;
00829 
00830     double vXsq = m_View[XDIR]*m_View[XDIR];
00831     double vYsq = m_View[YDIR]*m_View[YDIR];
00832     double vZsq = m_View[ZDIR]*m_View[ZDIR];
00833 
00834     if ((vXsq >= vYsq) && (vXsq >= vZsq)) {        // main viewing direction: X
00835         mainViewingDir = XDIR;
00836     } else if ((vYsq >= vXsq) && (vYsq >=vZsq)) {  // main viewing direction: Y
00837         mainViewingDir = YDIR;
00838     } else {                                       // main viewing direction Z
00839         mainViewingDir = ZDIR;
00840     }
00841 }
00842 
00843 //-----------------------------------------------------------------------
00844 //-- Determines the permutation matrix                                ---
00845 //-----------------------------------------------------------------------
00846 void vu111211A::computePermutationMatrix() 
00847 {
00848   // depending on the main viewing direction a permutation matrix is choosen
00849   switch(mainViewingDir) {
00850       case XDIR: {
00851           permMatrix[0][0] = 0; permMatrix[0][1] = 1; permMatrix[0][2] = 0; permMatrix[0][3] = 0;
00852           permMatrix[1][0] = 0; permMatrix[1][1] = 0; permMatrix[1][2] = 1; permMatrix[1][3] = 0;
00853           permMatrix[2][0] = 1; permMatrix[2][1] = 0; permMatrix[2][2] = 0; permMatrix[2][3] = 0;
00854           permMatrix[3][0] = 0; permMatrix[3][1] = 0; permMatrix[3][2] = 0; permMatrix[3][3] = 1;
00855           break;
00856       }
00857       case YDIR: {
00858           permMatrix[0][0] = 0; permMatrix[0][1] = 0; permMatrix[0][2] = 1; permMatrix[0][3] = 0;
00859           permMatrix[1][0] = 1; permMatrix[1][1] = 0; permMatrix[1][2] = 0; permMatrix[1][3] = 0;
00860           permMatrix[2][0] = 0; permMatrix[2][1] = 1; permMatrix[2][2] = 0; permMatrix[2][3] = 0;
00861           permMatrix[3][0] = 0; permMatrix[3][1] = 0; permMatrix[3][2] = 0; permMatrix[3][3] = 1;
00862           break;
00863       }
00864       case ZDIR: {
00865           permMatrix[0][0] = 1; permMatrix[0][1] = 0; permMatrix[0][2] = 0; permMatrix[0][3] = 0;
00866           permMatrix[1][0] = 0; permMatrix[1][1] = 1; permMatrix[1][2] = 0; permMatrix[1][3] = 0;
00867           permMatrix[2][0] = 0; permMatrix[2][1] = 0; permMatrix[2][2] = 1; permMatrix[2][3] = 0;
00868           permMatrix[3][0] = 0; permMatrix[3][1] = 0; permMatrix[3][2] = 0; permMatrix[3][3] = 1;
00869           break;
00870       }
00871   }
00872 }
00873 
00874 
00875 //-----------------------------------------------------------------------
00876 //-- Computes the viewing matrix according to the given view- and     ---
00877 //-- up-vectors.                                                      ---
00878 //-----------------------------------------------------------------------
00879 void vu111211A::computeViewMatrix(void)
00880 {
00881   vuVector puffer;
00882   vuVector cross;
00883   vuVector up = m_Up;
00884   vuVector v = m_View;
00885 
00886   double to_degree = 180/M_PI;  // for dealing with angles in degree
00887   double alpha = 0;
00888   double beta = 0;
00889 
00890   vuMatrix rotate_view;  // a rotation matrix to bring the 
00891                          // view-vector in correct position 
00892 
00893   vuMatrix rotate_up;    // a rotation matrix to bring the 
00894                          // up-vector in correct position
00895 
00896   // init the matrices with the identity
00897   viewMatrix.makeIdentity();
00898   rotate_view.makeIdentity();
00899   rotate_up.makeIdentity();
00900 
00901   // unitVecMain ist the unit vector pointing in the direction of the 
00902   // main viewing direction which is Z
00903   vuVector unitVecMain;
00904   unitVecMain[0] = 0;
00905   unitVecMain[1] = 0;
00906   unitVecMain[2] = 1;
00907   unitVecMain[3] = 1;
00908 
00909   // compute the standard object coordinates for view- and up-vector
00910   puffer = v;
00911   v = puffer * permMatrix;
00912   puffer = up;
00913   up = puffer * permMatrix;
00914 
00915   // compute the direction-vector of the axis for the following
00916   // rotation (= v x unitVecMain)
00917   cross = v.cross(unitVecMain);
00918   cross.makeUnit();
00919 
00920   // compute angle between view-vector and unit-vector of main viewing
00921   // direction and set step2-matrix to the according rotation around cross
00922   v.makeUnit();
00923   alpha = -to_degree * acos(v[ZDIR]);
00924   rotate_view.makeRotate(cross,alpha);
00925 
00926   // --> now step2 represents the matrix that turns the viewing vector 
00927   //     into exactly the same direction as the main viewing vector
00928 
00929   // compute the value of the up vector when it is transformed by steps 1 and 2
00930   puffer = up;
00931   up = puffer * rotate_view;
00932 
00933   // compute angle between up-vector (without main-component)
00934   // and set step2-matrix to the according rotation around the main viewing axis
00935   up[ZDIR] = 0;
00936   up.makeUnit();
00937 
00938   // compute the angle for rotation to bring the up vector in correct position
00939   beta = to_degree * acos(up[YDIR]);
00940   if (up[XDIR] >= 0) beta *= (-1);
00941   rotate_up.makeRotateZ(beta);
00942 
00943   // compute final view-matrix and return the value
00944   viewMatrix = rotate_view * rotate_up;      // in [Lacroute] this is M'
00945   worldMatrix = permMatrix * viewMatrix;
00946 
00947   if (viewingMode == VIEWING_MODE_PERSPECTIVE) {
00948       projMatrix[0][0] = eye_distance; projMatrix[0][1] = 0; projMatrix[0][2] = 0; projMatrix[0][3] = 0;
00949       projMatrix[1][0] = 0; projMatrix[1][1] = eye_distance; projMatrix[1][2] = 0; projMatrix[1][3] = 0;
00950       projMatrix[2][0] = 0; projMatrix[2][1] = 0; projMatrix[2][2] = 1; projMatrix[2][3] = 0;
00951       projMatrix[3][0] = 0; projMatrix[3][1] = 0; projMatrix[3][2] = 1; projMatrix[3][3] = eye_distance;
00952       viewMatrix = worldMatrix * projMatrix;
00953   }
00954 }
00955 
00956 //-----------------------------------------------------------------------
00957 //-- Determines width, height and depth of the Volume and saves these ---
00958 //-- values in volumeWidth, volumeHeight and volumeDepth.             ---
00959 //-----------------------------------------------------------------------
00960 void vu111211A::computeRunlengthSizes()
00961 {
00962     switch(mainViewingDir) {
00963         case XDIR: {
00964             volumeDepth = m_Dim1Size;
00965             volumeWidth = m_Dim2Size;
00966             volumeHeight = m_Dim3Size;
00967             break;
00968         }
00969         case YDIR: {
00970             volumeDepth = m_Dim2Size;  
00971             volumeWidth = m_Dim3Size;  
00972             volumeHeight = m_Dim1Size; 
00973             break;
00974         }
00975         case ZDIR: {
00976             volumeDepth = m_Dim3Size;
00977             volumeWidth = m_Dim1Size;
00978             volumeHeight = m_Dim2Size;
00979             break;
00980         }
00981     }
00982 }
00983 
00984 //-----------------------------------------------------------------------
00985 //-- Computes the Shear Factors si, sj and the translation ti, tj     ---
00986 //-- for orthogonal projection.                                       ---
00987 //-----------------------------------------------------------------------
00988 void vu111211A::computeShearAndWarpFactorsOrtho()
00989 {
00990     if ((viewMatrix[0][0]*viewMatrix[1][1]-viewMatrix[1][0]*viewMatrix[0][1]) == 0) si = 1; 
00991     else si = ((viewMatrix[1][1]*viewMatrix[0][2]-viewMatrix[0][1]*viewMatrix[1][2])/
00992                (viewMatrix[0][0]*viewMatrix[1][1]-viewMatrix[1][0]*viewMatrix[0][1]));
00993     if ((viewMatrix[0][0]*viewMatrix[1][1]-viewMatrix[1][0]*viewMatrix[0][1]) == 0) sj = 1;
00994     else sj = ((viewMatrix[0][0]*viewMatrix[1][2]-viewMatrix[1][0]*viewMatrix[0][2])/
00995                (viewMatrix[0][0]*viewMatrix[1][1]-viewMatrix[1][0]*viewMatrix[0][1]));
00996 
00997     if (si >= 0) ti = 0.0; else ti = (-si*volumeDepth);
00998     if (sj >= 0) tj = 0.0; else tj = (-sj*volumeDepth);
00999 }
01000 
01001 //-----------------------------------------------------------------------
01002 //-- Shear the planes of the volume for orthogonal projection         ---
01003 //-----------------------------------------------------------------------
01004 void vu111211A::shearOrtho()
01005 {
01006     // normalize the viewing vector
01007     m_View.makeUnit();
01008 
01009     double shearPos_i = 0;
01010     double shearPos_j = 0;
01011 
01012     if (orthoWarpOpenGL == 1) {
01013         shearPos_i = ti;
01014         shearPos_j = tj;
01015     } else {
01016         shearPos_i = maxSize - ((volumeWidth  + si * volumeDepth) / 2);
01017         shearPos_j = maxSize - ((volumeHeight + sj * volumeDepth) / 2);
01018     }
01019 
01020     // select the used dataset...
01021     switch(mainViewingDir) {
01022         case ZDIR: {
01023             curr = dataZ;
01024             break;
01025         }
01026         case YDIR: {
01027             curr = dataY;
01028             break;
01029         }
01030         case XDIR: {
01031             curr = dataX;
01032             break;
01033         }
01034     }
01035 
01036     if (m_View[mainViewingDir] >= 0) direction = 0; else direction = 1;
01037 
01038     for (dword i = 0; i < volumeDepth; ++i) {
01039         curr[i].dim1_pos = shearPos_i;
01040         curr[i].dim2_pos = shearPos_j;
01041         shearPos_i += si;
01042         shearPos_j += sj;
01043         if (shearPos_i < 0) shearPos_i = 0;
01044         if (shearPos_j < 0) shearPos_j = 0;
01045     }
01046 }
01047 
01048 //-----------------------------------------------------------------------
01049 //-- Computes the length of the view vector normalized in Z           ---
01050 //-----------------------------------------------------------------------
01051 float vu111211A::computeZNormalizedViewingVectorLength() {
01052     //precompute length of z-normalized viewing vector (later used for opacity correction)
01053     double z,x,y;
01054     m_View.makeUnit();
01055     if      (mainViewingDir == XDIR) { z = m_View[XDIR]; y = m_View[YDIR]; x = m_View[ZDIR]; }
01056     else if (mainViewingDir == YDIR) { z = m_View[YDIR]; y = m_View[ZDIR]; x = m_View[XDIR]; }
01057     else                          { z = m_View[ZDIR]; y = m_View[XDIR]; x = m_View[YDIR]; }
01058     return sqrt((x/z)*(x/z)+(y/z)*(y/z)+1);
01059 }
01060 
01061 //-----------------------------------------------------------------------
01062 //-- Reset the intermediate image                                     ---
01063 //-----------------------------------------------------------------------
01064 void vu111211A::resetIntermediateImage(int intermediateSize) {
01065     //reset intermediate image
01066     for (int i = 0; i < intermediateSize; ++i) {
01067         intermediate[i].trans = 1;
01068         intermediate[i].red = 0;
01069         intermediate[i].green = 0;
01070         intermediate[i].blue = 0;
01071         intermediate[i].offset = 0;
01072     }
01073 }
01074 
01075 //-----------------------------------------------------------------------
01076 //-- Pre-classify the voxels                                          ---
01077 //-----------------------------------------------------------------------
01078 void vu111211A::preClassification(RLEplane current, float viewVectorLength) {
01079 
01080     // do the shading and the opacity correction
01081     for (dword k = 0; k < current.size; ++k) {
01082         current.voxel[k].normal.makeUnit();
01083 
01084         float diffuse = m_View.dot(current.voxel[k].normal);
01085 
01086         float spec = 0;
01087 
01088         if (diffuse <= 0) { // <= 0
01089             current.voxel[k].red_shaded = 0.0;
01090             current.voxel[k].green_shaded = 0.0;
01091             current.voxel[k].blue_shaded = 0.0;
01092             current.voxel[k].opacity_corr = 0.0;
01093         } else {
01094             current.voxel[k].opacity_corr = 
01095                 1 - pow((1-current.voxel[k].opacity),viewVectorLength);
01096 
01097             if (specular == 1) {
01098                 spec = pow(diffuse,100);
01099             } else {
01100                 specular = 0;
01101             }
01102 
01103             float fac = current.voxel[k].opacity_corr * (0.1 + 0.9 * diffuse);
01104 
01105             current.voxel[k].red_shaded   = current.voxel[k].red   * fac;
01106             current.voxel[k].green_shaded = current.voxel[k].green * fac;
01107             current.voxel[k].blue_shaded  = current.voxel[k].blue  * fac;
01108 
01109             current.voxel[k].red_shaded   += spec;
01110             current.voxel[k].green_shaded += spec;
01111             current.voxel[k].blue_shaded  += spec;
01112             current.voxel[k].opacity_corr += spec;
01113 
01114             if (current.voxel[k].red_shaded   > 1.0) current.voxel[k].red_shaded   = 1.0;
01115             if (current.voxel[k].green_shaded > 1.0) current.voxel[k].green_shaded = 1.0;
01116             if (current.voxel[k].blue_shaded  > 1.0) current.voxel[k].blue_shaded  = 1.0;
01117             if (current.voxel[k].opacity_corr > 1.0) current.voxel[k].opacity_corr = 1.0;
01118         }
01119     }
01120 }
01121 
01122 //-----------------------------------------------------------------------
01123 //-- Init the start-position of the first pointer                     ---
01124 //-----------------------------------------------------------------------
01125 void vu111211A::initRunPtr1(RLEplane current, int &volumeIndex, 
01126                             int &voxelPtr1, int &runPtr1, byte &state1) {
01127     do {
01128         volumeIndex += current.runlength[runPtr1]; 
01129         if (state1 == 0) {
01130             state1 = 1;
01131         } else {
01132             state1 = 0;
01133             voxelPtr1 += current.runlength[runPtr1];
01134         }
01135         ++runPtr1;
01136     } while ((dword)volumeIndex < volumeWidth);
01137 }
01138 
01139 //-----------------------------------------------------------------------
01140 //-- Skip runs with zero length                                       ---
01141 //-----------------------------------------------------------------------
01142 void vu111211A::skipZeroLengthRuns(RLEplane slice, int &runPtr, 
01143                                    byte &count, byte &state) {
01144     while (slice.runlength[runPtr] == 0) {
01145         ++runPtr;
01146         count = 0;
01147         if (state == 0) state = 1; else state = 0;
01148     }
01149 }
01150 
01151 //-----------------------------------------------------------------------
01152 //-- Skip transparent voxels                                          ---
01153 //-----------------------------------------------------------------------
01154 int vu111211A::skipTransparentVoxels(RLEplane slice, 
01155                                      int runPtr1, byte &count1, 
01156                                      int runPtr2, byte &count2, 
01157                                      int &add, 
01158                                      int intermediateWidth, int &pos) {
01159 
01160     int skip;
01161 
01162     if (slice.runlength[runPtr1]-count1 < slice.runlength[runPtr2]-count2) {
01163         skip = slice.runlength[runPtr1]-count1-1;
01164     } else {
01165         skip = slice.runlength[runPtr2]-count2-1;
01166     }
01167 
01168     add += skip;
01169     if ((dword)add >= volumeWidth) {
01170         pos += (intermediateWidth - volumeWidth);
01171         add = add-volumeWidth;
01172         pos += skip-add;
01173     } else {
01174         pos += skip;
01175     }
01176 
01177     count1 += skip;
01178     count2 += skip;
01179 
01180     return skip;
01181 }
01182 
01183 //-----------------------------------------------------------------------
01184 //-- Interpolate a Pixel                                              ---
01185 //-----------------------------------------------------------------------
01186 void vu111211A::interpolatePixel(RLEplane slice,
01187                                  intermediatePixel &pixel,
01188                                  float w1, float w2, float w3, float w4, int volumeIndex,
01189                                  int &voxelPtr1, int runPtr1, int count1, int state1,
01190                                  int &voxelPtr2, int runPtr2, int count2, int state2) {
01191 
01192     double red     = 0;
01193     double green   = 0;
01194     double blue    = 0;
01195     double opacity = 0;
01196 
01197     if ((dword)voxelPtr1 < slice.size) {
01198 
01199         //does state1 point to a non-transparent voxel?
01200         if (state1 == 1) {
01201 
01202             red     += w1*slice.voxel[voxelPtr1].red_shaded;
01203             green   += w1*slice.voxel[voxelPtr1].green_shaded;
01204             blue    += w1*slice.voxel[voxelPtr1].blue_shaded;
01205             opacity += w1*slice.voxel[voxelPtr1].opacity_corr;
01206 
01207             if ((slice.runlength[runPtr1]-count1 > 1) &&
01208                 ((dword)voxelPtr1+1 < slice.size)) {
01209 
01210                 red     += w2*slice.voxel[voxelPtr1+1].red_shaded;
01211                 green   += w2*slice.voxel[voxelPtr1+1].green_shaded;
01212                 blue    += w2*slice.voxel[voxelPtr1+1].blue_shaded;
01213                 opacity += w2*slice.voxel[voxelPtr1+1].opacity_corr;
01214 
01215             }
01216             ++voxelPtr1;
01217         } else {
01218             if ((slice.runlength[runPtr1]-count1 == 1) &&   //non-transp. run ends after ptr1?
01219                 (slice.runlength[runPtr1+1]!=0) &&          //no zero-length run after current?
01220                 ((volumeIndex+1) % volumeWidth)!=0) {   //no line break?
01221 
01222                 red     += w2*slice.voxel[voxelPtr1].red_shaded;  // vorher voxelPtr1+1
01223                 green   += w2*slice.voxel[voxelPtr1].green_shaded;
01224                 blue    += w2*slice.voxel[voxelPtr1].blue_shaded;
01225                 opacity += w2*slice.voxel[voxelPtr1].opacity_corr;
01226             }
01227         }
01228     }
01229 
01230     if ((dword)voxelPtr2 < slice.size) {
01231         if (state2 == 1) {
01232 
01233             red     += w3*slice.voxel[voxelPtr2].red_shaded;
01234             green   += w3*slice.voxel[voxelPtr2].green_shaded;
01235             blue    += w3*slice.voxel[voxelPtr2].blue_shaded;
01236             opacity += w3*slice.voxel[voxelPtr2].opacity_corr;
01237 
01238             if ((slice.runlength[runPtr2]-count2 > 1) && 
01239                 ((dword)voxelPtr2+1 < slice.size)) {
01240 
01241                 red     += w4*slice.voxel[voxelPtr2+1].red_shaded;
01242                 green   += w4*slice.voxel[voxelPtr2+1].green_shaded;
01243                 blue    += w4*slice.voxel[voxelPtr2+1].blue_shaded;
01244                 opacity += w4*slice.voxel[voxelPtr2+1].opacity_corr;
01245 
01246             }
01247             ++voxelPtr2;
01248         } else {
01249             if ((slice.runlength[runPtr2]-count2 == 1) &&  //non-transp. run ends after ptr1?
01250                 (slice.runlength[runPtr2+1]!=0) &&         //no zero-length run after current?
01251                 ((volumeIndex+1) % volumeWidth)!=0) {  //no line break?
01252 
01253                 red     += w4*slice.voxel[voxelPtr2].red_shaded;  //vorher voxelPtr2+1
01254                 green   += w4*slice.voxel[voxelPtr2].green_shaded;
01255                 blue    += w4*slice.voxel[voxelPtr2].blue_shaded;
01256                 opacity += w4*slice.voxel[voxelPtr2].opacity_corr;
01257 
01258             }
01259         }
01260     }
01261 
01262     if (opacity > 0) {
01263 
01264         pixel.offset = 1;
01265 
01266         double fac = pixel.trans;
01267 
01268         pixel.red   += red * fac;
01269         pixel.green += green * fac;
01270         pixel.blue  += blue * fac;
01271 
01272         pixel.trans *= (1-opacity);
01273     }
01274 }
01275 
01276 //-----------------------------------------------------------------------
01277 //-- Skip opaque pixel                                                ---
01278 //-----------------------------------------------------------------------
01279 void vu111211A::skipOpaquePixel(byte state1,int &voxelPtr1,byte state2,int &voxelPtr2) {
01280 
01281     if (state1 == 1) {
01282         ++voxelPtr1;
01283     }
01284 
01285     if (state2 == 1) {
01286         ++voxelPtr2;
01287     }
01288 }
01289 
01290 //-----------------------------------------------------------------------
01291 //-- Goto next voxel                                                  ---
01292 //-----------------------------------------------------------------------
01293 void vu111211A::volumeNext(RLEplane slice, int &volumeIndex, int &add, 
01294                            int intermediateWidth, int &pos,
01295                            int &runPtr1, byte &count1, byte &state1,
01296                            int &runPtr2, byte &count2, byte &state2) {
01297 
01298     //one step foreward through the volume
01299     ++volumeIndex;
01300 
01301     ++add;
01302     if ((dword)add >= volumeWidth) {
01303         pos += (intermediateWidth - volumeWidth);
01304         add = add-volumeWidth;
01305         ++pos;
01306     } else {
01307         ++pos;
01308     }
01309 
01310     ++count1;
01311     ++count2;
01312 
01313     if (count1 >= slice.runlength[runPtr1]) {
01314         ++runPtr1;
01315         count1 = 0;
01316         if (state1 == 1) state1 = 0; else state1 = 1;
01317     }
01318 
01319     if (count2 >= slice.runlength[runPtr2]) {
01320         ++runPtr2;
01321         count2 = 0;
01322         if (state2 == 1) state2 = 0; else state2 = 1;
01323     }
01324 }
01325 
01326 //-----------------------------------------------------------------------
01327 //-- compute the inverse warp matrix for orthogonal viewing           ---
01328 //-----------------------------------------------------------------------
01329 void vu111211A::computeInvWarpOrtho() 
01330 {
01331     // the 2D-warp-Matrix
01332     float m00 = viewMatrix[0][0];
01333     float m01 = viewMatrix[0][1];
01334     float m10 = viewMatrix[1][0];
01335     float m11 = viewMatrix[1][1];
01336 
01337     // cmompute the inverse of the 2D-warp-Matrix
01338     float invDet = 1 / (m00 * m11 - (m01 * m10));
01339 
01340     invWarpOrtho00 = invDet * m11;
01341     invWarpOrtho01 = - invDet * m01;
01342     invWarpOrtho10 = - invDet * m10;
01343     invWarpOrtho11 = invDet * m00;
01344 }
01345 
01346 //-----------------------------------------------------------------------
01347 //-- inverse orthogonal warping of a single point                     ---
01348 //-----------------------------------------------------------------------
01349 void vu111211A::warpOrthoInv(float x, float y, float &x_result, float &y_result)
01350 {
01351     // apply the inverse warp to the given coordinates
01352     x_result = invWarpOrtho00 * x + invWarpOrtho01 * y;
01353     y_result = invWarpOrtho10 * x + invWarpOrtho11 * y;
01354 }
01355 
01356 //-----------------------------------------------------------------------
01357 //-- orthogonal warping                                               ---
01358 //-----------------------------------------------------------------------
01359 void vu111211A::warpOrtho()
01360 {
01361     int glSize = glImageWidth * glImageHeight;
01362     int glHalf = glImageWidth / 2;
01363 
01364     int intermediateWidth = maxSize * 2;
01365 
01366     float int_x = 0.0f;
01367     float int_y = 0.0f;
01368 
01369     float gl_x = 0.0f;
01370     float gl_y = 0.0f;
01371 
01372     float maxSizef = (float)maxSize-1;
01373 
01374     int int_index = 0;
01375     int gl_index  = 0;
01376     while(gl_index < glSize) {
01377 
01378         gl_x = (gl_index % glImageWidth) - glHalf;
01379         gl_y = (gl_index / glImageWidth) - glHalf;
01380 
01381         warpOrthoInv(gl_x, gl_y, int_x, int_y);
01382 
01383         float d1 = (int_x - floor(int_x));
01384         float d2 = (int_y - floor(int_y));
01385 
01386         float w1 = 1-d1-d2+d1*d2;
01387         float w2 = d1-d1*d2;
01388         float w3 = d2-d1*d2;
01389         float w4 = d1*d2;
01390 
01391         int_index = (int)(floor(int_x + maxSizef) + 
01392                           floor(int_y + maxSizef) * intermediateWidth);
01393 
01394 
01395         if ((int_x > -maxSizef) && (int_x < maxSizef - 1) &&
01396             (int_y > -maxSizef) && (int_y < maxSizef - 1)) {
01397 
01398 
01399             float red = intermediate[int_index].red*w1 + 
01400                     intermediate[int_index+1].red*w2 +
01401                     intermediate[int_index+intermediateWidth].red*w3 +
01402                     intermediate[int_index+intermediateWidth+1].red*w4;
01403             float green = intermediate[int_index].green*w1 +
01404                     intermediate[int_index+1].green*w2 +
01405                     intermediate[int_index+intermediateWidth].green*w3 +
01406                     intermediate[int_index+intermediateWidth+1].green*w4;
01407             float blue = intermediate[int_index].blue*w1 + 
01408                     intermediate[int_index+1].blue*w2 +
01409                     intermediate[int_index+intermediateWidth].blue*w3 +
01410                     intermediate[int_index+intermediateWidth+1].blue*w4;
01411 
01412 
01413             glImage[gl_index*4+0] = (GLubyte)(red*255);
01414             glImage[gl_index*4+1] = (GLubyte)(green*255);
01415             glImage[gl_index*4+2] = (GLubyte)(blue*255);
01416             glImage[gl_index*4+3] = 255;
01417         } else {
01418             glImage[gl_index*4+0] = 0;
01419             glImage[gl_index*4+1] = 0;
01420             glImage[gl_index*4+2] = 0;
01421             glImage[gl_index*4+3] = 255;
01422         }
01423 
01424         gl_index++;
01425     }
01426 }
01427 
01428 
01429 //-----------------------------------------------------------------------
01430 //-- Produces the intermediate image for orthogonal projection        ---
01431 //-----------------------------------------------------------------------
01432 void vu111211A::makeIntermediateImageOrtho()
01433 {
01434     int from = 0;
01435     int to = 0;
01436 
01437     int pos = 0;
01438     int add = 0;
01439 
01440     int intermediateWidth = maxSize*2;
01441     int intermediateSize = intermediateWidth*intermediateWidth;
01442 
01443     int step = 0;
01444 
01445     double l = computeZNormalizedViewingVectorLength();
01446 
01447     resetIntermediateImage(intermediateSize);
01448 
01449     if (direction == 0) { to = volumeDepth-1; step = +1; }
01450     else { from = volumeDepth-1; step = -1; }
01451 
01452     int runIndex = 0;
01453     byte state1 = 0;
01454     byte state2 = 0;
01455     byte count1 = 0;
01456     byte count2 = 0;
01457     int voxelIndex = 0;
01458     int volumeIndex = 0;
01459 
01460     //for every plane in the runlength encoded volume...
01461     for (int i = from; i != to + step; i += step) {
01462 
01463         runIndex = 0;
01464         state1 = 0;
01465         state2 = 0;
01466         count1 = 0;
01467         count2 = 0;
01468         voxelIndex = 0;
01469         volumeIndex = 0;
01470 
01471         //index of the voxels of each of the two voxel scanlines
01472         int voxelPtr1 = 0;
01473         int voxelPtr2 = 0;
01474 
01475         //index of the current runs
01476         int runPtr1 = 0;
01477         int runPtr2 = 0;
01478 
01479         //compute the rounded total shearing (rounding dim1_pos and dim2_pos)
01480         int shearX = (int)(curr[i].dim1_pos);
01481         int shearY = (int)(curr[i].dim2_pos);
01482 
01483         //compute the bilinear interpolation-coefficients
01484         float s1 = 1 - (curr[i].dim1_pos - shearX);
01485         float s2 = 1 - (curr[i].dim2_pos - shearY);
01486 
01487         float w1 = s2-s1*s2;
01488         float w2 = s1*s2;
01489         float w3 = 1-s1-s2+s1*s2;
01490         float w4 = s1-s1*s2;
01491 
01492         vuVector normal;
01493 
01494         //if current slice is not empty
01495         if (curr[i].size > 0) {
01496 
01497             // pre-classify the current slice
01498             preClassification(curr[i],l);
01499 
01500             // voxel-pointer 2 is already in right position; 
01501             // find value for voxel-pointer 1
01502             initRunPtr1(curr[i], volumeIndex, voxelPtr1, runPtr1, state1);
01503 
01504             volumeIndex = 0;
01505 
01506             // position of the pixel that is affected by the
01507             // four surrounding voxels
01508             pos = intermediateWidth * shearY + shearX;
01509             add = 0;
01510 
01511             //trace through the slice
01512             while((dword)voxelPtr2 < curr[i].size) {
01513 
01514                 //if neccessary, skip zero-length runs
01515                 skipZeroLengthRuns(curr[i], runPtr1, count1, state1);
01516                 skipZeroLengthRuns(curr[i], runPtr2, count2, state2);
01517 
01518                 //if inside a transparent run, skip voxels
01519                 if ((state1 == 0) && (curr[i].runlength[runPtr1]-count1 > 1) && 
01520                     (state2 == 0) && (curr[i].runlength[runPtr2]-count2 > 1)) {
01521 
01522                     volumeIndex += skipTransparentVoxels(curr[i], 
01523                                                          runPtr1, 
01524                                                          count1, 
01525                                                          runPtr2, 
01526                                                          count2, 
01527                                                          add, 
01528                                                          intermediateWidth,
01529                                                          pos);
01530 
01531                 } else {
01532 
01533                     // only check this pixel, if it still is transparent
01534                     if (intermediate[pos].trans > MIN_TRANSPARENCY_MAKING_INTERMEDIATE) {
01535 
01536                         interpolatePixel(curr[i], intermediate[pos],
01537                                          w1, w2, w3, w4, volumeIndex,
01538                                          voxelPtr1, runPtr1, count1, state1,
01539                                          voxelPtr2, runPtr2, count2, state2);
01540 
01541                     } else {
01542 
01543                         skipOpaquePixel(state1,voxelPtr1,state2,voxelPtr2);
01544 
01545                     }
01546 
01547                     //one step foreward through the volume
01548                     volumeNext(curr[i], volumeIndex, add, intermediateWidth, pos, 
01549                                runPtr1, count1, state1, 
01550                                runPtr2, count2, state2);
01551 
01552                 }
01553             }
01554         }
01555     }
01556 }
01557 
01558 //-----------------------------------------------------------------------
01559 //-- Convert format of intermediate image to OpenGL texture format    ---
01560 //-----------------------------------------------------------------------
01561 void vu111211A::changeFormatIntermediate() {
01562     int w = (int)(si*volumeDepth);
01563     if (w < 0) w *= -1;
01564     w += volumeWidth;
01565 
01566     int h = (int)(sj*volumeDepth);
01567     if (h < 0) h *= -1;
01568     h += volumeHeight;
01569 
01570     int diffX = (int)((glImageWidth-w)/2.0);
01571     int diffY = (int)((glImageHeight-h)/2.0);
01572 
01573     for (word i = 0; i < glImageWidth; ++i) {
01574         for (word j = 0; j < glImageHeight; ++j) {
01575             dword pos_gl = (i + j*glImageWidth)*4;
01576             if ((i-diffX >= 0) && (j-diffY >= 0) && (i-diffX < w) && (j-diffY < h)) {
01577 
01578                 dword pos_int = (i-diffX) + (j-diffY)*maxSize*2;
01579 
01580                 glImage[pos_gl+0] = (GLubyte)(intermediate[pos_int].red*255);
01581                 glImage[pos_gl+1] = (GLubyte)(intermediate[pos_int].green*255);
01582                 glImage[pos_gl+2] = (GLubyte)(intermediate[pos_int].blue*255);
01583                 glImage[pos_gl+3] = 255;
01584 
01585             } else {
01586                 glImage[pos_gl+0] = (GLubyte)0;
01587                 glImage[pos_gl+1] = (GLubyte)0;
01588                 glImage[pos_gl+2] = (GLubyte)0;
01589                 glImage[pos_gl+3] = 255;
01590             }
01591         }
01592     }
01593 }
01594 
01595 //-----------------------------------------------------------------------
01596 //-- Warp and draw image (warp is performed by OpenGL)                ---
01597 //-----------------------------------------------------------------------
01598 void vu111211A::drawWarpOpenGL() {
01599     //the positions of the edges of the OpenGL-polygon
01600     float pos[4][2];
01601 
01602     //the 2D warping matrix
01603     float m00 = viewMatrix[0][0];
01604     float m01 = viewMatrix[0][1];
01605     float m10 = viewMatrix[1][0];
01606     float m11 = viewMatrix[1][1];
01607 
01608     //compute the warped positions of the edges of the OpenGL-polygon
01609     pos[0][XDIR] = 0 * m00 + 0 * m01;
01610     pos[0][YDIR] = 0 * m10 + 0 * m11;
01611 
01612     pos[1][XDIR] = 1 * m00 + 0 * m01;
01613     pos[1][YDIR] = 1 * m10 + 0 * m11;
01614 
01615     pos[2][XDIR] = 1 * m00 + 1 * m01;
01616     pos[2][YDIR] = 1 * m10 + 1 * m11;
01617 
01618     pos[3][XDIR] = 0 * m00 + 1 * m01;
01619     pos[3][YDIR] = 0 * m10 + 1 * m11;
01620 
01621     //compute the translation of the OpenGL-plane to place 
01622     //it in the center of the screen
01623     float t[2];
01624     t[XDIR] = 0.5 * m00 + 0.5 * m01;
01625     t[YDIR] = 0.5 * m10 + 0.5 * m11;
01626 
01627     //clear screen, enable texturing...
01628     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
01629     glEnable(GL_TEXTURE_2D);
01630     glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
01631 
01632     //set parameters for texture mapping
01633     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
01634     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
01635     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
01636     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
01637 
01638     //create new texture
01639     glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, glImageWidth,
01640                  glImageHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, glImage);
01641 
01642     //work with projection matrix
01643     glMatrixMode(GL_PROJECTION);
01644     glLoadIdentity();
01645 
01646     // the size of the OpenGL plane with the texture upon it
01647     float planeSize = 40*zoomValue*zoomValue;  
01648 
01649     glOrtho(-10,10,-10,10,-1000,1000);
01650         
01651     glMatrixMode(GL_MODELVIEW);
01652     glLoadIdentity();
01653 
01654     glBegin(GL_QUADS);
01655 
01656     glTexCoord2f(0, 0); 
01657     glVertex3f(planeSize*(pos[0][XDIR]-t[XDIR]), planeSize*(pos[0][YDIR]-t[YDIR]), 0.0);
01658     glTexCoord2f(1, 0); 
01659     glVertex3f(planeSize*(pos[1][XDIR]-t[XDIR]), planeSize*(pos[1][YDIR]-t[YDIR]), 0.0);
01660     glTexCoord2f(1, 1); 
01661     glVertex3f(planeSize*(pos[2][XDIR]-t[XDIR]), planeSize*(pos[2][YDIR]-t[YDIR]), 0.0);
01662     glTexCoord2f(0, 1); 
01663     glVertex3f(planeSize*(pos[3][XDIR]-t[XDIR]), planeSize*(pos[3][YDIR]-t[YDIR]), 0.0);
01664 
01665     glEnd();
01666 
01667     glFlush();
01668 
01669     glDisable(GL_TEXTURE_2D);
01670 }
01671 
01672 //-----------------------------------------------------------------------
01673 //-- get the value of the voxel at the given coordinates              ---
01674 //-----------------------------------------------------------------------
01675 byte vu111211A::getSample(dword x, dword y, dword z)
01676 {
01677     dword index = x + y * m_Dim1Size + z * m_Dim1Size * m_Dim2Size;
01678     byte result = 0;
01679     if (index < m_Dim1Size*m_Dim2Size*m_Dim3Size) {
01680         result = m_Data[index];
01681     } 
01682     return result;
01683 }
01684 
01685 //-----------------------------------------------------------------------
01686 //-- Produces the intermediate image for orthogonal projection using  ---
01687 //-- Fast Classification                                              ---
01688 //-----------------------------------------------------------------------
01689 void vu111211A::makeIntermediateImageOrthoFastClass()
01690 {
01691     dword from = 0;
01692     dword to = 0;
01693 
01694     dword pos = 0;
01695 
01696     dword row = 0;
01697     dword col = 0;
01698 
01699     dword intermediateWidth = maxSize*2;
01700     dword intermediateSize = intermediateWidth*intermediateWidth;
01701 
01702     int step = 0;
01703 
01704     double l = computeZNormalizedViewingVectorLength();
01705 
01706     //initialize intermediate image
01707     for (dword i = 0; i < intermediateSize; ++i) {
01708         intermediate[i].trans = 1;
01709         intermediate[i].red = 0;
01710         intermediate[i].green = 0;
01711         intermediate[i].blue = 0;
01712         intermediate[i].offset = 0;
01713     }
01714 
01715     //Select start-, end-position and step-size
01716     if (direction == 0) { to = volumeDepth-1; step = +1; }
01717     else { from = volumeDepth-1; step = -1; }
01718 
01719     //for every plane in the runlength encoded volume...
01720     for (dword i = from; i != to + step; i += step) {
01721 
01722         //precompute the rounded total shearing (rounding dim1_pos and dim2_pos)
01723         dword shearX = (dword)(curr[i].dim1_pos);
01724         dword shearY = (dword)(curr[i].dim2_pos);
01725 
01726         //precompute the bilinear interpolation-coefficients
01727         double s1 = 1 - (curr[i].dim1_pos - shearX);
01728         double s2 = 1 - (curr[i].dim2_pos - shearY);
01729 
01730         double w[4];
01731         w[0] = s2-s1*s2;
01732         w[1] = s1*s2;
01733         w[2] = 1-s1-s2+s1*s2;
01734         w[3] = s1-s1*s2;
01735 
01736         //position in intermediate image
01737         pos = intermediateWidth * shearY + shearX;
01738 
01739         switch(mainViewingDir) {
01740             case XDIR: {
01741                 row = m_Dim2Size;
01742                 col = m_Dim3Size;
01743                 break;
01744             }
01745             case YDIR: {
01746                 row = m_Dim3Size;
01747                 col = m_Dim1Size;
01748                 break;
01749             }
01750             case ZDIR: {
01751                 row = m_Dim1Size;
01752                 col = m_Dim2Size;
01753                 break;
01754             }
01755         }
01756 
01757         //The coordinates of all 4 sample points
01758         dword xVol[4] = {0, 0, 0, 0};
01759         dword yVol[4] = {0, 0, 0, 0};
01760         dword zVol[4] = {0, 0, 0, 0};
01761 
01762         dword j = 0;
01763 
01764         while (j < col) {
01765 
01766             dword k = 0;
01767 
01768             while (k < row) {
01769 
01770                 switch(mainViewingDir) {
01771                     case XDIR: {
01772                         xVol[2] = i;
01773                         yVol[2] = k;
01774                         zVol[2] = j;
01775 
01776                         xVol[3] = i;
01777                         yVol[3] = k + 1;
01778                         zVol[3] = j;
01779 
01780                         xVol[0] = i;
01781                         yVol[0] = k;
01782                         zVol[0] = j + 1;
01783 
01784                         xVol[1] = i;
01785                         yVol[1] = k + 1;
01786                         zVol[1] = j + 1;
01787 
01788                         break;
01789                     }
01790                     case YDIR: {
01791                         xVol[2] = j;
01792                         yVol[2] = i;
01793                         zVol[2] = k;
01794 
01795                         xVol[3] = j;
01796                         yVol[3] = i;
01797                         zVol[3] = k + 1;
01798 
01799                         xVol[0] = j + 1;
01800                         yVol[0] = i;
01801                         zVol[0] = k;
01802 
01803                         xVol[1] = j + 1;
01804                         yVol[1] = i;
01805                         zVol[1] = k + 1;
01806 
01807                         break;
01808                     }
01809                     case ZDIR: {
01810                         xVol[2] = k;
01811                         yVol[2] = j;
01812                         zVol[2] = i;
01813 
01814                         xVol[3] = k + 1;
01815                         yVol[3] = j;
01816                         zVol[3] = i;
01817 
01818                         xVol[0] = k;
01819                         yVol[0] = j + 1;
01820                         zVol[0] = i;
01821 
01822                         xVol[1] = k + 1;
01823                         yVol[1] = j + 1;
01824                         zVol[1] = i;
01825 
01826                         break;
01827                     }
01828                 }
01829 
01830 
01831                 if (intermediate[pos].trans > MIN_TRANSPARENCY_MAKING_INTERMEDIATE) {
01832 
01833                     float red     = 0;
01834                     float green   = 0;
01835                     float blue    = 0;
01836                     float opacity = 0;
01837 
01838                     //the normal vector
01839                     vuVector normal;
01840                     normal[XDIR] = 0;
01841                     normal[YDIR] = 0;
01842                     normal[ZDIR] = 0;
01843 
01844                     for (int m = 0; m < 4; ++m) {
01845                         //get the value of the voxel
01846                         byte sample = getSample(xVol[m],yVol[m],zVol[m]);
01847 
01848                         //get the normal vector of the voxel
01849                         vuVector n;
01850                         dword indexNormal = (xVol[m] + yVol[m] * m_Dim1Size + 
01851                                              zVol[m] * m_Dim1Size * m_Dim2Size) * 3;
01852                         n[XDIR] = m_Normals[indexNormal+0];
01853                         n[YDIR] = m_Normals[indexNormal+1];
01854                         n[ZDIR] = m_Normals[indexNormal+2];
01855                         n.makeUnit();
01856 
01857                         //weigh the normal vector
01858                         normal[XDIR] += w[m] * n[XDIR];
01859                         normal[YDIR] += w[m] * n[YDIR];
01860                         normal[ZDIR] += w[m] * n[ZDIR];
01861 
01862                         red     += w[m] * m_TFunc[sample][0];
01863                         green   += w[m] * m_TFunc[sample][1];
01864                         blue    += w[m] * m_TFunc[sample][2];
01865                         opacity += w[m] * m_TFunc[sample][3];
01866                     }
01867 
01868                     opacity = 1 - pow((1-opacity),(float)l);
01869 
01870                     intermediate[pos].offset = 1;
01871 
01872                     //compute shading
01873                     normal.makeUnit();
01874                     float diffuse = normal.dot(m_View);
01875 
01876                     if (diffuse > 0) {
01877 
01878                         double fac = intermediate[pos].trans * opacity * (0.1 + 0.9 * diffuse);
01879 
01880                         intermediate[pos].red   += red   * fac;
01881                         intermediate[pos].green += green * fac;
01882                         intermediate[pos].blue  += blue  * fac;
01883                         intermediate[pos].trans *= (1-opacity);
01884 
01885                         if (specular == 1) {
01886                             float spec = pow(diffuse,150);
01887                             intermediate[pos].red     += spec * opacity;
01888                             intermediate[pos].green   += spec * opacity;
01889                             intermediate[pos].blue    += spec * opacity;
01890                             intermediate[pos].trans   -= spec * opacity;
01891 
01892                             if (intermediate[pos].red     > 1.0) intermediate[pos].red     = 1.0;
01893                             if (intermediate[pos].green   > 1.0) intermediate[pos].green   = 1.0;
01894                             if (intermediate[pos].blue    > 1.0) intermediate[pos].blue    = 1.0;
01895                             if (intermediate[pos].trans   < 0.0) intermediate[pos].trans   = 0.0;
01896                         }
01897                     }
01898                 }
01899 
01900                 int skip1 = fastClassificationObj->skip(xVol[1],yVol[1],zVol[1],mainViewingDir);
01901                 int skip2 = fastClassificationObj->skip(xVol[3],yVol[3],zVol[3],mainViewingDir);
01902 
01903                 if (skip2 < skip1) skip1 = skip2;
01904 
01905                 pos += skip1;
01906                 k   += skip1;
01907 
01908             }
01909 
01910             pos += (intermediateWidth - volumeWidth);
01911 
01912             ++j;
01913         }
01914     }
01915 
01916     //-- change format of intermediate Image ------------------
01917 
01918     int w = (int)(si*volumeDepth);
01919     if (w < 0) w *= -1;
01920     w += volumeWidth;
01921 
01922     int h = (int)(sj*volumeDepth);
01923     if (h < 0) h *= -1;
01924     h += volumeHeight;
01925 
01926     int diffX = (int)((glImageWidth-w)/2.0);
01927     int diffY = (int)((glImageHeight-h)/2.0);
01928 
01929     for (word i = 0; i < glImageWidth; ++i) {
01930         for (word j = 0; j < glImageHeight; ++j) {
01931             dword pos_gl = (i + j*glImageWidth)*4;
01932             if ((i-diffX >= 0) && (j-diffY >= 0) && (i-diffX < w) && (j-diffY < h)) {
01933 
01934                 dword pos_int = (i-diffX) + (j-diffY)*maxSize*2;
01935 
01936                 glImage[pos_gl+0] = (GLubyte)(intermediate[pos_int].red*255);
01937                 glImage[pos_gl+1] = (GLubyte)(intermediate[pos_int].green*255);
01938                 glImage[pos_gl+2] = (GLubyte)(intermediate[pos_int].blue*255);
01939                 glImage[pos_gl+3] = 255;
01940 
01941             } else {
01942                 glImage[pos_gl+0] = (GLubyte)0;
01943                 glImage[pos_gl+1] = (GLubyte)0;
01944                 glImage[pos_gl+2] = (GLubyte)0;
01945                 glImage[pos_gl+3] = 255;
01946             }
01947         }
01948     }
01949 }
01950 
01951 int vu111211A::getMaxSize() {
01952     return maxSize;
01953 }
01954 
01955 float vu111211A::getMinEyeDistance() {
01956     return maxSize * 1.0f;
01957 }
01958 
01959 void vu111211A::setEyeDistance(float distance) {
01960     eye_distance = distance;
01961 }
01962 
01963 //-----------------------------------------------------------------------
01964 //-- Computes the position of the eye point in object space           ---
01965 //-----------------------------------------------------------------------
01966 void vu111211A::computeEyePoint() 
01967 {
01968     vuVector eye_w;
01969     eye_w[XDIR] = 0.0f;
01970     eye_w[YDIR] = 0.0f;
01971     eye_w[ZDIR] = - eye_distance;
01972     eye_w[WDIR] = 1.0f;
01973 
01974     vuMatrix invWorldMatrix = worldMatrix.inverse();
01975 
01976     eye_o = eye_w * invWorldMatrix * permMatrix;
01977 }
01978 
01979 //-----------------------------------------------------------------------
01980 //-- Computes the Shear Factors si, sj and the Warp Factors ti, tj    ---
01981 //-- for perspective projection.                                      ---
01982 //-----------------------------------------------------------------------
01983 void vu111211A::computeShearPerspective()
01984 {
01985     //init the shearMatrix with the identity
01986     shearMatrix.makeIdentity();
01987 
01988     //change the three values in the matrix to get the final shear matrix
01989     shearMatrix[0][2] = - eye_o[XDIR] / eye_o[ZDIR];
01990     shearMatrix[1][2] = - eye_o[YDIR] / eye_o[ZDIR];
01991     shearMatrix[3][2] = - eye_o[WDIR] / eye_o[ZDIR];
01992 }
01993 
01994 
01995 //-----------------------------------------------------------------------
01996 //-- Shear the planes of the volume for perspective projection        ---
01997 //-----------------------------------------------------------------------
01998 void vu111211A::computePlanePositionScale(float &pos_i, 
01999                                           float &pos_j,
02000                                           float &scale,
02001                                           int plane)
02002 {
02003     int halfX = volumeWidth  / 2;
02004     int halfY = volumeHeight / 2;
02005 
02006     vuVector lu = vuVector(-halfX, -halfY, (float)plane, 1.0f);  // left upper point of slice
02007     vuVector ru = vuVector( halfX, -halfY, (float)plane, 1.0f);   // right upper point of slice
02008     vuVector ll = vuVector(-halfX,  halfY, (float)plane, 1.0f);   // left lower point of slice
02009     vuVector rl = vuVector( halfX,  halfY, (float)plane, 1.0f);    // right lower point of slice
02010 
02011     vuVector lu_t = lu * shearMatrix;
02012     vuVector ru_t = ru * shearMatrix;
02013     vuVector ll_t = ll * shearMatrix;
02014     vuVector rl_t = rl * shearMatrix;
02015 
02016     // homogenous coordinates! -> divide by 4th coordinate
02017     for (int i = 0; i < 4; ++i) {
02018         lu_t[i] = lu_t[i] / lu_t[WDIR];
02019         ru_t[i] = ru_t[i] / ru_t[WDIR];
02020         ll_t[i] = ll_t[i] / ll_t[WDIR];
02021         rl_t[i] = rl_t[i] / rl_t[WDIR];
02022     }
02023 
02024     if (scale_total < 0) {
02025         scale_total = (lu_t[XDIR] - rl_t[XDIR]) / (lu[XDIR] - rl[XDIR]);
02026     }
02027 
02028     scale = ((lu_t[XDIR] - ru_t[XDIR]) / (lu[XDIR] - ru[XDIR])) / scale_total;
02029 
02030     pos_i = lu_t[XDIR] / scale_total + maxSize;
02031     pos_j = lu_t[YDIR] / scale_total + maxSize;
02032 }
02033 
02034 //-----------------------------------------------------------------------
02035 //-- Shear the planes of the volume for perspective projection        ---
02036 //-----------------------------------------------------------------------
02037 void vu111211A::shearPerspective()
02038 {
02039     // normalize the viewing vector
02040     m_View.makeUnit();
02041 
02042     // in which direction and from where shall the volume be traversed:
02043     int step = 0;
02044     int start = 0;
02045     int end = 0;
02046 
02047 // !!!!
02048 //    dword intermediateWidth = maxSize*2;
02049 //    dword intermediateSize = intermediateWidth*intermediateWidth;
02050 
02051 //    float shearPos_i = intermediateWidth / 2;
02052 //    float shearPos_j = intermediateWidth / 2;
02053 
02054     // select the used dataset...
02055     switch(mainViewingDir) {
02056         case XDIR: {
02057             curr = dataX;
02058             break;
02059         }
02060         case YDIR: {
02061             curr = dataY;
02062             break;
02063         }
02064         case ZDIR: {
02065             curr = dataZ;
02066             break;
02067         }
02068     }
02069 
02070     if (m_View[mainViewingDir] >= 0) {
02071         direction = 0;
02072         step = 1;
02073         start = 0;
02074         end = volumeDepth;
02075     } else {
02076         direction = 1;
02077         step = -1;
02078         start = volumeDepth-1;
02079         end = -1;
02080     }
02081 
02082     int z = 0;
02083 
02084     if (direction == 0) {
02085         z = -(volumeDepth / 2);
02086     } else {
02087         z = volumeDepth / 2;
02088     }
02089 
02090     // needs to be assigned a negative value -> computePlanePositionScale
02091     // computes the total scale-value
02092     scale_total = -100.0f;
02093 
02094     for (int i = start; i != end; i += step) {
02095         computePlanePositionScale(curr[i].dim1_pos,curr[i].dim2_pos,curr[i].scale,z);
02096         if (direction == 0) z += 1; else z -= 1;
02097     }
02098 }
02099 
02100 
02101 //-----------------------------------------------------------------------
02102 //-- Compute the warp- and the inverse warp-matrix for                ---
02103 //-- perspective viewing (scaling of the most front plane is          ---
02104 //-- included in the computation)                                     ---
02105 //-----------------------------------------------------------------------
02106 void vu111211A::computeWarpPerspective()
02107 {
02108     vuMatrix scaleMatrix;
02109     scaleMatrix.makeIdentity();
02110     scaleMatrix.makeScale(scale_total,scale_total,1.0f);
02111 
02112     warpMatrix = scaleMatrix * shearMatrix.inverse() * permMatrix.inverse() * viewMatrix;
02113 
02114     invWarpMatrix = warpMatrix.inverse();
02115 }
02116 
02117 //-----------------------------------------------------------------------
02118 //-- inverse perspective warping of a single point                    ---
02119 //-----------------------------------------------------------------------
02120 void vu111211A::warpPerspectiveInv(float x, 
02121                                    float y, 
02122                                    float &x_result, 
02123                                    float &y_result)
02124 {
02125     vuVector v;
02126     v[XDIR] = x;
02127     v[YDIR] = y;
02128     v[ZDIR] = 0.0f;
02129     v[WDIR] = 1.0f;
02130 
02131     // apply the inverse warp to the given coordinates
02132     vuVector result = v * invWarpMatrix;
02133 
02134     // the result must be divided by the 4th coordinate (homogenous coordinates!)
02135     x_result = result[XDIR] / result[WDIR];
02136     y_result = result[YDIR] / result[WDIR];
02137 }
02138 
02139 //-----------------------------------------------------------------------
02140 //-- perspective warping                                              ---
02141 //-----------------------------------------------------------------------
02142 void vu111211A::warpPerspective()
02143 {
02144     int glSize = glImageWidth * glImageHeight;
02145     int glHalf = glImageWidth / 2;
02146 
02147     int intermediateWidth = maxSize * 2;
02149 
02150     float int_x = 0.0f;
02151     float int_y = 0.0f;
02152 
02153     float gl_x = 0.0f;
02154     float gl_y = 0.0f;
02155 
02156     float maxSizef = (float)maxSize;
02157 
02158     int int_index = 0;
02159     int gl_index  = 0;
02160     while(gl_index < glSize) {
02161 
02162         gl_x = (gl_index % glImageWidth) - glHalf;
02163         gl_y = (gl_index / glImageWidth) - glHalf;
02164 
02165         warpPerspectiveInv(gl_x, gl_y, int_x, int_y);
02166 
02167         float d1 = (int_x - floor(int_x));
02168         float d2 = (int_y - floor(int_y));
02169 
02170         float w1 = 1-d1-d2+d1*d2;
02171         float w2 = d1-d1*d2;
02172         float w3 = d2-d1*d2;
02173         float w4 = d1*d2;
02174 
02175         int_index = (int)(floor(int_x + maxSizef) + floor(int_y + maxSizef) * intermediateWidth);
02176 
02177         if ((int_x > -maxSizef) && (int_x < maxSizef-1) &&
02178             (int_y > -maxSizef) && (int_y < maxSizef-1)) {
02179 
02180             float red = intermediate[int_index].red*w1 + 
02181                     intermediate[int_index+1].red*w2 +
02182                     intermediate[int_index+intermediateWidth].red*w3 +
02183                     intermediate[int_index+intermediateWidth+1].red*w4;
02184             float green = intermediate[int_index].green*w1 + 
02185                     intermediate[int_index+1].green*w2 +
02186                     intermediate[int_index+intermediateWidth].green*w3 +
02187                     intermediate[int_index+intermediateWidth+1].green*w4;
02188             float blue = intermediate[int_index].blue*w1 + 
02189                     intermediate[int_index+1].blue*w2 +
02190                     intermediate[int_index+intermediateWidth].blue*w3 +
02191                     intermediate[int_index+intermediateWidth+1].blue*w4;
02192 
02193             glImage[gl_index*4+0] = (GLubyte)(red*255);
02194             glImage[gl_index*4+1] = (GLubyte)(green*255);
02195             glImage[gl_index*4+2] = (GLubyte)(blue*255);
02196             glImage[gl_index*4+3] = 255;
02197         } else {
02198             glImage[gl_index*4+0] = 0;
02199             glImage[gl_index*4+1] = 0;
02200             glImage[gl_index*4+2] = 0;
02201             glImage[gl_index*4+3] = 255;
02202         }
02203 
02204         gl_index++;
02205     }
02206 }
02207 
02208 //-----------------------------------------------------------------------
02209 //-- Produces the intermediate image for perspective projection       ---
02210 //-----------------------------------------------------------------------
02211 void vu111211A::makeIntermediateImagePerspective()
02212 {
02213     dword from = 0;
02214     dword to = 0;
02215 
02216     dword pos = 0;
02217 
02218     dword intermediateWidth = maxSize*2;
02219     dword intermediateSize = intermediateWidth*intermediateWidth;
02220 
02221     int step = 0;
02222 
02223     double l = computeZNormalizedViewingVectorLength();
02224 
02225     //initialize intermediate image
02226     for (dword i = 0; i < intermediateSize; ++i) {
02227         intermediate[i].trans = 1;
02228         intermediate[i].red = 0;
02229         intermediate[i].green = 0;
02230         intermediate[i].blue = 0;
02231         intermediate[i].offset = 0;
02232     }
02233 
02234     if (direction == 0) { to = volumeDepth-1; step = +1; }
02235     else { from = volumeDepth-1; step = -1; }
02236 
02237     dword runIndex = 0;
02238     byte state1 = 0;
02239     byte state2 = 0;
02240     byte count1 = 0;
02241     byte count2 = 0;
02242     dword voxelIndex = 0;
02243     dword indexRunlength = 0;
02244 
02245     //for every plane in the runlength encoded volume...
02246     for (dword i = from; i != to + step; i += step) {
02247 //    for (dword i = 10; i == 10; i += step) {
02248 
02249         runIndex = 0;
02250         state1 = 0;
02251         state2 = 0;
02252         count1 = 0;
02253         count2 = 0;
02254         voxelIndex = 0;
02255         indexRunlength = 0;
02256 
02257         //index of the voxels of each of the two voxel scanlines
02258         dword voxelPtr1 = 0;
02259         dword voxelPtr2 = 0;
02260 
02261         //index of the current runs
02262         dword runPtr1 = 0;
02263         dword runPtr2 = 0;
02264 
02265         //precompute the rounded total shearing (rounding dim1_pos and dim2_pos)
02266 // !!!! dword shearX = (dword)(curr[i].dim1_pos);
02267 //      dword shearY = (dword)(curr[i].dim2_pos);
02268 
02269         vuVector normal;
02270 
02271         //value of sampled point resample = w1*p1 + w2*p2 + w3*p3 + w4*p4
02272         // wi precomputed weights for every voxel;
02273         // pi points:
02274         //               p3    p4
02275         //                  sp
02276         //               p1    p2
02277 
02278         if (curr[i].size > 0) {
02279 
02280 //---------------------------------------
02281             // do the shading and the opacity correction
02282             for (dword k = 0; k < curr[i].size; ++k) {
02283                 curr[i].voxel[k].normal.makeUnit();
02284 
02285                 float diffuse = m_View.dot(curr[i].voxel[k].normal);
02286 
02287                 float spec = 0;
02288 
02289                 if (diffuse <= 0) { // <= 0
02290                     curr[i].voxel[k].red_shaded = 0.0;
02291                     curr[i].voxel[k].green_shaded = 0.0;
02292                     curr[i].voxel[k].blue_shaded = 0.0;
02293                     curr[i].voxel[k].opacity_corr = 0.0;
02294                 } else {
02295                     curr[i].voxel[k].opacity_corr = 1 - pow((1-curr[i].voxel[k].opacity),(float)l);
02296 
02297                     float fac = curr[i].voxel[k].opacity_corr * (0.1 + 0.9 * diffuse);
02298 
02299                     curr[i].voxel[k].red_shaded = curr[i].voxel[k].red * fac;
02300                     curr[i].voxel[k].green_shaded = curr[i].voxel[k].green * fac;
02301                     curr[i].voxel[k].blue_shaded = curr[i].voxel[k].blue * fac;
02302 
02303                     // if neccessary, add specular light
02304                     if (specular == 1) {
02305                         spec = pow(diffuse,150);
02306 
02307                         curr[i].voxel[k].red_shaded   += spec;
02308                         curr[i].voxel[k].green_shaded += spec;
02309                         curr[i].voxel[k].blue_shaded  += spec;
02310                         curr[i].voxel[k].opacity_corr += spec;
02311 
02312                         if (curr[i].voxel[k].red_shaded > 1.0) curr[i].voxel[k].red_shaded = 1.0;
02313                         if (curr[i].voxel[k].green_shaded > 1.0) curr[i].voxel[k].green_shaded = 1.0; 
02314                         if (curr[i].voxel[k].blue_shaded > 1.0) curr[i].voxel[k].blue_shaded = 1.0;
02315                         if (curr[i].voxel[k].opacity_corr > 1.0) curr[i].voxel[k].opacity_corr = 1.0;
02316                     }
02317                 }
02318             }
02319 
02320 //---------------------------------------
02321 
02322             //voxel-pointer 2 is already in right position; find value for voxel-pointer 1
02323             do {
02324                 indexRunlength += curr[i].runlength[runPtr1];
02325                 if (state1 == 0) {
02326                     state1 = 1;
02327                 } else {
02328                     state1 = 0;
02329                     voxelPtr1 += curr[i].runlength[runPtr1];
02330                 }
02331                 ++runPtr1;
02332             } while (indexRunlength < volumeWidth);
02333 
02334             indexRunlength = 0;
02335 
02336             dword skip;
02337 
02338             //trace through the plane
02339             while(voxelPtr2 < curr[i].size) {
02340 
02341                 //if neccessary, skip zero-length runs
02342                 while (curr[i].runlength[runPtr1] == 0) { 
02343                     ++runPtr1;
02344                     count1 = 0;
02345                     if (state1 == 0) state1 = 1; else state1 = 0;
02346                 }
02347                 while (curr[i].runlength[runPtr2] == 0) {
02348                     ++runPtr2;
02349                     count2 = 0;
02350                     if (state2 == 0) state2 = 1; else state2 = 0;
02351                 }
02352 
02353                 //if inside a transparent run, skip voxels
02354                 if ((state1 == 0) && (curr[i].runlength[runPtr1]-count1 > 1) && 
02355                     (state2 == 0) && (curr[i].runlength[runPtr2]-count2 > 1)) {
02356 
02357                     if (curr[i].runlength[runPtr1]-count1 < curr[i].runlength[runPtr2]-count2) {
02358                         skip = curr[i].runlength[runPtr1]-count1-1;
02359                     } else {
02360                         skip = curr[i].runlength[runPtr2]-count2-1;
02361                     }
02362 
02363                     indexRunlength += skip;
02364 
02365                     count1 += skip;
02366                     count2 += skip;
02367                 } else {
02368 
02369                     word xpos = indexRunlength % volumeWidth;
02370                     word ypos = indexRunlength / volumeWidth;
02371 
02372                     // positions of the left upper and right lower sample for the intermediate image
02373                     float sx1 = curr[i].dim1_pos + (float)xpos * curr[i].scale;
02374                     float sy1 = curr[i].dim2_pos + (float)ypos * curr[i].scale;
02375                     float sx4 = curr[i].dim1_pos + (float)(xpos+1) * curr[i].scale;
02376                     float sy4 = curr[i].dim2_pos + (float)(ypos+1) * curr[i].scale;
02377 
02378                     pos = intermediateWidth * ((dword)sy4) + ((dword)sx4);
02379 
02380                     // only check this pixel, if it still is transparent
02381                     if ((intermediate[pos].trans > MIN_TRANSPARENCY_MAKING_INTERMEDIATE) &&
02382                         (floor(sx1) != floor(sx4)) && (floor(sy1) != floor(sy4))) {
02383 
02384                         double red = 0;
02385                         double green = 0;
02386                         double blue = 0;
02387                         double opacity = 0;
02388 
02389                         //precompute the bilinear interpolation-coefficients
02390                         double s1 = 1 - (sx1 - (int)sx1);
02391                         double s2 = 1 - (sy1 - (int)sy1);
02392 
02393                         s1 /= curr[i].scale;
02394                         s2 /= curr[i].scale;
02395 
02396                         double w1 = s2-s1*s2;
02397                         double w2 = s1*s2;
02398                         double w3 = 1-s1-s2+s1*s2;
02399                         double w4 = s1-s1*s2;
02400 
02401                         if (voxelPtr1 < curr[i].size) {
02402                             if (state1 == 1) {
02403 
02404                                 red     += w1*curr[i].voxel[voxelPtr1].red_shaded;
02405                                 green   += w1*curr[i].voxel[voxelPtr1].green_shaded;
02406                                 blue    += w1*curr[i].voxel[voxelPtr1].blue_shaded;
02407                                 opacity += w1*curr[i].voxel[voxelPtr1].opacity_corr;
02408 
02409                                 if ((curr[i].runlength[runPtr1]-count1 > 1) &&
02410                                     (voxelPtr1+1 < curr[i].size)) {
02411 
02412                                     red     += w2*curr[i].voxel[voxelPtr1+1].red_shaded;
02413                                     green   += w2*curr[i].voxel[voxelPtr1+1].green_shaded;
02414                                     blue    += w2*curr[i].voxel[voxelPtr1+1].blue_shaded;
02415                                     opacity += w2*curr[i].voxel[voxelPtr1+1].opacity_corr;
02416 
02417                                 }
02418                                 ++voxelPtr1;
02419                             } else {
02420                                 if ((curr[i].runlength[runPtr1]-count1 == 1) &&   //non-transp. run ends after ptr1?
02421                                     (curr[i].runlength[runPtr1+1]!=0) &&          //no zero-length run after current?
02422                                     ((indexRunlength+1) % volumeWidth)!=0) {   //no line break?
02423 
02424                                     red     += w2*curr[i].voxel[voxelPtr1].red_shaded;  // vorher voxelPtr1+1
02425                                     green   += w2*curr[i].voxel[voxelPtr1].green_shaded;
02426                                     blue    += w2*curr[i].voxel[voxelPtr1].blue_shaded;
02427                                     opacity += w2*curr[i].voxel[voxelPtr1].opacity_corr;
02428 
02429                                 }
02430                             }
02431                         }
02432 
02433                         if (voxelPtr2 < curr[i].size) {
02434                             if (state2 == 1) {
02435 
02436                                 red     += w3*curr[i].voxel[voxelPtr2].red_shaded;
02437                                 green   += w3*curr[i].voxel[voxelPtr2].green_shaded;
02438                                 blue    += w3*curr[i].voxel[voxelPtr2].blue_shaded;
02439                                 opacity += w3*curr[i].voxel[voxelPtr2].opacity_corr;
02440 
02441                                 if ((curr[i].runlength[runPtr2]-count2 > 1) &&
02442                                     (voxelPtr2+1 < curr[i].size)) {
02443 
02444                                     red     += w4*curr[i].voxel[voxelPtr2+1].red_shaded;
02445                                     green   += w4*curr[i].voxel[voxelPtr2+1].green_shaded;
02446                                     blue    += w4*curr[i].voxel[voxelPtr2+1].blue_shaded;
02447                                     opacity += w4*curr[i].voxel[voxelPtr2+1].opacity_corr;
02448 
02449                                 }
02450                                 ++voxelPtr2;
02451                             } else {
02452                                 if ((curr[i].runlength[runPtr2]-count2 == 1) &&  //non-transp. run ends after ptr1?
02453                                     (curr[i].runlength[runPtr2+1]!=0) &&         //no zero-length run after current?
02454                                     ((indexRunlength+1) % volumeWidth)!=0) {  //no line break?
02455                                     
02456                                     red     += w4*curr[i].voxel[voxelPtr2].red_shaded;  //vorher voxelPtr2+1
02457                                     green   += w4*curr[i].voxel[voxelPtr2].green_shaded;
02458                                     blue    += w4*curr[i].voxel[voxelPtr2].blue_shaded;
02459                                     opacity += w4*curr[i].voxel[voxelPtr2].opacity_corr;
02460                                 }
02461                             }
02462                         }
02463 
02464 
02465                         if (opacity > 0) {
02466 
02467                             intermediate[pos].offset = 1;
02468 
02469                             double fac = intermediate[pos].trans;
02470 
02471                             intermediate[pos].red   += red * fac;
02472                             intermediate[pos].green += green * fac;
02473                             intermediate[pos].blue  += blue * fac;
02474 
02475                             intermediate[pos].trans *= (1-opacity);
02476                         }
02477 
02478                     } else {
02479 
02480                         if (state1 == 1) {
02481                             ++voxelPtr1;
02482                         }
02483 
02484                         if (state2 == 1) {
02485                             ++voxelPtr2;
02486                         } 
02487 
02488                     }
02489 
02490                     //one step foreward through the volume
02491                     ++indexRunlength;
02492 
02493                     ++count1;
02494                     ++count2;
02495 
02496                     if (count1 >= curr[i].runlength[runPtr1]) {
02497                         ++runPtr1;
02498                         count1 = 0;
02499                         if (state1 == 1) state1 = 0; else state1 = 1;
02500                     }
02501 
02502                     if (count2 >= curr[i].runlength[runPtr2]) {
02503                         ++runPtr2;
02504                         count2 = 0;
02505                         if (state2 == 1) state2 = 0; else state2 = 1;
02506                     }
02507                 }
02508             }
02509         }
02510     }
02511 }
02512 
02513 //-----------------------------------------------------------------------
02514 //-- Draws the perspective intermediate Image using OpenGL --------------
02515 //-----------------------------------------------------------------------
02516 void vu111211A::drawOpenGL() 
02517 {
02518     //clear screen, enable texturing...
02519     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
02520     glEnable(GL_TEXTURE_2D);
02521     glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
02522 
02523     //set parameters for texture mapping
02524     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
02525     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
02526     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
02527     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
02528 
02529     //create new texture
02530     glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, glImageWidth,
02531                  glImageHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, glImage);
02532 
02533     //work with projection matrix
02534     glMatrixMode(GL_PROJECTION);
02535     glLoadIdentity();
02536 
02537     // the size of the OpenGL plane with the texture upon it
02538     float planeSize = 40*zoomValue*zoomValue;  
02539 
02540     glOrtho(-10,10,-10,10,-1000,1000);
02541 
02542     glMatrixMode(GL_MODELVIEW);
02543     glLoadIdentity();
02544 
02545     glBegin(GL_QUADS);
02546 
02547     glTexCoord2f(0, 0); glVertex3f(planeSize*(-1), planeSize*(-1), 0.0);
02548     glTexCoord2f(1, 0); glVertex3f(planeSize*(1), planeSize*(-1), 0.0);
02549     glTexCoord2f(1, 1); glVertex3f(planeSize*(1), planeSize*(1), 0.0);
02550     glTexCoord2f(0, 1); glVertex3f(planeSize*(-1), planeSize*(1), 0.0);
02551 
02552     glEnd();
02553 
02554     glFlush();
02555 
02556     glDisable(GL_TEXTURE_2D);
02557 }
02558 
02559 //-----------------------------------------------------------------------
02560 //-- Initializes OpenGL -------------------------------------------------
02561 //-----------------------------------------------------------------------
02562 void vu111211A::initOpenGL(void)
02563 {
02564     // selecting background color
02565     glClearColor(0.0, 0.0, 0.0, 0.0);
02566 
02567     glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
02568 
02569     glGenTextures(1, &m_GLShearWarp);
02570 
02571     glBindTexture(GL_TEXTURE_2D, m_GLShearWarp);
02572 
02573     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
02574     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
02575     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
02576     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
02577 
02578     glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, glImageWidth,
02579                  glImageHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, glImage);
02580 }
02581 
02582 
02583 //-----------------------------------------------------------------------
02584 //-- Renders the volume -------------------------------------------------
02585 //-----------------------------------------------------------------------
02586 void vu111211A::render() 
02587 {
02588     computeMainViewingDirection();
02589     computePermutationMatrix();
02590     computeViewMatrix();
02591     computeRunlengthSizes();
02592 
02593     switch(viewingMode) {
02594         case VIEWING_MODE_ORTHO: {
02595             computeShearAndWarpFactorsOrtho();
02596             shearOrtho();
02597             computeInvWarpOrtho();
02598             if (fastClassification == 1) {
02599 
02600                 //Save the transfer function in an array
02601                 float *tf = new float[256];
02602                 for (int i = 0; i < 256; ++i) {
02603                     tf[i] = m_TFunc[i][3];
02604                 }
02605                 fastClassificationObj->buildSummedAreaTable(tf,THRESHOLD_RUNLENGTH);
02606                 fastClassificationObj->classify();
02607 
02608                 makeIntermediateImageOrthoFastClass();
02609             } else {
02610                 makeIntermediateImageOrtho();
02611             }
02612             if (orthoWarpOpenGL == 0) {
02613                 warpOrtho();
02614                 drawOpenGL();
02615             } else {
02616                 changeFormatIntermediate();
02617                 drawWarpOpenGL();
02618             }
02619             break;
02620         }
02621         case VIEWING_MODE_PERSPECTIVE: {
02622             computeEyePoint();
02623             computeShearPerspective();
02624             shearPerspective();
02625             computeWarpPerspective();
02626             makeIntermediateImagePerspective();
02627             warpPerspective();
02628             drawOpenGL();
02629             break;
02630         }
02631     }
02632 }
02633 
02634 
02635 
02636 
02637 
02638 
02639 
02640 

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