---

VolumetricData.h

Go to the documentation of this file.

/***************************************************************************
 *cr                                                                       
 *cr            (C) Copyright 1995-2019 The Board of Trustees of the           
 *cr                        University of Illinois                       
 *cr                         All Rights Reserved                        
 *cr                                                                   
 ***************************************************************************/

/***************************************************************************
 * RCS INFORMATION:
 *
 *      $RCSfile: VolumetricData.h,v $
 *      $Author: johns $        $Locker:  $             $State: Exp $
 *      $Revision: 1.57 $       $Date: 2020/10/22 03:43:24 $
 *
 ***************************************************************************/
#ifndef VOLUMETRICDATA_H
#define VOLUMETRICDATA_H

#include <stddef.h>

class VolumetricData {
public:
  double origin[3];        
  double xaxis[3];         
  double yaxis[3];         
  double zaxis[3];         
  int xsize, ysize, zsize; 
  char *name;              
  float *data;             


  VolumetricData(const char *name, const float *origin, 
                 const float *xaxis, const float *yaxis, const float *zaxis,
                 int xs, int ys, int zs, float *dataptr);

  VolumetricData(const char *name, const double *origin, 
                 const double *xaxis, const double *yaxis, const double *zaxis,
                 int xs, int ys, int zs, float *dataptr);

  ~VolumetricData();

  ptrdiff_t gridsize() const { return ptrdiff_t(xsize)*ptrdiff_t(ysize)*ptrdiff_t(zsize); }

  void datarange(float &min, float &max);

  void set_name(const char* name);


  float mean();

  float sigma();

  double integral();

  void cell_lengths(float *xl, float *yl, float *zl) const;

  void cell_axes(float *xax, float *yax, float *zax) const;
  void cell_axes(double *xax, double *yax, double *zax) const;

  void cell_dirs(float *xax, float *yax, float *zax) const;

  double cell_volume() const;

  void set_volume_axes(float *xax, float *yax, float *zax) {
    for (int idx = 0; idx < 3; idx++) {
      xaxis[idx] = xax[idx];
      yaxis[idx] = yax[idx];
      zaxis[idx] = zax[idx];
    }

    // force complete destruction, deallocation, allocation, and
    // recomputation of the gradient since the map voxel scale
    // has changed and no longer matches the old gradient
    // XXX we might be able to simply apply the correct scale
    //     to the existing gradient map, but its simple to just invalidate
    invalidate_gradient();
  }

  void set_volume_axes(double *xax, double *yax, double *zax) {
    for (int idx = 0; idx < 3; idx++) {
      xaxis[idx] = xax[idx];
      yaxis[idx] = yax[idx];
      zaxis[idx] = zax[idx];
    }

    // force complete destruction, deallocation, allocation, and
    // recomputation of the gradient since the map voxel scale
    // has changed and no longer matches the old gradient
    // XXX we might be able to simply apply the correct scale
    //     to the existing gradient map, but its simple to just invalidate
    invalidate_gradient();
  }

  void scale_volume(double scalex, double scaley, double scalez) {
    for (int idx = 0; idx < 3; idx++) {
      xaxis[idx] *= scalex;
      yaxis[idx] *= scaley;
      zaxis[idx] *= scalez;
    }

    // force complete destruction, deallocation, allocation, and
    // recomputation of the gradient since the map voxel scale
    // has changed and no longer matches the old gradient
    // XXX we might be able to simply apply the correct scale
    //     to the existing gradient map, but its simple to just invalidate
    invalidate_gradient();
  }

  void set_volume_origin(float *org) {
    for (int idx = 0; idx < 3; idx++) {
      origin[idx] = org[idx];
    }
    // no cached stats/values to invalidate at present...
  }

  void set_volume_origin(double *org) {
    for (int idx = 0; idx < 3; idx++) {
      origin[idx] = org[idx];
    }
    // no cached stats/values to invalidate at present...
  }


  void voxel_coord_from_cartesian_coord(const float *carcoord, float *voxcoord, int shiftflag) const;

  ptrdiff_t voxel_index_from_coord(float xpos, float ypos, float zpos) const;

  inline float voxel_value(int x, int y, int z) const {
    return data[z*ptrdiff_t(xsize)*ptrdiff_t(ysize) + y*ptrdiff_t(xsize) + x];
  }

  float voxel_value_safe(int x, int y, int z) const;

  float voxel_value_interpolate(float xv, float yv, float zv) const;

  float voxel_value_from_coord(float xpos, float ypos, float zpos) const;
  float voxel_value_interpolate_from_coord(float xpos, float ypos, float zpos) const;
  float voxel_value_from_coord_safe(float xpos, float ypos, float zpos) const;
  float voxel_value_interpolate_from_coord_safe(float xpos, float ypos, float zpos) const;


