00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 00004 *cr University of Illinois 00005 *cr All Rights Reserved 00006 *cr 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * RCS INFORMATION: 00011 * 00012 * $RCSfile: VolumetricData.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.71 $ $Date: 2022年03月25日 05:50:31 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * 00019 * Base class for storing volumetric data and associated gradient data 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> // FLT_MAX etc 00032 00033 #ifndef NAN //not a number 00034 const float NAN = sqrtf(-1.f); //need some kind of portable NAN definition 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(); // min/max will be computed on next request 00053 invalidate_mean(); // mean will be computed on first request 00054 invalidate_sigma(); // sigma will be computed on first request 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(); // min/max will be computed on next request 00080 invalidate_mean(); // mean will be computed on first request 00081 invalidate_sigma(); // sigma will be computed on first request 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 // set scaling factors correctly 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 // set scaling factors correctly 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 // set scaling factors correctly 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 // shift coordinates by a half-voxel so that integer truncation works 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 //create transformation matrix... 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 // calculate grid coordinate, shifting by a half-voxel 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 // calculate grid coordinate, shifting by a half-voxel 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 // calculate grid coordinate, shifting by a half-voxel 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 // calculate cell side lengths 00457 cell_lengths(&xl, &yl, &zl); 00458 00459 // gradient axis scaling values and averaging factors, to correctly 00460 // calculate the gradient for volumes with irregular cell spacing 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 // Calculate the volume gradient at each grid cell. 00489 // Gradients are now stored unnormalized, since we need them in pure 00490 // form in order to draw field lines etc. Shading code will now have 00491 // to do renormalization for itself on-the-fly. 00492 00493 // XXX this gradient is only correct for orthogonal grids, since 00494 // we're using the array index offsets to calculate the gradient 00495 // rather than voxel coordinate offsets. This will have to be 00496 // re-worked for non-orthogonal datasets. 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 // Pad each side of the volmap's grid with zeros. 00600 // Negative padding results in cropping/trimming of the map 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 // Since an arbitrary part of the original map may have been 00649 // cropped out, we have to recompute the min/max voxel values 00650 // from scratch. 00651 // 00652 // If we know that the map was only padded and not cropped, 00653 // we could avoid this and instead just update datamin/datamax 00654 // for the new zero-valued voxels that have been added to the edges. 00655 invalidate_minmax(); // min/max will be computed on next request 00656 00657 // Both mean and sigma have to be recomputed from scratch when cropping. 00658 // They could be scaled in the case of padding since we know how 00659 // many new zero-valued voxels are added. 00660 // For now we always recompute them. 00661 invalidate_mean(); // mean will be computed on next request 00662 invalidate_sigma(); // sigma will be computed on next request 00663 00664 // force complete destruction, deallocation, allocation, and 00665 // recomputation of the gradient since the actual map dimensions 00666 // have changed and no longer match the old gradient 00667 invalidate_gradient(); 00668 } 00669 00670 00671 // Crop the map based on minmax values given in coordinate space. If 00672 // the 'cropping box' exceeds the map boundaries, the map is padded 00673 // with zeroes. 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 // Right now, only works for orthogonal cells. 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 // the pad() method will update datain/datamax for us 00689 pad(padxm, padxp, padym, padyp, padzm, padzp); 00690 } 00691 00692 00693 void VolumetricData::clamp(float min_value, float max_value) { 00694 // clamp the voxel values themselves 00695 ptrdiff_t gsz = gridsize(); 00696 for (ptrdiff_t i=0; i<gsz; i++) { 00697 #if 1 00698 // far more vectorizable, but always writes to memory 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 // slow arithmetic, but only writes to memory if clamping occurs 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 // update datamin/datamax as a result of the value clamping operation 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(); // mean will be computed on next request 00717 invalidate_sigma(); // sigma will be computed on next request 00718 00719 // force regeneration of volume gradient to match clamped voxels 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 // update min/max as a result of the value scaling operation 00731 if (ff >= 0.0) { 00732 cached_min *= ff; 00733 cached_max *= ff; 00734 } else { 00735 // max/min are swapped when negative scale factors are applied 00736 float tmp = cached_min; 00737 cached_min = cached_max * ff; 00738 cached_max = tmp * ff; 00739 } 00740 00741 invalidate_mean(); // mean will be computed on next request 00742 invalidate_sigma(); // sigma will be computed on next request 00743 00744 // force regeneration of volume gradient to match scaled voxels 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 // update datamin/datamax as a result of the scalar addition operation 00756 cached_min += ff; 00757 cached_max += ff; 00758 00759 invalidate_mean(); // mean will be computed on next request 00760 invalidate_sigma(); // sigma will be computed on next request 00761 00762 // adding a scalar has no impact on the volume gradient 00763 // so we leave the existing gradient in place, unmodified 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 // update min/max as a result of the rescaling operation 00778 cached_min = min_val; 00779 cached_max = max_val; 00780 00781 invalidate_mean(); // mean will be computed on next request 00782 invalidate_sigma(); // sigma will be computed on next request 00783 00784 // force regeneration of volume gradient to match scaled voxels 00785 invalidate_gradient(); 00786 } 00787 00788 00789 // Downsample grid size by factor of 2 00790 // Changed from volutil so it doesn't support PMF; unclear 00791 // whether that is still important. 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(); // min/max will be computed on next request 00836 invalidate_mean(); // mean will be computed on next request 00837 invalidate_sigma(); // sigma will be computed on next request 00838 00839 // force regeneration of volume gradient to match resampled voxels 00840 invalidate_gradient(); 00841 } 00842 00843 00844 // Supersample grid size by factor of 2 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 // Copy original voxels to the matching voxels in the new finer grid 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 // this only copies the matching (even) voxels over, the 00861 // odd ones must be added through interpolation later 00862 for (gx=0; gx<xsize; gx++) { 00863 data_new[2*gx + newyzplane] = data[gx + oldyzplane]; 00864 } 00865 } 00866 } 00867 00868 // Perform cubic interpolation for the rest of the (odd) voxels had not 00869 // yet been assigned above... 00870 00871 // x direction 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 // compute odd voxels through cubic interpolation 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 // borders 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 // compute odd voxels through cubic interpolation 00894 // gx = 0 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 // gx = xsize-2 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 // y direction 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 // borders 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 // gy = 0 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 // gy = ysize-2 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 // z direction 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 // borders 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 // gz = 0 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 // gz = zsize-2 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 // XXX in principle we might expect that supersampling a map 00995 // would have no impact on the min/max and only a small 00996 // affect on mean/sigma, but for paranoia's sake we will 00997 // recompute them all for the time being... 00998 invalidate_minmax(); // min/max will be computed on next request 00999 invalidate_mean(); // mean will be computed on next request 01000 invalidate_sigma(); // sigma will be computed on next request 01001 01002 // force regeneration of volume gradient to match resampled voxels 01003 invalidate_gradient(); 01004 } 01005 01006 01007 // Transform map to a sigma scale, so that isovalues in VMD correspond 01008 // to number of sigmas above the mean 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(); // min/max will be computed on next request 01021 invalidate_mean(); // mean will be computed on next request 01022 invalidate_sigma(); // sigma will be computed on next request 01023 01024 // force regeneration of volume gradient to match rescaled voxels 01025 invalidate_gradient(); 01026 } 01027 01028 01029 // Make a binary mask map for values above a threshold 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(); // min/max will be computed on next request 01038 invalidate_mean(); // mean will be computed on next request 01039 invalidate_sigma(); // sigma will be computed on next request 01040 01041 // force regeneration of volume gradient to match masked voxels 01042 invalidate_gradient(); 01043 } 01044 01045 01046 // 01047 // XXX will this work for an arbitrarily large sigma??? 01048 // XXX does it need to? 01049 // XXX needs to be fixed+tested for multi-billion voxel maps 01050 // 01051 void VolumetricData::gaussian_blur(double sigma) { 01052 /*#if defined(VMDCUDA) 01053 bool cuda = true; 01054 #else 01055 bool cuda = false; 01056 #endif 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(); // min/max will be computed on next request 01066 invalidate_mean(); // mean will be computed on next request 01067 invalidate_sigma(); // sigma will be computed on next request 01068 01069 // force regeneration of volume gradient to match blurred voxels 01070 invalidate_gradient(); 01071 } 01072 01073 void VolumetricData::mdff_potential(float threshold) { 01074 //calculate new max 01075 float threshinvert = -threshold; 01076 01077 //calculate new min 01078 float mininvert = -cached_max; 01079 01080 //range of the thresholded and inverted (scaleby -1) data 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 // clamp the voxel values themselves 01087 if (data[i] < threshold) 01088 data[i] = threshold; 01089 01090 //scaleby -1 01091 data[i] = -data[i]; 01092 01093 //rescale voxel value range between 0 and 1 01094 data[i] = vrangeratio*(data[i] - mininvert); 01095 } 01096 01097 // update min/max as a result of the rescaling operation 01098 cached_min = 0; 01099 cached_max = 1; 01100 01101 invalidate_mean(); // mean will be computed on next request 01102 invalidate_sigma(); // sigma will be computed on next request 01103 01104 // force regeneration of volume gradient to match scaled voxels 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(); // min/max/mean all need updating 01119 } else if (!minmax_isvalid) { 01120 compute_minmax(); // only min/max need updating 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 // use fast 16-byte-aligned min/max/mean routine 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 // use fast 16-byte-aligned min/max routine 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(); // min/max/mean all need updating 01162 } else if (!mean_isvalid) { 01163 compute_mean(); // only mean requires updating 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 // XXX This is a slow and lazy implementation that 01201 // loses precision if we get large magnitude 01202 // values early-on. 01203 // If there is a single NaN value amidst the data 01204 // the returned mean will be a NaN. 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 // XXX This is a slow and lazy implementation that 01219 // loses precision if we get large magnitude 01220 // values early-on. 01221 // If there is a single NaN value amidst the data 01222 // the returned mean will be a NaN. 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