Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

VolumetricData.C

Go to the documentation of this file.
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 

Generated on Tue Dec 3 02:46:14 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002