00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <math.h>
00026 #include <string.h>
00027 #include "VolumetricData.h"
00028 #include "Matrix4.h"
00029 #include "utilities.h"
00030 #include "Segmentation.h"
00031 #include <float.h>
00032
00033 #ifndef NAN //not a number
00034 const float NAN = sqrtf(-1.f);
00035 #endif
00036
00037 #define MIN(X,Y) (((X)<(Y))? (X) : (Y))
00038 #define MAX(X,Y) (((X)>(Y))? (X) : (Y))
00039
00041 VolumetricData::VolumetricData(const char *dataname, const float *o,
00042 const float *xa, const float *ya, const float *za,
00043 int xs, int ys, int zs, float *dataptr) {
00044 name = stringdup(dataname);
00045 data = dataptr;
00046 xsize = xs;
00047 ysize = ys;
00048 zsize = zs;
00049
00050 gradient = NULL;
00051 gradient_isvalid = false;
00052 invalidate_minmax();
00053 invalidate_mean();
00054 invalidate_sigma();
00055
00056 for (int i=0; i<3; i++) {
00057 origin[i] = (double)o[i];
00058 xaxis[i] = (double)xa[i];
00059 yaxis[i] = (double)ya[i];
00060 zaxis[i] = (double)za[i];
00061 }
00062
00063 compute_minmax();
00064 }
00065
00066
00068 VolumetricData::VolumetricData(const char *dataname, const double *o,
00069 const double *xa, const double *ya, const double *za,
00070 int xs, int ys, int zs, float *dataptr) {
00071 name = stringdup(dataname);
00072 data = dataptr;
00073 xsize = xs;
00074 ysize = ys;
00075 zsize = zs;
00076
00077 gradient = NULL;
00078 gradient_isvalid = false;
00079 invalidate_minmax();
00080 invalidate_mean();
00081 invalidate_sigma();
00082
00083 for (int i=0; i<3; i++) {
00084 origin[i] = o[i];
00085 xaxis[i] = xa[i];
00086 yaxis[i] = ya[i];
00087 zaxis[i] = za[i];
00088 }
00089
00090 compute_minmax();
00091 }
00092
00093
00095 VolumetricData::~VolumetricData() {
00096 delete [] name;
00097 delete [] data;
00098 delete [] gradient;
00099 }
00100
00101
00104 void VolumetricData::set_name(const char *newname) {
00105 if (name) delete[] name;
00106 name = new char[strlen(newname)+1];
00107 strcpy(name, newname);
00108 }
00109
00110
00112 void VolumetricData::cell_lengths(float *xl, float *yl, float *zl) const {
00113 float xsinv, ysinv, zsinv;
00114
00115
00116 if (xsize > 1)
00117 xsinv = 1.0f / (xsize - 1.0f);
00118 else
00119 xsinv = 1.0f;
00120
00121 if (ysize > 1)
00122 ysinv = 1.0f / (ysize - 1.0f);
00123 else
00124 ysinv = 1.0f;
00125
00126 if (zsize > 1)
00127 zsinv = 1.0f / (zsize - 1.0f);
00128 else
00129 zsinv = 1.0f;
00130
00131 *xl = sqrtf(float(dot_prod(xaxis, xaxis))) * xsinv;
00132 *yl = sqrtf(float(dot_prod(yaxis, yaxis))) * ysinv;
00133 *zl = sqrtf(float(dot_prod(zaxis, zaxis))) * zsinv;
00134 }
00135
00136
00138 void VolumetricData::cell_axes(float *xax, float *yax, float *zax) const {
00139 float xsinv, ysinv, zsinv;
00140
00141
00142 if (xsize > 1)
00143 xsinv = 1.0f / (xsize - 1.0f);
00144 else
00145 xsinv = 1.0f;
00146
00147 if (ysize > 1)
00148 ysinv = 1.0f / (ysize - 1.0f);
00149 else
00150 ysinv = 1.0f;
00151
00152 if (zsize > 1)
00153 zsinv = 1.0f / (zsize - 1.0f);
00154 else
00155 zsinv = 1.0f;
00156
00157 xax[0] = float(xaxis[0] * xsinv);
00158 xax[1] = float(xaxis[1] * xsinv);
00159 xax[2] = float(xaxis[2] * xsinv);
00160
00161 yax[0] = float(yaxis[0] * ysinv);
00162 yax[1] = float(yaxis[1] * ysinv);
00163 yax[2] = float(yaxis[2] * ysinv);
00164
00165 zax[0] = float(zaxis[0] * zsinv);
00166 zax[1] = float(zaxis[1] * zsinv);
00167 zax[2] = float(zaxis[2] * zsinv);
00168 }
00169
00170
00172 void VolumetricData::cell_axes(double *xax, double *yax, double *zax) const {
00173 double xsinv, ysinv, zsinv;
00174
00175
00176 if (xsize > 1)
00177 xsinv = 1.0 / (xsize - 1.0);
00178 else
00179 xsinv = 1.0;
00180
00181 if (ysize > 1)
00182 ysinv = 1.0 / (ysize - 1.0);
00183 else
00184 ysinv = 1.0;
00185
00186 if (zsize > 1)
00187 zsinv = 1.0 / (zsize - 1.0);
00188 else
00189 zsinv = 1.0;
00190
00191 xax[0] = xaxis[0] * xsinv;
00192 xax[1] = xaxis[1] * xsinv;
00193 xax[2] = xaxis[2] * xsinv;
00194
00195 yax[0] = yaxis[0] * ysinv;
00196 yax[1] = yaxis[1] * ysinv;
00197 yax[2] = yaxis[2] * ysinv;
00198
00199 zax[0] = zaxis[0] * zsinv;
00200 zax[1] = zaxis[1] * zsinv;
00201 zax[2] = zaxis[2] * zsinv;
00202 }
00203
00204
00206 void VolumetricData::cell_dirs(float *xad, float *yad, float *zad) const {
00207 float xl, yl, zl;
00208
00209 cell_lengths(&xl, &yl, &zl);
00210 cell_axes(xad, yad, zad);
00211
00212 xad[0] /= xl;
00213 xad[1] /= xl;
00214 xad[2] /= xl;
00215
00216 yad[0] /= yl;
00217 yad[1] /= yl;
00218 yad[2] /= yl;
00219
00220 zad[0] /= zl;
00221 zad[1] /= zl;
00222 zad[2] /= zl;
00223 }
00224
00225
00227 double VolumetricData::cell_volume() const {
00228 double c[3][3];
00229 cell_axes(c[0], c[1], c[2]);
00230
00231 double vol =
00232 ( c[0][0] * (c[1][1]*c[2][2] - c[2][1]*c[1][2]))
00233 - (c[1][0] * (c[0][1]*c[2][2] - c[2][1]*c[0][2]))
00234 + (c[2][0] * (c[0][1]*c[1][2] - c[1][1]*c[0][2]));
00235
00236 return vol;
00237 }
00238
00239
00242 void VolumetricData::voxel_coord_from_cartesian_coord(const float *carcoord, float *voxcoord, int shiftflag) const {
00243 float min_coord[3];
00244
00245 if (shiftflag) {
00246
00247 for (int i=0; i<3; i++)
00248 min_coord[i] = (float) origin[i] - 0.5f*float(xaxis[i]/(xsize-1) + yaxis[i]/(ysize-1) + zaxis[i]/(zsize-1));
00249 } else {
00250 for (int i=0; i<3; i++)
00251 min_coord[i] = (float) origin[i];
00252 }
00253
00254
00255 float matval[16];
00256 matval[0 ] = (float) xaxis[0]/(xsize-1);
00257 matval[1 ] = (float) yaxis[0]/(ysize-1);
00258 matval[2 ] = (float) zaxis[0]/(zsize-1);
00259 matval[3 ] = 0.f;
00260 matval[4 ] = (float) xaxis[1]/(xsize-1);
00261 matval[5 ] = (float) yaxis[1]/(ysize-1);
00262 matval[6 ] = (float) zaxis[1]/(zsize-1);
00263 matval[7 ] = 0.f;
00264 matval[8 ] = (float) xaxis[2]/(xsize-1);
00265 matval[9 ] = (float) yaxis[2]/(ysize-1);
00266 matval[10] = (float) zaxis[2]/(zsize-1);
00267 matval[11] = 0.f;
00268 matval[12] = (float) min_coord[0];
00269 matval[13] = (float) min_coord[1];
00270 matval[14] = (float) min_coord[2];
00271 matval[15] = 1.f;
00272 Matrix4 matrix(matval);
00273 matrix.inverse();
00274
00275 matrix.multpoint3d(carcoord, voxcoord);
00276 }
00277
00278
00280 ptrdiff_t VolumetricData::voxel_index_from_coord(float xpos, float ypos, float zpos) const {
00281 float realcoord[3];
00282 float gridpoint[3];
00283 realcoord[0] = xpos;
00284 realcoord[1] = ypos;
00285 realcoord[2] = zpos;
00286
00287
00288 voxel_coord_from_cartesian_coord(realcoord, gridpoint, 1);
00289
00290 int gx = (int) gridpoint[0];
00291 int gy = (int) gridpoint[1];
00292 int gz = (int) gridpoint[2];
00293 if (gx < 0 || gx >= xsize) return -1;
00294 if (gy < 0 || gy >= ysize) return -1;
00295 if (gz < 0 || gz >= zsize) return -1;
00296 return gz*ptrdiff_t(xsize)*ptrdiff_t(ysize) + gy*ptrdiff_t(xsize) + gx;
00297 }
00298
00299
00301 float VolumetricData::voxel_value_safe(int x, int y, int z) const {
00302 ptrdiff_t xx, yy, zz;
00303 xx = (x > 0) ? ((x < xsize) ? x : xsize-1) : 0;
00304 yy = (y > 0) ? ((y < ysize) ? y : ysize-1) : 0;
00305 zz = (z > 0) ? ((z < zsize) ? z : zsize-1) : 0;
00306 ptrdiff_t index = zz*xsize*ysize + yy*xsize + xx;
00307 return data[index];
00308 }
00309
00310
00312 float VolumetricData::voxel_value_interpolate(float xv, float yv, float zv) const {
00313 int x = (int) xv;
00314 int y = (int) yv;
00315 int z = (int) zv;
00316 float xf = xv - x;
00317 float yf = yv - y;
00318 float zf = zv - z;
00319 float xlerps[4];
00320 float ylerps[2];
00321 float tmp;
00322
00323 tmp = voxel_value_safe(x, y, z);
00324 xlerps[0] = tmp + xf*(voxel_value_safe(x+1, y, z) - tmp);
00325
00326 tmp = voxel_value_safe(x, y+1, z);
00327 xlerps[1] = tmp + xf*(voxel_value_safe(x+1, y+1, z) - tmp);
00328
00329 tmp = voxel_value_safe(x, y, z+1);
00330 xlerps[2] = tmp + xf*(voxel_value_safe(x+1, y, z+1) - tmp);
00331
00332 tmp = voxel_value_safe(x, y+1, z+1);
00333 xlerps[3] = tmp + xf*(voxel_value_safe(x+1, y+1, z+1) - tmp);
00334
00335 ylerps[0] = xlerps[0] + yf*(xlerps[1] - xlerps[0]);
00336 ylerps[1] = xlerps[2] + yf*(xlerps[3] - xlerps[2]);
00337
00338 return ylerps[0] + zf*(ylerps[1] - ylerps[0]);
00339 }
00340
00341
00343 float VolumetricData::voxel_value_from_coord(float xpos, float ypos, float zpos) const {
00344 ptrdiff_t ind = voxel_index_from_coord(xpos, ypos, zpos);
00345 if (ind > 0)
00346 return data[ind];
00347 else
00348 return NAN;
00349 }
00350
00351
00353 float VolumetricData::voxel_value_interpolate_from_coord(float xpos, float ypos, float zpos) const {;
00354 float realcoord[3];
00355 float gridpoint[3];
00356 realcoord[0] = xpos;
00357 realcoord[1] = ypos;
00358 realcoord[2] = zpos;
00359
00360
00361 voxel_coord_from_cartesian_coord(realcoord, gridpoint, 0);
00362
00363 int gx = (int) gridpoint[0];
00364 int gy = (int) gridpoint[1];
00365 int gz = (int) gridpoint[2];
00366 if (gx < 0 || gx >= xsize) return NAN;
00367 if (gy < 0 || gy >= ysize) return NAN;
00368 if (gz < 0 || gz >= zsize) return NAN;
00369 return voxel_value_interpolate(gridpoint[0], gridpoint[1], gridpoint[2]);
00370 }
00371
00372
00375 float VolumetricData::voxel_value_from_coord_safe(float xpos, float ypos, float zpos) const {
00376 ptrdiff_t ind = voxel_index_from_coord(xpos, ypos, zpos);
00377 if (ind > 0)
00378 return data[ind];
00379 else
00380 return 0;
00381 }
00382
00383
00386 float VolumetricData::voxel_value_interpolate_from_coord_safe(float xpos, float ypos, float zpos) const {;
00387 float realcoord[3];
00388 float gridpoint[3];
00389 realcoord[0] = xpos;
00390 realcoord[1] = ypos;
00391 realcoord[2] = zpos;
00392
00393
00394 voxel_coord_from_cartesian_coord(realcoord, gridpoint, 0);
00395
00396 int gx = (int) gridpoint[0];
00397 int gy = (int) gridpoint[1];
00398 int gz = (int) gridpoint[2];
00399 if (gx < 0 || gx >= xsize) return 0;
00400 if (gy < 0 || gy >= ysize) return 0;
00401 if (gz < 0 || gz >= zsize) return 0;
00402 return voxel_value_interpolate(gridpoint[0], gridpoint[1], gridpoint[2]);
00403 }
00404
00405
00406 void VolumetricData::invalidate_gradient() {
00407 if (gradient) {
00408 delete [] gradient;
00409 gradient = NULL;
00410 gradient_isvalid = false;
00411 }
00412 }
00413
00414
00415 const float * VolumetricData::access_volume_gradient() {
00416 if ((gradient == NULL) || (!gradient_isvalid)) {
00417 compute_volume_gradient();
00418 }
00419 return gradient;
00420 }
00421
00422
00424 void VolumetricData::set_volume_gradient(float *grad) {
00425 if (gradient) {
00426 delete [] gradient;
00427 gradient = NULL;
00428 gradient_isvalid = false;
00429 }
00430
00431 if (!gradient) {
00432 ptrdiff_t gsz = ptrdiff_t(xsize) * ptrdiff_t(ysize) * ptrdiff_t(zsize) * 3L;
00433 gradient = new float[gsz];
00434 memcpy(gradient, grad, gsz*sizeof(float));
00435 }
00436
00437 gradient_isvalid = true;
00438 }
00439
00440
00442 void VolumetricData::compute_volume_gradient(void) {
00443 int xi, yi, zi;
00444 float xs, ys, zs;
00445 float xl, yl, zl;
00446
00447 if (!gradient) {
00448 ptrdiff_t gsz = this->gridsize() * 3L;
00449 #if 0
00450 printf("xsz: %d ysz: %d zsz: %d\n", xsize, ysize, zsize);
00451 printf("gradient map sz: %ld, MB: %ld\n", gsz, gsz*sizeof(float)/(1024*1024));
00452 #endif
00453 gradient = new float[gsz];
00454 }
00455
00456
00457 cell_lengths(&xl, &yl, &zl);
00458
00459
00460
00461 xs = -0.5f / xl;
00462 ys = -0.5f / yl;
00463 zs = -0.5f / zl;
00464
00465 ptrdiff_t zplanesz = xsize * ysize;
00466 for (zi=0; zi<zsize; zi++) {
00467 int zm = clamp_int(zi - 1, 0, zsize - 1);
00468 int zp = clamp_int(zi + 1, 0, zsize - 1);
00469
00470 ptrdiff_t zmplanesz = zm*zplanesz;
00471 ptrdiff_t ziplanesz = zi*zplanesz;
00472 ptrdiff_t zpplanesz = zp*zplanesz;
00473 for (yi=0; yi<ysize; yi++) {
00474 int ym = clamp_int(yi - 1, 0, ysize - 1);
00475 int yp = clamp_int(yi + 1, 0, ysize - 1);
00476
00477 int yirowsz = yi*xsize;
00478 ptrdiff_t row = ziplanesz + yirowsz;
00479 ptrdiff_t rowym = ziplanesz + ym*xsize;
00480 ptrdiff_t rowyp = ziplanesz + yp*xsize;
00481 ptrdiff_t rowzm = zmplanesz + yirowsz;
00482 ptrdiff_t rowzp = zpplanesz + yirowsz;
00483 for (xi=0; xi<xsize; xi++) {
00484 ptrdiff_t gradidx = (row + xi) * 3L;
00485 int xm = clamp_int(xi - 1, 0, xsize - 1);
00486 int xp = clamp_int(xi + 1, 0, xsize - 1);
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497 #if 1
00498 gradient[gradidx ] = (data[row + xp] - data[row + xm]) * xs;
00499 gradient[gradidx + 1] = (data[rowyp+xi] - data[rowym+xi]) * ys;
00500 gradient[gradidx + 2] = (data[rowzp+xi] - data[rowzm+xi]) * zs;
00501 #else
00502 gradient[gradidx ] =
00503 (voxel_value(xp, yi, zi) - voxel_value(xm, yi, zi)) * xs;
00504 gradient[gradidx + 1] =
00505 (voxel_value(xi, yp, zi) - voxel_value(xi, ym, zi)) * ys;
00506 gradient[gradidx + 2] =
00507 (voxel_value(xi, yi, zp) - voxel_value(xi, yi, zm)) * zs;
00508 #endif
00509 }
00510 }
00511 }
00512
00513 gradient_isvalid = true;
00514 }
00515
00516
00518 void VolumetricData::voxel_gradient_safe(int x, int y, int z, float *grad) const {
00519 int xx, yy, zz;
00520 xx = (x > 0) ? ((x < xsize) ? x : xsize-1) : 0;
00521 yy = (y > 0) ? ((y < ysize) ? y : ysize-1) : 0;
00522 zz = (z > 0) ? ((z < zsize) ? z : zsize-1) : 0;
00523 ptrdiff_t index = zz*xsize*ysize + yy*xsize + xx;
00524 grad[0] = gradient[index*3L ];
00525 grad[1] = gradient[index*3L + 1];
00526 grad[2] = gradient[index*3L + 2];
00527 }
00528
00529
00531 void VolumetricData::voxel_gradient_interpolate(const float *voxcoord, float *grad) const {
00532 float gtmpa[3], gtmpb[3];
00533 int x = (int) voxcoord[0];
00534 int y = (int) voxcoord[1];
00535 int z = (int) voxcoord[2];
00536 float xf = voxcoord[0] - x;
00537 float yf = voxcoord[1] - y;
00538 float zf = voxcoord[2] - z;
00539 float xlerpsa[3];
00540 float xlerpsb[3];
00541 float xlerpsc[3];
00542 float xlerpsd[3];
00543 float ylerpsa[3];
00544 float ylerpsb[3];
00545
00546 voxel_gradient_safe(x, y, z, gtmpa);
00547 voxel_gradient_safe(x+1, y, z, gtmpb);
00548 vec_lerp(xlerpsa, gtmpa, gtmpb, xf);
00549
00550 voxel_gradient_safe(x, y+1, z, gtmpa);
00551 voxel_gradient_safe(x+1, y+1, z, gtmpb);
00552 vec_lerp(xlerpsb, gtmpa, gtmpb, xf);
00553
00554 voxel_gradient_safe(x, y, z+1, gtmpa);
00555 voxel_gradient_safe(x+1, y, z+1, gtmpb);
00556 vec_lerp(xlerpsc, gtmpa, gtmpb, xf);
00557
00558 voxel_gradient_safe(x, y+1, z+1, gtmpa);
00559 voxel_gradient_safe(x+1, y+1, z+1, gtmpb);
00560 vec_lerp(xlerpsd, gtmpa, gtmpb, xf);
00561
00562 vec_lerp(ylerpsa, xlerpsa, xlerpsb, yf);
00563 vec_lerp(ylerpsb, xlerpsc, xlerpsd, yf);
00564
00565 vec_lerp(grad, ylerpsa, ylerpsb, zf);
00566 }
00567
00568
00570 void VolumetricData::voxel_gradient_from_coord(const float *coord, float *grad) const {
00571 int vx, vy, vz;
00572 float voxcoord[3];
00573 voxel_coord_from_cartesian_coord(coord, voxcoord, 1);
00574 vx = (int) voxcoord[0];
00575 vy = (int) voxcoord[1];
00576 vz = (int) voxcoord[2];
00577 voxel_gradient_safe(vx, vy, vz, grad);
00578 }
00579
00580
00582 void VolumetricData::voxel_gradient_interpolate_from_coord(const float *coord, float *grad) const {
00583 float voxcoord[3];
00584 voxel_coord_from_cartesian_coord(coord, voxcoord, 0);
00585
00586 if (voxcoord[0] < 0 || voxcoord[0] >= (xsize-1) ||
00587 voxcoord[1] < 0 || voxcoord[1] >= (ysize-1) ||
00588 voxcoord[2] < 0 || voxcoord[2] >= (zsize-1)) {
00589 grad[0] = NAN;
00590 grad[1] = NAN;
00591 grad[2] = NAN;
00592 return;
00593 }
00594
00595 voxel_gradient_interpolate(voxcoord, grad);
00596 }
00597
00598
00599
00600
00601 void VolumetricData::pad(int padxm, int padxp, int padym, int padyp, int padzm, int padzp) {
00602 double xdelta[3] = {(xaxis[0] / (xsize - 1)), (xaxis[1] / (xsize - 1)), (xaxis[2] / (xsize - 1))};
00603 double ydelta[3] = {(yaxis[0] / (ysize - 1)), (yaxis[1] / (ysize - 1)), (yaxis[2] / (ysize - 1))};
00604 double zdelta[3] = {(zaxis[0] / (zsize - 1)), (zaxis[1] / (zsize - 1)), (zaxis[2] / (zsize - 1))};
00605
00606 int xsize_new = MAX(1, xsize + padxm + padxp);
00607 int ysize_new = MAX(1, ysize + padym + padyp);
00608 int zsize_new = MAX(1, zsize + padzm + padzp);
00609
00610 ptrdiff_t newgridsize = ptrdiff_t(xsize_new)*ptrdiff_t(ysize_new)*ptrdiff_t(zsize_new);
00611 float *data_new = new float[newgridsize];
00612 memset(data_new, 0, newgridsize*sizeof(float));
00613
00614 int startx = MAX(0, padxm);
00615 int starty = MAX(0, padym);
00616 int startz = MAX(0, padzm);
00617 int endx = MIN(xsize_new, xsize+padxm);
00618 int endy = MIN(ysize_new, ysize+padym);
00619 int endz = MIN(zsize_new, zsize+padzm);
00620
00621 int gx, gy, gz;
00622 for (gz=startz; gz<endz; gz++) {
00623 for (gy=starty; gy<endy; gy++) {
00624 ptrdiff_t oldyzaddrminpadxm = (gy-padym)*ptrdiff_t(xsize) +
00625 (gz-padzm)*ptrdiff_t(xsize)*ptrdiff_t(ysize) - padxm;
00626 ptrdiff_t newyzaddr = gy*ptrdiff_t(xsize_new) + gz*ptrdiff_t(xsize_new)*ptrdiff_t(ysize_new);
00627 for (gx=startx; gx<endx; gx++) {
00628 data_new[gx + newyzaddr] = data[gx + oldyzaddrminpadxm];
00629 }
00630 }
00631 }
00632
00633 delete data;
00634 data = data_new;
00635
00636 xsize = xsize_new;
00637 ysize = ysize_new;
00638 zsize = zsize_new;
00639
00640 vec_scaled_add(xaxis, float(padxm+padxp), xdelta);
00641 vec_scaled_add(yaxis, float(padym+padyp), ydelta);
00642 vec_scaled_add(zaxis, float(padzm+padzp), zdelta);
00643
00644 vec_scaled_add(origin, float(-padxm), xdelta);
00645 vec_scaled_add(origin, float(-padym), ydelta);
00646 vec_scaled_add(origin, float(-padzm), zdelta);
00647
00648
00649
00650
00651
00652
00653
00654
00655 invalidate_minmax();
00656
00657
00658
00659
00660
00661 invalidate_mean();
00662 invalidate_sigma();
00663
00664
00665
00666
00667 invalidate_gradient();
00668 }
00669
00670
00671
00672
00673
00674 void VolumetricData::crop(double crop_minx, double crop_miny, double crop_minz, double crop_maxx, double crop_maxy, double crop_maxz) {
00675 double xdelta[3] = {(xaxis[0] / (xsize - 1)), (xaxis[1] / (xsize - 1)), (xaxis[2] / (xsize - 1))};
00676 double ydelta[3] = {(yaxis[0] / (ysize - 1)), (yaxis[1] / (ysize - 1)), (yaxis[2] / (ysize - 1))};
00677 double zdelta[3] = {(zaxis[0] / (zsize - 1)), (zaxis[1] / (zsize - 1)), (zaxis[2] / (zsize - 1))};
00678
00679
00680 int padxm = int((origin[0] - crop_minx)/xdelta[0]);
00681 int padym = int((origin[1] - crop_miny)/ydelta[1]);
00682 int padzm = int((origin[2] - crop_minz)/zdelta[2]);
00683
00684 int padxp = int((crop_maxx - origin[0] - xaxis[0])/xdelta[0]);
00685 int padyp = int((crop_maxy - origin[1] - yaxis[1])/ydelta[1]);
00686 int padzp = int((crop_maxz - origin[2] - zaxis[2])/zdelta[2]);
00687
00688
00689 pad(padxm, padxp, padym, padyp, padzm, padzp);
00690 }
00691
00692
00693 void VolumetricData::clamp(float min_value, float max_value) {
00694
00695 ptrdiff_t gsz = gridsize();
00696 for (ptrdiff_t i=0; i<gsz; i++) {
00697 #if 1
00698
00699 const float d = data[i];
00700 float t = (d < min_value) ? min_value : d;
00701 data[i] = (t > max_value) ? max_value : t;
00702 #else
00703
00704 if (data[i] < min_value) data[i] = min_value;
00705 else if (data[i] > max_value) data[i] = max_value;
00706 #endif
00707 }
00708
00709
00710 if (cached_min < min_value)
00711 cached_min = min_value;
00712
00713 if (cached_max > max_value)
00714 cached_max = max_value;
00715
00716 invalidate_mean();
00717 invalidate_sigma();
00718
00719
00720 invalidate_gradient();
00721 }
00722
00723
00724 void VolumetricData::scale_by(float ff) {
00725 ptrdiff_t gsz = gridsize();
00726 for (ptrdiff_t i=0; i<gsz; i++) {
00727 data[i] *= ff;
00728 }
00729
00730
00731 if (ff >= 0.0) {
00732 cached_min *= ff;
00733 cached_max *= ff;
00734 } else {
00735
00736 float tmp = cached_min;
00737 cached_min = cached_max * ff;
00738 cached_max = tmp * ff;
00739 }
00740
00741 invalidate_mean();
00742 invalidate_sigma();
00743
00744
00745 invalidate_gradient();
00746 }
00747
00748
00749 void VolumetricData::scalar_add(float ff) {
00750 ptrdiff_t gsz = gridsize();
00751 for (ptrdiff_t i=0; i<gsz; i++){
00752 data[i] += ff;
00753 }
00754
00755
00756 cached_min += ff;
00757 cached_max += ff;
00758
00759 invalidate_mean();
00760 invalidate_sigma();
00761
00762
00763
00764 }
00765
00766
00767 void VolumetricData::rescale_voxel_value_range(float min_val, float max_val) {
00768 float newvrange = max_val - min_val;
00769 float oldvrange = cached_max - cached_min;
00770 float vrangeratio = newvrange / oldvrange;
00771
00772 ptrdiff_t gsz = gridsize();
00773 for (ptrdiff_t i=0; i<gsz; i++) {
00774 data[i] = min_val + vrangeratio*(data[i] - cached_min);
00775 }
00776
00777
00778 cached_min = min_val;
00779 cached_max = max_val;
00780
00781 invalidate_mean();
00782 invalidate_sigma();
00783
00784
00785 invalidate_gradient();
00786 }
00787
00788
00789
00790
00791
00792 void VolumetricData::downsample() {
00793 int xsize_new = xsize/2;
00794 int ysize_new = ysize/2;
00795 int zsize_new = zsize/2;
00796 float *data_new = new float[ptrdiff_t(xsize_new)*ptrdiff_t(ysize_new)*ptrdiff_t(zsize_new)];
00797
00798 int index_shift[8] = {0, 1, xsize, xsize+1, xsize*ysize, xsize*ysize + 1, xsize*ysize + xsize, xsize*ysize + xsize + 1};
00799
00800 int gx, gy, gz, j;
00801 const double eighth = 1.0/8.0;
00802 for (gz=0; gz<zsize_new; gz++) {
00803 for (gy=0; gy<ysize_new; gy++) {
00804 ptrdiff_t oldyzaddr = gy*ptrdiff_t(xsize) + gz*ptrdiff_t(xsize)*ptrdiff_t(ysize);
00805 ptrdiff_t newyzaddr = gy*ptrdiff_t(xsize_new) + gz*ptrdiff_t(xsize_new)*ptrdiff_t(ysize_new);
00806 for (gx=0; gx<xsize_new; gx++) {
00807 int n = 2*(gx + oldyzaddr);
00808 int n_new = gx + newyzaddr;
00809 double Z=0.0;
00810 for (j=0; j<8; j++)
00811 Z += data[n+index_shift[j]];
00812
00813 data_new[n_new] = float(Z * eighth);
00814 }
00815 }
00816 }
00817
00818 xsize = xsize_new;
00819 ysize = ysize_new;
00820 zsize = zsize_new;
00821
00822 double xscale = 0.5*(xsize)/(xsize/2);
00823 double yscale = 0.5*(ysize)/(ysize/2);
00824 double zscale = 0.5*(zsize)/(zsize/2);
00825 int i;
00826 for (i=0; i<3; i++) {
00827 xaxis[i] *= xscale;
00828 yaxis[i] *= yscale;
00829 zaxis[i] *= zscale;
00830 }
00831
00832 delete[] data;
00833 data = data_new;
00834
00835 invalidate_minmax();
00836 invalidate_mean();
00837 invalidate_sigma();
00838
00839
00840 invalidate_gradient();
00841 }
00842
00843
00844
00845 void VolumetricData::supersample() {
00846 int gx, gy, gz;
00847 int xsize_new = xsize*2-1;
00848 int ysize_new = ysize*2-1;
00849 int zsize_new = zsize*2-1;
00850 ptrdiff_t xysize = ptrdiff_t(xsize)*ptrdiff_t(ysize);
00851 ptrdiff_t xysize_new = ptrdiff_t(xsize_new)*ptrdiff_t(ysize_new);
00852 float *data_new = new float[ptrdiff_t(xsize_new)*ptrdiff_t(ysize_new)*ptrdiff_t(zsize_new)];
00853
00854
00855 for (gz=0; gz<zsize; gz++) {
00856 for (gy=0; gy<ysize; gy++) {
00857 ptrdiff_t oldyzplane = ptrdiff_t(gy)*ptrdiff_t(xsize) + ptrdiff_t(gz)*xysize;
00858 ptrdiff_t newyzplane = 2L*ptrdiff_t(gy)*ptrdiff_t(xsize_new) + 2L*ptrdiff_t(gz)*xysize_new;
00859
00860
00861
00862 for (gx=0; gx<xsize; gx++) {
00863 data_new[2*gx + newyzplane] = data[gx + oldyzplane];
00864 }
00865 }
00866 }
00867
00868
00869
00870
00871
00872 for (gz=0; gz<zsize; gz++) {
00873 for (gy=0; gy<ysize; gy++) {
00874 ptrdiff_t newyzplane = 2L*ptrdiff_t(gy)*ptrdiff_t(xsize_new) + 2L*ptrdiff_t(gz)*xysize_new;
00875 for (gx=1; gx<xsize-2; gx++) {
00876
00877
00878 data_new[2*gx+1 + newyzplane] =
00879 cubic_interp(data_new[(2*gx-2) + newyzplane],
00880 data_new[(2*gx) + newyzplane],
00881 data_new[(2*gx+2) + newyzplane],
00882 data_new[(2*gx+4) + newyzplane],
00883 0.5);
00884 }
00885 }
00886 }
00887
00888
00889 for (gz=0; gz<zsize; gz++) {
00890 for (gy=0; gy<ysize; gy++) {
00891 ptrdiff_t newyzplane = 2L*ptrdiff_t(gy)*ptrdiff_t(xsize_new) + 2L*ptrdiff_t(gz)*xysize_new;
00892
00893
00894
00895 data_new[1 + newyzplane] =
00896 cubic_interp(data_new[0 + newyzplane],
00897 data_new[0 + newyzplane],
00898 data_new[2 + newyzplane],
00899 data_new[4 + newyzplane],
00900 0.5);
00901
00902
00903 data_new[2*(xsize-2)+1 + newyzplane] =
00904 cubic_interp(data_new[2*(xsize-2)-2 + newyzplane],
00905 data_new[2*(xsize-2) + newyzplane],
00906 data_new[2*(xsize-2)+2 + newyzplane],
00907 data_new[2*(xsize-2)+2 + newyzplane],
00908 0.5);
00909 }
00910 }
00911
00912
00913 for (gz=0; gz<zsize; gz++) {
00914 for (gy=1; gy<ysize-2; gy++) {
00915 ptrdiff_t newzplane = 2L*ptrdiff_t(gz)*xysize_new;
00916 for (gx=0; gx<xsize_new; gx++) {
00917 data_new[gx + (2*gy+1)*xsize_new + 2*gz*xysize_new] =
00918 cubic_interp(data_new[gx + (2*gy-2)*ptrdiff_t(xsize_new) + newzplane],
00919 data_new[gx + (2*gy )*ptrdiff_t(xsize_new) + newzplane],
00920 data_new[gx + (2*gy+2)*ptrdiff_t(xsize_new) + newzplane],
00921 data_new[gx + (2*gy+4)*ptrdiff_t(xsize_new) + newzplane],
00922 0.5);
00923 }
00924 }
00925 }
00926
00927
00928 for (gz=0; gz<zsize; gz++) {
00929 ptrdiff_t newzplane = 2L*ptrdiff_t(gz)*xysize_new;
00930 for (gx=0; gx<xsize_new; gx++) {
00931
00932 data_new[gx + 1*xsize_new + newzplane] = \
00933 cubic_interp(data_new[gx + 0*xsize_new + newzplane],
00934 data_new[gx + 0*xsize_new + newzplane],
00935 data_new[gx + 2*xsize_new + newzplane],
00936 data_new[gx + 4*xsize_new + newzplane],
00937 0.5);
00938
00939
00940 data_new[gx + (2*(ysize-2)+1)*xsize_new + 2*gz*xysize_new] = \
00941 cubic_interp(data_new[gx + (2*(ysize-2)-2)*xsize_new + newzplane],
00942 data_new[gx + 2*(ysize-2)*xsize_new + newzplane],
00943 data_new[gx + (2*(ysize-2)+2)*xsize_new + newzplane],
00944 data_new[gx + (2*(ysize-2)+2)*xsize_new + newzplane],
00945 0.5);
00946 }
00947 }
00948
00949
00950 for (gy=0; gy<ysize_new; gy++) {
00951 for (gz=1; gz<zsize-2; gz++) {
00952 ptrdiff_t newyzplane = gy*ptrdiff_t(xsize_new) + (2L*gz+1L)*ptrdiff_t(xysize_new);
00953 for (gx=0; gx<xsize_new; gx++) {
00954 ptrdiff_t newxplusyrow = gx + gy*ptrdiff_t(xsize_new);
00955 data_new[gx + newyzplane] =
00956 cubic_interp(data_new[newxplusyrow + (2*gz-2)*xysize_new],
00957 data_new[newxplusyrow + (2*gz )*xysize_new],
00958 data_new[newxplusyrow + (2*gz+2)*xysize_new],
00959 data_new[newxplusyrow + (2*gz+4)*xysize_new],
00960 0.5);
00961 }
00962 }
00963 }
00964
00965
00966 for (gy=0; gy<ysize_new; gy++) {
00967 for (gx=0; gx<xsize_new; gx++) {
00968 ptrdiff_t newxplusyrow = gx + gy*ptrdiff_t(xsize_new);
00969
00970
00971 data_new[gx + gy*xsize_new + 1*xysize_new] = \
00972 cubic_interp(data_new[newxplusyrow + 0*xysize_new],
00973 data_new[newxplusyrow + 0*xysize_new],
00974 data_new[newxplusyrow + 2*xysize_new],
00975 data_new[newxplusyrow + 4*xysize_new],
00976 0.5);
00977
00978 data_new[gx + gy*xsize_new + (2*(zsize-2)+1)*xysize_new] = \
00979 cubic_interp(data_new[newxplusyrow + (2*(zsize-2)-2)*xysize_new],
00980 data_new[newxplusyrow + 2*(zsize-2)*xysize_new],
00981 data_new[newxplusyrow + (2*(zsize-2)+2)*xysize_new],
00982 data_new[newxplusyrow + (2*(zsize-2)+2)*xysize_new],
00983 0.5);
00984 }
00985 }
00986
00987 xsize = xsize_new;
00988 ysize = ysize_new;
00989 zsize = zsize_new;
00990
00991 delete[] data;
00992 data = data_new;
00993
00994
00995
00996
00997
00998 invalidate_minmax();
00999 invalidate_mean();
01000 invalidate_sigma();
01001
01002
01003 invalidate_gradient();
01004 }
01005
01006
01007
01008
01009 void VolumetricData::sigma_scale() {
01010 float oldmean = mean();
01011 float oldsigma = sigma();
01012 float oldsigma_inv = 1.0f / oldsigma;
01013
01014 ptrdiff_t gsz = gridsize();
01015 for (ptrdiff_t i=0; i<gsz; i++) {
01016 data[i] -= oldmean;
01017 data[i] *= oldsigma_inv;
01018 }
01019
01020 invalidate_minmax();
01021 invalidate_mean();
01022 invalidate_sigma();
01023
01024
01025 invalidate_gradient();
01026 }
01027
01028
01029
01030 void VolumetricData::binmask(float threshold) {
01031 ptrdiff_t gsz = gridsize();
01032 for (ptrdiff_t i=0; i<gsz; i++) {
01033 float tmp=data[i];
01034 data[i] = (tmp > threshold) ? 1.0f : 0.0f;
01035 }
01036
01037 invalidate_minmax();
01038 invalidate_mean();
01039 invalidate_sigma();
01040
01041
01042 invalidate_gradient();
01043 }
01044
01045
01046
01047
01048
01049
01050
01051 void VolumetricData::gaussian_blur(double sigma) {
01052
01053
01054
01055
01056
01057
01058 bool cuda = false;
01059 GaussianBlur<float>* gaussian_f = new GaussianBlur<float>(data, xsize, ysize, zsize, cuda);
01060 gaussian_f->blur(float(sigma));
01061
01062 delete[] data;
01063 data = gaussian_f->get_image();
01064
01065 invalidate_minmax();
01066 invalidate_mean();
01067 invalidate_sigma();
01068
01069
01070 invalidate_gradient();
01071 }
01072
01073 void VolumetricData::mdff_potential(float threshold) {
01074
01075 float threshinvert = -threshold;
01076
01077
01078 float mininvert = -cached_max;
01079
01080
01081 float oldvrange = threshinvert - mininvert;
01082 float vrangeratio = 1 / oldvrange;
01083
01084 ptrdiff_t gsz = gridsize();
01085 for (ptrdiff_t i=0; i<gsz; i++) {
01086
01087 if (data[i] < threshold)
01088 data[i] = threshold;
01089
01090
01091 data[i] = -data[i];
01092
01093
01094 data[i] = vrangeratio*(data[i] - mininvert);
01095 }
01096
01097
01098 cached_min = 0;
01099 cached_max = 1;
01100
01101 invalidate_mean();
01102 invalidate_sigma();
01103
01104
01105 invalidate_gradient();
01106 }
01107
01108
01109 void VolumetricData::invalidate_minmax() {
01110 minmax_isvalid = false;
01111 cached_min = 0;
01112 cached_max = 0;
01113 }
01114
01115
01116 void VolumetricData::datarange(float &min, float &max) {
01117 if (!minmax_isvalid && !mean_isvalid) {
01118 compute_minmaxmean();
01119 } else if (!minmax_isvalid) {
01120 compute_minmax();
01121 }
01122
01123 min = cached_min;
01124 max = cached_max;
01125 }
01126
01127
01128 void VolumetricData::compute_minmaxmean() {
01129 cached_min = 0;
01130 cached_max = 0;
01131 cached_mean = 0;
01132
01133 ptrdiff_t gsz = gridsize();
01134 if (gsz > 0) {
01135
01136 minmaxmean_1fv_aligned(data, gsz, &cached_min, &cached_max, &cached_mean);
01137 }
01138
01139 mean_isvalid = true;
01140 minmax_isvalid = true;
01141 }
01142
01143
01144 void VolumetricData::compute_minmax() {
01145 cached_min = 0;
01146 cached_max = 0;
01147
01148 ptrdiff_t gsz = gridsize();
01149 if (gsz > 0) {
01150
01151 minmax_1fv_aligned(data, gsz, &cached_min, &cached_max);
01152 }
01153
01154 minmax_isvalid = true;
01155 }
01156
01157
01158
01159 float VolumetricData::mean() {
01160 if (!mean_isvalid && !minmax_isvalid) {
01161 compute_minmaxmean();
01162 } else if (!mean_isvalid) {
01163 compute_mean();
01164 }
01165
01166 return cached_mean;
01167 }
01168
01169
01170 void VolumetricData::invalidate_mean() {
01171 mean_isvalid = false;
01172 cached_mean = 0;
01173 }
01174
01175
01176 float VolumetricData::sigma() {
01177 if (!sigma_isvalid)
01178 compute_sigma();
01179
01180 return cached_sigma;
01181 }
01182
01183
01184 void VolumetricData::invalidate_sigma() {
01185 sigma_isvalid = false;
01186 cached_sigma = 0;
01187 }
01188
01189
01190 double VolumetricData::integral() {
01191 double vol = cell_volume();
01192 return mean() * vol * xsize * ysize * zsize;
01193 }
01194
01195
01196 void VolumetricData::compute_mean() {
01197 double mean=0.0;
01198 ptrdiff_t gsz = gridsize();
01199
01200
01201
01202
01203
01204
01205 for (ptrdiff_t i=0; i<gsz; i++)
01206 mean += data[i];
01207 mean /= gsz;
01208 cached_mean = float(mean);
01209 mean_isvalid = true;
01210 }
01211
01212
01213 void VolumetricData::compute_sigma() {
01214 double sigma = 0.0;
01215 ptrdiff_t gsz = gridsize();
01216 float mymean = mean();
01217
01218
01219
01220
01221
01222
01223 for (ptrdiff_t i=0; i<gsz; i++) {
01224 float delta = data[i] - mymean;
01225 sigma += delta*delta;
01226 }
01227
01228 sigma /= gsz;
01229 sigma = sqrt(sigma);
01230 cached_sigma = float(sigma);
01231 sigma_isvalid = true;
01232 }
01233
01234