  const float *access_volume_gradient();

  void set_volume_gradient(float *gradient);

  void compute_volume_gradient(void);


  void voxel_gradient_fast(int x, int y, int z, float *grad) const {
    ptrdiff_t index = (z*ptrdiff_t(xsize)*ptrdiff_t(ysize) + y*ptrdiff_t(xsize) + x) * 3L;
    grad[0] = gradient[index    ];
    grad[1] = gradient[index + 1];
    grad[2] = gradient[index + 2];
  }

  void voxel_gradient_safe(int x, int y, int z, float *grad) const;

  void voxel_gradient_interpolate(const float *voxcoord, float *gradient) const;

  void voxel_gradient_from_coord(const float *coord, float *gradient) const;
  void voxel_gradient_interpolate_from_coord(const float *coord, float *gradient) const;

  inline void voxel_coord(int x, int y, int z, 
                          float &gx, float &gy, float &gz) const {
    float xdelta[3], ydelta[3], zdelta[3];
    cell_axes(xdelta, ydelta, zdelta);
    
    gx = float(origin[0] + (x * xdelta[0]) + (y * ydelta[0]) + (z * zdelta[0]));
    gy = float(origin[1] + (x * xdelta[1]) + (y * ydelta[1]) + (z * zdelta[1]));
    gz = float(origin[2] + (x * xdelta[2]) + (y * ydelta[2]) + (z * zdelta[2]));
  }

  inline void voxel_coord(ptrdiff_t i, float &x, float &y, float &z) const {
    float xdelta[3], ydelta[3], zdelta[3];
    cell_axes(xdelta, ydelta, zdelta);
    
    ptrdiff_t gz = i / (ptrdiff_t(ysize)*ptrdiff_t(xsize));
    ptrdiff_t gy = (i / xsize) % ysize;
    ptrdiff_t gx = i % xsize;

    x = float(origin[0] + (gx * xdelta[0]) + (gy * ydelta[0]) + (gz * zdelta[0]));
    y = float(origin[1] + (gx * xdelta[1]) + (gy * ydelta[1]) + (gz * zdelta[1]));
    z = float(origin[2] + (gx * xdelta[2]) + (gy * ydelta[2]) + (gz * zdelta[2]));
  }

  //
  //--unary ops--
  // 
  
  void pad(int padxm, int padxp, int padym, int padyp, int padzm, int padzp);
  
  void crop(double crop_minx, double crop_miny, double crop_minz, double crop_maxx, double crop_maxy, double crop_maxz);

  void clamp(float min_value, float max_value);

  void scale_by(float ff);

  void scalar_add(float ff);

  void rescale_voxel_value_range(float min_value, float max_value);

  void downsample();

  void supersample();

  void sigma_scale();

  void binmask(float threshold=0.0f);
  
  void gaussian_blur(double sigma);
  
  void mdff_potential(float threshold);

private:
  float *gradient;            
  bool gradient_isvalid;      
  void invalidate_gradient(); 

  bool minmax_isvalid;        
  float cached_min;           
  float cached_max;           

  bool mean_isvalid;          
  float cached_mean;          
 
  bool sigma_isvalid;         
  float cached_sigma;         

  void compute_minmaxmean();  
  void invalidate_minmax();   
  void compute_minmax();      

  void invalidate_mean();     
  void compute_mean();        

  void invalidate_sigma();    
  void compute_sigma();       


  inline float cubic_interp(float y0, float y1, float y2, float y3, float mu) const {
    float mu2 = mu*mu;
    float a0 = y3 - y2 - y0 + y1;
    float a1 = y0 - y1 - a0;
    float a2 = y2 - y0;
    float a3 = y1;

    return (a0*mu*mu2+a1*mu2+a2*mu+a3);
  }
  
};


//
// Fast and loose accessor macros, don't use unless you have to
// 

#define VOXEL_GRADIENT_FAST_IDX(gradientmap, index, newgrad) \
  { (newgrad)[0] = (gradientmap)[index    ]; \
    (newgrad)[1] = (gradientmap)[index + 1]; \
    (newgrad)[2] = (gradientmap)[index + 2]; \
  }

#define VOXEL_GRADIENT_FAST(gradientmap, planesz, rowsz, x, y, z, newgrad) \
  { ptrdiff_t index = ((z)*(planesz) + (y)*(rowsz) + (x)) * 3L; \
    (newgrad)[0] = (gradientmap)[index    ]; \
    (newgrad)[1] = (gradientmap)[index + 1]; \
    (newgrad)[2] = (gradientmap)[index + 2]; \
  }

#endif // VOLUMETRICDATA_H

---