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

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

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