00001
00002
00003
00004
00005
00006 #include "CVolumeStatistics.h"
00007 #include <QProgressDialog>
00008
00009
00010
00011
00012 CVector3f CVolumeStatistics::scale = CVector3f();
00013
00014
00015
00016
00017
00018 CVolumeStatistics::CVolumeStatistics(CVolumeFormat & pFormat)
00019 {
00020 maxVolumeDataV = 0;
00021 minVolumeDataV = USHRT_MAX;
00022
00023
00024 if(!ACCHist(&pFormat))
00025 {
00026
00027 Hist(&pFormat);
00028 }
00029 }
00030
00031
00032 CVolumeStatistics::~CVolumeStatistics()
00033 {
00034
00035 }
00036
00037
00038 void CVolumeStatistics::Hist(CVolumeFormat * const pFormat)
00039 {
00040
00041 if(pFormat->getBitDepth() == 8)
00042 {
00043 uchar nmbOfBuckets = UCHAR_MAX;
00044 uchar * voxelData = (uchar*)pFormat->getVolumeData();
00045 uint nmbOfVoxels = pFormat->getVolumeDataSize();
00046
00047
00048 originalH.data = new float[nmbOfBuckets];
00049 ZeroMemory(originalH.data,nmbOfBuckets*sizeof(originalH.data[0]));
00050
00051
00052 originalH.nmbBuckets = nmbOfBuckets;
00053
00054
00055 for(int i = 0 ; i < nmbOfVoxels ; i++)
00056 {
00057
00058 originalH.data[voxelData[i]]++;
00059
00060
00061 if(voxelData[i] > maxVolumeDataV) maxVolumeDataV = voxelData[i];
00062 if(voxelData[i] < minVolumeDataV) minVolumeDataV = voxelData[i];
00063 }
00064
00065
00066 SetMinMaxBucket(originalH);
00067 }
00068
00069
00070 if(pFormat->getBitDepth() == 12 || pFormat->getBitDepth() == 16)
00071 {
00072 ushort nmbOfBuckets = USHRT_MAX;
00073 ushort * voxelData = (ushort *)pFormat->getVolumeData();
00074 uint nmbOfVoxels = pFormat->getVolumeDataSize()/2;
00075
00076
00077 originalH.data = new float[nmbOfBuckets];
00078 ZeroMemory(originalH.data,nmbOfBuckets*sizeof(originalH.data[0]));
00079
00080
00081 originalH.nmbBuckets = nmbOfBuckets;
00082
00083
00084 for(int i = 0 ; i < nmbOfVoxels ; i++)
00085 {
00086
00087 ++originalH.data[voxelData[i]];
00088
00089
00090 if(voxelData[i] > maxVolumeDataV) maxVolumeDataV = voxelData[i];
00091 if(voxelData[i] < minVolumeDataV) minVolumeDataV = voxelData[i];
00092 }
00093
00094
00095 SetMinMaxBucket(originalH);
00096 }
00097 }
00098
00099
00100 bool CVolumeStatistics::ACCHist(CVolumeFormat * const pFormat)
00101 {
00102
00103 if(pFormat->getBitDepth() == 8)
00104 {
00105 Hist(pFormat);
00106 }
00107
00108
00109 if(pFormat->getBitDepth() == 12 || pFormat->getBitDepth() == 16)
00110 {
00111 ushort nmbOfBuckets = USHRT_MAX;
00112 ushort * voxelData = (ushort *)pFormat->getVolumeData();
00113 uint nmbOfVoxels = pFormat->getVolumeDataSize()/2;
00114
00115
00116 originalH.data = new float[nmbOfBuckets];
00117 ZeroMemory(originalH.data,nmbOfBuckets*sizeof(originalH.data[0]));
00118
00119
00120 originalH.nmbBuckets = nmbOfBuckets;
00121
00122 int X = pFormat->getSizeX();
00123 int Y = pFormat->getSizeY();
00124 int Z = pFormat->getSizeZ();
00125
00126
00127 QProgressDialog progress;
00128 progress.setLabelText("Computing ACC Histogram");
00129 progress.setRange(0,Z);
00130 progress.setModal(true);
00131
00132
00133 for(int i = 0 ; i < nmbOfVoxels ; i++)
00134 {
00135 if(voxelData[i] > maxVolumeDataV) maxVolumeDataV = voxelData[i];
00136 if(voxelData[i] < minVolumeDataV) minVolumeDataV = voxelData[i];
00137 }
00138
00139
00140
00141
00142 float tmp = pow(float(maxVolumeDataV - minVolumeDataV),2);
00143 originalGradH.nmbBuckets = (int)sqrt(tmp+tmp+tmp)+1;
00144
00145
00146 originalGradH.data = new float[originalGradH.nmbBuckets];
00147 ZeroMemory(originalGradH.data,originalGradH.nmbBuckets*sizeof(originalGradH.data[0]));
00148
00149
00150 scale.m_X = pFormat->scSX; scale.m_Y = pFormat->scSY; scale.m_Z = pFormat->scSZ;
00151
00152
00153 if(scale.m_X == 1.0 && scale.m_Y == 1.0 && scale.m_Z == 1.0)
00154 {
00155 GradientPtr = Gradient;
00156 }
00157 else
00158 {
00159 GradientPtr = GradientLI;
00160 }
00161
00162
00163
00164
00165 int lO = X*Y;
00166 int rO = Y;
00167
00168
00169 ushort cube[8] = {0};
00170
00171 int cubeMin,cubeMax;
00172
00173
00174 float gradMag[8] = {0};
00175
00176 int gradMin,gradMax;
00177
00178
00179 int z1,z0,z2,y1,y0,y2,x1,x0,x2,zlO,z1lO,z0lO,z2lO,yrO,y1rO,y0rO,y2rO;
00180
00181
00182 for(int z = 1 ; z < Z-2; z++)
00183 {
00184
00185 if(progress.wasCanceled())
00186 {
00187
00188 delete [] originalH.data;
00189
00190 ZeroMemory(originalGradH.data,originalGradH.nmbBuckets*sizeof(originalGradH.data[0]));
00191
00192 return false;
00193 }
00194
00195 progress.setValue(z);
00196
00197 for(int y = 1 ; y < Y-2; y++)
00198 {
00199 qApp->processEvents();
00200
00201 for(int x = 1 ; x < X-2; x++)
00202 {
00203 x0 = x-1; x1 = x+1 ; x2 = x+2;
00204 y0 = y-1; y1 = y+1 ; y2 = y+2;
00205 z0 = z-1; z1 = z+1 ; z2 = z+2;
00206
00207 z0lO = z0*lO;
00208 zlO = z*lO;
00209 z1lO = z1*lO;
00210 z2lO = z2*lO;
00211
00212 y0rO = y0*rO;
00213 yrO = y*rO;
00214 y1rO = y1*rO;
00215 y2rO = y2*rO;
00216
00217 cube[0] = voxelData[zlO + yrO + x];
00218 cube[1] = voxelData[zlO + yrO + x1];
00219 cube[2] = voxelData[zlO + y1rO + x];
00220 cube[3] = voxelData[zlO + y1rO + x1];
00221 cube[4] = voxelData[z1lO + yrO + x];
00222 cube[5] = voxelData[z1lO + yrO + x1];
00223 cube[6] = voxelData[z1lO + y1rO + x];
00224 cube[7] = voxelData[z1lO + y1rO + x1];
00225
00226
00227 gradMag[0] = (*GradientPtr)(voxelData[zlO + yrO + x0],cube[1],
00228 voxelData[zlO + y0rO + x ],cube[2],
00229 voxelData[z0lO + yrO + x ],cube[4],cube[0]);
00230
00231
00232 gradMag[1] = (*GradientPtr)(cube[0],voxelData[zlO + yrO + x2],
00233 voxelData[zlO + y0rO + x1],cube[3],
00234 voxelData[z0lO + yrO + x1],cube[5],cube[1]);
00235
00236
00237 gradMag[2] = (*GradientPtr)(voxelData[zlO + y1rO + x0],cube[3],
00238 cube[0],voxelData[zlO + y2rO + x],
00239 voxelData[z0lO + y1rO + x],cube[6],cube[2]);
00240
00241
00242 gradMag[3] = (*GradientPtr)(cube[2],voxelData[zlO + y1rO + x2],
00243 cube[1],voxelData[zlO + y2rO + x1],
00244 voxelData[z0lO + y1rO + x1],cube[7],cube[3]);
00245
00246
00247 gradMag[4] = (*GradientPtr)(voxelData[z1lO + yrO + x0],cube[5],
00248 voxelData[z1lO + y0rO + x ],cube[6],
00249 cube[0],voxelData[z2lO + yrO + x],cube[4]);
00250
00251
00252 gradMag[5] = (*GradientPtr)(cube[4],voxelData[z1lO + yrO + x2],
00253 voxelData[z1lO + y0rO + x1],cube[7],
00254 cube[1],voxelData[z2lO + yrO + x1],cube[5]);
00255
00256
00257 gradMag[6] = (*GradientPtr)(voxelData[z1lO + y1rO + x0],cube[7],
00258 cube[4],voxelData[z1lO + y2rO + x],
00259 cube[2],voxelData[z2lO + y1rO + x],cube[6]);
00260
00261
00262 gradMag[7] = (*GradientPtr)(cube[6],voxelData[z1lO + y1rO + x2],
00263 cube[5],voxelData[zlO + y2rO + x1],
00264 cube[3],voxelData[z2lO + y1rO + x1],cube[7]);
00265
00266 cubeMin = USHRT_MAX;
00267 cubeMax = 0;
00268 gradMin = originalGradH.nmbBuckets;
00269 gradMax = 0;
00270
00271 for(int i = 0 ; i < 8 ; i++)
00272 {
00273 if(cube[i] > cubeMax) cubeMax = cube[i];
00274 if(cube[i] < cubeMin) cubeMin = cube[i];
00275 if(gradMag[i] > gradMax) gradMax = gradMag[i];
00276 if(gradMag[i] < gradMin) gradMin = gradMag[i];
00277 }
00278
00279 for(int i = cubeMin ; i <= cubeMax ; i++)
00280 {
00281 originalH.data[i]++;
00282 }
00283 for(int i = (int)gradMin ; i <= (int)gradMax ; i++)
00284 {
00285 originalGradH.data[i]++;
00286 }
00287 }
00288 }
00289 }
00290
00291
00292 SetMinMaxBucket(originalH);
00293 SetMinMaxBucket(originalGradH);
00294 }
00295 return true;
00296 }
00297
00298
00299 SHistogramData & CVolumeStatistics::GetHistogram(int quantisation, char normalizeMethod)
00300 {
00301 if(quantizedH.nmbBuckets != quantisation)
00302 {
00303
00304 SAFE_DELETE_ARRAY(quantizedH.data);
00305
00306
00307 quantizedH.data = new float[quantisation];
00308 ZeroMemory(quantizedH.data,quantisation*sizeof(quantizedH.data[0]));
00309
00310
00311 quantizedH.nmbBuckets = quantisation;
00312
00313
00314 double step = quantisation / (double)originalH.nmbBuckets;
00315 for(int i = 0 ; i < originalH.nmbBuckets ; i++){
00316 quantizedH.data[int(step*i)] += originalH.data[i];
00317 }
00318
00319 SetMinMaxBucket(quantizedH);
00320
00321
00322 if(normQuantizedH.nmbBuckets != quantisation)
00323 {
00324
00325 SAFE_DELETE_ARRAY(normQuantizedH.data);
00326
00327
00328 normQuantizedH.data = new float[quantisation];
00329
00330
00331 normQuantizedH.nmbBuckets = quantisation;
00332 }
00333
00334 if(normalizeMethod != NORMALIZE_NONE)
00335 {
00336
00337 Normalize(quantizedH,normQuantizedH,normalizeMethod);
00338 normQuantizedH.normalizeMethod = normalizeMethod;
00339 }
00340 }
00341
00342 if(normQuantizedH.normalizeMethod != normalizeMethod)
00343 {
00344
00345 Normalize(quantizedH,normQuantizedH,normalizeMethod);
00346 normQuantizedH.normalizeMethod = normalizeMethod;
00347 }
00348
00349 if(normalizeMethod == NORMALIZE_NONE){
00350 return quantizedH;
00351 }else{
00352 return normQuantizedH;
00353 }
00354 }
00355
00356
00357 SHistogramData & CVolumeStatistics::GetGradHistogram(int quantisation, char normalizeMethod)
00358 {
00359 if(quantizedGradH.nmbBuckets != quantisation)
00360 {
00361
00362 SAFE_DELETE_ARRAY(quantizedGradH.data);
00363
00364
00365 quantizedGradH.data = new float[quantisation];
00366 ZeroMemory(quantizedGradH.data,quantisation*sizeof(quantizedGradH.data[0]));
00367
00368
00369 quantizedGradH.nmbBuckets = quantisation;
00370
00371
00372 double step = quantisation / (double)originalGradH.nmbBuckets;
00373 for(int i = 0 ; i < originalGradH.nmbBuckets ; i++){
00374 quantizedGradH.data[int(step*i)] += originalGradH.data[i];
00375 }
00376
00377 SetMinMaxBucket(quantizedGradH);
00378
00379
00380 if(normQuantizedGradH.nmbBuckets != quantisation)
00381 {
00382
00383 SAFE_DELETE_ARRAY(normQuantizedGradH.data);
00384
00385
00386 normQuantizedGradH.data = new float[quantisation];
00387
00388
00389 normQuantizedGradH.nmbBuckets = quantisation;
00390 }
00391
00392 if(normalizeMethod != NORMALIZE_NONE)
00393 {
00394
00395 Normalize(quantizedGradH,normQuantizedGradH,normalizeMethod);
00396 normQuantizedGradH.normalizeMethod = normalizeMethod;
00397 }
00398 }
00399
00400 if(normQuantizedGradH.normalizeMethod != normalizeMethod)
00401 {
00402
00403 Normalize(quantizedGradH,normQuantizedGradH,normalizeMethod);
00404 normQuantizedGradH.normalizeMethod = normalizeMethod;
00405 }
00406
00407 if(normalizeMethod == NORMALIZE_NONE){
00408 return quantizedGradH;
00409 }else{
00410 return normQuantizedGradH;
00411 }
00412 }
00413
00414
00415
00416
00417
00418 void CVolumeStatistics::SetMinMaxBucket(SHistogramData & hist)
00419 {
00420
00421 for(int i = 0 ; i < hist.nmbBuckets ; i++)
00422 {
00423 if(hist.data[i] > hist.maxBucket) hist.maxBucket = hist.data[i];
00424 if(hist.data[i] < hist.minBucket) hist.minBucket = hist.data[i];
00425 }
00426 }
00427
00428 void CVolumeStatistics::Normalize(SHistogramData & quantizedH, SHistogramData & normQuantizedH, char normalizeMethod)
00429 {
00430 switch(normalizeMethod)
00431 {
00432 case NORMALIZE_LINEAR:
00433 {
00434
00435 float range = quantizedH.maxBucket - quantizedH.minBucket;
00436 for(int i = 0 ; i < normQuantizedH.nmbBuckets ; i++){
00437 normQuantizedH.data[i] = (quantizedH.data[i] - quantizedH.minBucket) / range;
00438 }
00439 }
00440 break;
00441 case NORMALIZE_LOG:
00442 {
00443
00444 for(int i = 0 ; i < normQuantizedH.nmbBuckets ; i++){
00445 normQuantizedH.data[i] = log(quantizedH.data[i]+1);
00446 if(normQuantizedH.data[i] > normQuantizedH.maxBucket) normQuantizedH.maxBucket = normQuantizedH.data[i];
00447 if(normQuantizedH.data[i] < normQuantizedH.minBucket) normQuantizedH.minBucket = normQuantizedH.data[i];
00448 }
00449
00450 float rangeL = normQuantizedH.maxBucket - normQuantizedH.minBucket;
00451 for(int i = 0 ; i < normQuantizedH.nmbBuckets ; i++){
00452 normQuantizedH.data[i] = (normQuantizedH.data[i] - normQuantizedH.minBucket) / rangeL;
00453 }
00454 }
00455 break;
00456 }
00457 }
00458
00459 float CVolumeStatistics::Gradient(int x0, int x1, int y0, int y1, int z0, int z1,int center)
00460 {
00461 return sqrt(pow(double(x1-x0),2)+pow(double(y1-y0),2)+pow(double(z1-z0),2));
00462 }
00463
00464 float CVolumeStatistics::GradientLI(int x0, int x1, int y0, int y1, int z0, int z1,int center)
00465 {
00466 float i1,i0,dx,dy,dz;
00467
00468 i1 = scale.m_X*x1 - (1-scale.m_X)*center;
00469 i0 = scale.m_X*x0 - (1-scale.m_X)*center;
00470 dx = (i1-i0)*(i1-i0);
00471
00472 i1 = scale.m_Y*y1 - (1-scale.m_Y)*center;
00473 i0 = scale.m_Y*y0 - (1-scale.m_Y)*center;
00474 dy = (i1-i0)*(i1-i0);
00475
00476 i1 = scale.m_Z*z1 - (1-scale.m_Z)*center;
00477 i0 = scale.m_Z*z0 - (1-scale.m_Z)*center;
00478 dz = (i1-i0)*(i1-i0);
00479
00480 return sqrt((float)dx+dy+dz);
00481 }