00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00012
00013 #include "octree.h"
00014 #include "math.h"
00015
00016
00017
00018
00019 Octree::Octree(byte *volumeData, int dim1Size, int dim2Size, int dim3Size)
00020 {
00021 width = dim1Size;
00022 height = dim2Size;
00023 depth = dim3Size;
00024 data = volumeData;
00025
00026 step_dim3 = dim1Size * dim2Size;
00027
00028 int max = width;
00029 if (height > max) max = height;
00030 if (depth > max) max = depth;
00031
00032 int treeDepth = (int)(log((double)max)/log(2.0f)) - 1;
00033 if (treeDepth < 1) treeDepth = 1;
00034
00035 root = build(0,0,0,width-1,height-1,depth-1,treeDepth);
00036 }
00037
00038
00039
00040
00041 Octree::~Octree()
00042 {
00043 remove(root);
00044 }
00045
00046
00047
00048
00049 void Octree::remove(Node *tree)
00050 {
00051 for (int i = 0; i < 8; ++i) {
00052 if (tree->children[i] != 0) remove(tree->children[i]);
00053 }
00054
00055 delete tree;
00056 tree = NULL;
00057 }
00058
00059
00060
00061
00062 Node *Octree::build(int minX, int minY, int minZ,
00063 int maxX, int maxY, int maxZ, int count)
00064 {
00065 Node *result;
00066
00067 int partX = 1;
00068 int partY = 1;
00069 int partZ = 1;
00070
00071 if (minX > maxX) minX = maxX;
00072 if (minY > maxY) minY = maxY;
00073 if (minZ > maxZ) minZ = maxZ;
00074
00075 result = new Node();
00076
00077 result->maxX = maxX;
00078 result->maxY = maxY;
00079 result->maxZ = maxZ;
00080 result->minX = minX;
00081 result->minY = minY;
00082 result->minZ = minZ;
00083
00084 if (count == 0) {
00085
00086 for (int i = 0; i < 8; ++i) {
00087 result->children[i] = 0;
00088 }
00089 result->noChildren = 1;
00090
00091 result->min = 255;
00092 result->max = 0;
00093
00094 for (int x = minX; x <= maxX; ++x) {
00095 for (int y = minY; y <= maxY; ++y) {
00096 for (int z = minZ; z <= maxZ; ++z) {
00097 byte voxel = data[x+y*width+z*step_dim3];
00098 if (voxel > result->max) result->max = voxel;
00099 if (voxel < result->min) result->min = voxel;
00100 }
00101 }
00102 }
00103
00104 } else {
00105 if ((minX == maxX) && (minY == maxY) && (minZ == maxZ)) {
00106
00107
00108 for (int i = 0; i < 8; ++i) {
00109 result->children[i] = 0;
00110 }
00111 result->noChildren = 1;
00112
00113 result->min = result->max = data[minX+minY*width+minZ*step_dim3];
00114
00115 } else {
00116
00117 int midX = (minX + maxX) / 2;
00118 int midY = (minY + maxY) / 2;
00119 int midZ = (minZ + maxZ) / 2;
00120
00121 if (minX == maxX) partX = 0;
00122 if (minY == maxY) partY = 0;
00123 if (minZ == maxZ) partZ = 0;
00124
00125 for (int i = 0; i < 8; ++i) result->children[i] = 0;
00126
00127 result->children[0] = build(minX, minY, minZ, midX, midY, midZ, count-1);
00128 if (partX == 1)
00129 result->children[1] = build(midX+1, minY, minZ, maxX, midY, midZ, count-1);
00130 else result->children[1] = 0;
00131 if (partY == 1)
00132 result->children[2] = build(minX, midY+1, minZ, midX, maxY, midZ, count-1);
00133 else result->children[2] = 0;
00134 if ((partX == 1) && (partY == 1))
00135 result->children[3] = build(midX+1, midY+1, minZ, maxX, maxY, midZ, count-1);
00136 else result->children[3] = 0;
00137 if (partZ == 1)
00138 result->children[4] = build(minX, minY, midZ+1, midX, midY, maxZ, count-1);
00139 else result->children[4] = 0;
00140 if ((partX == 1) && (partZ == 1))
00141 result->children[5] = build(midX+1, minY, midZ+1, maxX, midY, maxZ, count-1);
00142 else result->children[5] = 0;
00143 if ((partY == 1) && (partZ == 1))
00144 result->children[6] = build(minX, midY+1, midZ+1, midX, maxY, maxZ, count-1);
00145 else result->children[6] = 0;
00146 if ((partX == 1) && (partY == 1) && (partZ == 1))
00147 result->children[7] = build(midX+1, midY+1, midZ+1, maxX, maxY, maxZ, count-1);
00148 else result->children[7] = 0;
00149
00150 result->noChildren = 0;
00151
00152 result->min = 255;
00153 result->max = 0;
00154
00155 for (int i = 0; i < 8; ++i) {
00156 if (result->children[i] != 0) {
00157 if (result->children[i]->min < result->min) result->min = result->children[i]->min;
00158 if (result->children[i]->max > result->max) result->max = result->children[i]->max;
00159 }
00160 }
00161 }
00162 }
00163
00164 return result;
00165 }
00166
00167
00168
00169
00170 void Octree::classify(SummedAreaTable *sat)
00171 {
00172 printf("Starting classification");
00173 classify(sat,root);
00174 printf("ready.");
00175 }
00176
00177
00178
00179
00180 void Octree::classify(SummedAreaTable *sat, Node *tree)
00181 {
00182 if (sat->transparent(tree->min,tree->max) == 1) {
00183 tree->classification = ALL_TRANSPARENT;
00184 } else {
00185 tree->classification = ALL_NON_TRANSPARENT;
00186 for (int i = 0; i < 8; ++i) {
00187 if (tree->children[i] != 0) {
00188 tree->classification = COMBINATION;
00189 classify(sat,tree->children[i]);
00190 }
00191 }
00192 }
00193 }
00194
00195
00196
00197
00198 int Octree::skip(int x, int y, int z, int mainViewingDir)
00199 {
00200 return skip(x, y, z, mainViewingDir, root);
00201 }
00202
00203
00204
00205
00206 int Octree::skip(int x, int y, int z, int mainViewingDir, Node *tree)
00207 {
00208 int result = 1;
00209
00210 if (tree->classification == ALL_TRANSPARENT) {
00211 switch(mainViewingDir) {
00212 case XDIR: {
00213 result = tree->maxY - y;
00214 break;
00215 }
00216 case YDIR: {
00217 result = tree->maxZ - z;
00218 break;
00219 }
00220 case ZDIR: {
00221 result = tree->maxX - x;
00222 break;
00223 }
00224 }
00225 } else {
00226
00227 if (tree->classification == COMBINATION) {
00228
00229 int midX = (tree->minX + tree->maxX) / 2;
00230 int midY = (tree->minY + tree->maxY) / 2;
00231 int midZ = (tree->minZ + tree->maxZ) / 2;
00232
00233 if (x <= midX) {
00234 if (y <= midY) {
00235 if (z <= midZ) {
00236 result = skip(x, y, z, mainViewingDir, tree->children[0]);
00237 } else {
00238 result = skip(x, y, z, mainViewingDir, tree->children[4]);
00239 }
00240 } else {
00241 if (z <= midZ) {
00242 result = skip(x, y, z, mainViewingDir, tree->children[2]);
00243 } else {
00244 result = skip(x, y, z, mainViewingDir, tree->children[6]);
00245 }
00246 }
00247 } else {
00248 if (y <= midY) {
00249 if (z <= midZ) {
00250 result = skip(x, y, z, mainViewingDir, tree->children[1]);
00251 } else {
00252 result = skip(x, y, z, mainViewingDir, tree->children[5]);
00253 }
00254 } else {
00255 if (z <= midZ) {
00256 result = skip(x, y, z, mainViewingDir, tree->children[3]);
00257 } else {
00258 result = skip(x, y, z, mainViewingDir, tree->children[7]);
00259 }
00260 }
00261 }
00262 } else {
00263 result = 1;
00264 }
00265 }
00266
00267 if (result <= 0) result = 1;
00268
00269 return result;
00270 }
00271
00272
00273
00274
00275 SummedAreaTable::SummedAreaTable()
00276 {
00277 table = new int[256];
00278 }
00279
00280
00281
00282
00283 SummedAreaTable::~SummedAreaTable()
00284 {
00285 if (table) delete [] table;
00286 table = NULL;
00287 }
00288
00289
00290
00291
00292 void SummedAreaTable::build(float *transferFunction, float threshold_runlength)
00293 {
00294 byte current = 0;
00295 for (int i = 0; i < 256; ++i) {
00296 if (transferFunction[i] > threshold_runlength) {
00297 current++;
00298 }
00299 table[i] = current;
00300 }
00301 }
00302
00303
00304
00305
00306 int SummedAreaTable::transparent(byte min, byte max)
00307 {
00308 int result;
00309
00310 if (min == max) {
00311 if (max == 0) {
00312 if (table[max] >= 1) result = 0; else result = 1;
00313 } else {
00314 if (table[max] > table[max-1]) result = 0; else result = 1;
00315 }
00316 } else {
00317 if (table[max] == table[min]) {
00318 result = 1;
00319 } else {
00320 result = 0;
00321 }
00322 }
00323
00324 return result;
00325 }
00326
00327
00328
00329
00330 FastClassification::FastClassification(byte *volumeData, int dim1Size, int dim2Size, int dim3Size)
00331 {
00332 sTable = new SummedAreaTable();
00333 octree = new Octree(volumeData, dim1Size, dim2Size, dim3Size);
00334
00335 width = dim1Size;
00336 height = dim2Size;
00337 depth = dim3Size;
00338
00339 data = volumeData;
00340 }
00341
00342
00343
00344
00345 FastClassification::~FastClassification()
00346 {
00347 if (octree) delete(octree);
00348 if (sTable) delete(sTable);
00349 octree = NULL;
00350 sTable = NULL;
00351 }
00352
00353
00354
00355
00356 void FastClassification::classify()
00357 {
00358 octree->classify(sTable);
00359 }
00360
00361
00362
00363
00364 int FastClassification::skip(int x, int y, int z, int mainViewingDirection)
00365 {
00366 int result;
00367
00368 result = octree->skip(x, y, z, mainViewingDirection);
00369
00370 return result;
00371 }
00372
00373
00374
00375
00376 void FastClassification::buildSummedAreaTable(float *transferFunction, float threshold_runlength)
00377 {
00378 sTable->build(transferFunction,threshold_runlength);
00379 }
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393