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