Load data from MetaImage file, translate grayscale values to densities. More...

#include <users_modules/bone_remodelling/src/DensityMaps.hpp>

Collaboration diagram for BoneRemodeling::DataFromMetaIO:

Public Member Functions
	DataFromMetaIO (MetaImage &metaFile)

PetscErrorCode	fromOptions ()

PetscErrorCode	getDataForGiveIndex (const vector< int > &vec_idx, vector< double > &vec_data)

PetscErrorCode	getInidcesForGivenCoordinate (const double x, const double y, const double z, const double dx, const double dy, const double dz, vector< int > &vec_ix, vector< int > &vec_iy, vector< int > &vec_iz, vector< int > &vec_idx, vector< double > &vec_dist)

double	getAverage (const vector< double > &vec_data)
	Averages data within given vector. More...

PetscErrorCode	mapDensityLinear (double *paramDensity, const vector< double > &vec_data, vector< double > &vec_density)
	Transform data to densities. More...

PetscErrorCode	mapDensityLinear (const vector< double > &vec_data, vector< double > &vec_density)

double	getElasticMod (double *param, double rho)
	Transform densities to Youngs modulus. More...

double	getElasticMod (double rho)

double	medianFilter (const vector< double > &vec_density)
	function returns median of data from given vector More...

double	mapGaussSmooth (const double sigma, const vector< double > &vec_density, const vector< double > &vec_dist)
	function smooths values in a given vector depending on distance from center (vec_dist) More...

double	mapGaussSmooth (const vector< double > &vec_density, const vector< double > &vec_dist)

Public Attributes
MetaImage &	metaFile

double	sIgma

double	cUbe_size

int	fIlter

double	sCale

double	dIst_tol

double	tHreshold

double	paramElastic [2]

double	paramDensity [2]

double	lAmbda

vector< double >	kErnel

Detailed Description

Load data from MetaImage file, translate grayscale values to densities.

Definition at line 205 of file DensityMaps.hpp.

Constructor & Destructor Documentation

◆ DataFromMetaIO()

BoneRemodeling::DataFromMetaIO::DataFromMetaIO ( MetaImage & metaFile )

inline

Definition at line 208 of file DensityMaps.hpp.

                                       :
    metaFile(metaFile) {
    }

Member Function Documentation

◆ fromOptions()

PetscErrorCode BoneRemodeling::DataFromMetaIO::fromOptions ( )

inline

Definition at line 223 of file DensityMaps.hpp.

                                 {
      PetscFunctionBegin;
      ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Get parameters for DataFromMetaIO","none"); CHKERRQ(ierr);
 
      PetscBool flg = PETSC_TRUE;
 
      sIgma = 1;
      ierr = PetscOptionsScalar(
        "-density_map_sigma",
        "get sigma for gauss weighted average","",
        sIgma,&sIgma,PETSC_NULL
      ); CHKERRQ(ierr);
 
      cUbe_size = 0;
      ierr = PetscOptionsScalar(
        "-cube_size",
        "get cube size for the averaging","",
        0,&cUbe_size,PETSC_NULL
      ); CHKERRQ(ierr);
 
      dIst_tol = 0;
      ierr = PetscOptionsScalar(
        "-dist_tol",
        "get distance tolerance from mesh surface","",
        0,&dIst_tol,PETSC_NULL
      ); CHKERRQ(ierr);
 
      sCale = 1;
      ierr = PetscOptionsScalar(
        "-scale",
        "scale the distance between voxels (eg. from mm to m)","",
        sCale,&sCale,PETSC_NULL
      ); CHKERRQ(ierr);
 
      tHreshold = -1024;
      ierr = PetscOptionsScalar(
        "-threshold",
        "minimum value of ct threshold","",
        tHreshold,&tHreshold,PETSC_NULL
      ); CHKERRQ(ierr);
 
      fIlter = 0;
      ierr = PetscOptionsInt(
        "-filter",
        "number of the filter","",
        fIlter,&fIlter,PETSC_NULL
      ); CHKERRQ(ierr);
 
      lAmbda = 1;
      ierr = PetscOptionsScalar(
        "-lambda",
        "get lambda coefficient for smoothing","",
        lAmbda,&lAmbda,PETSC_NULL
      ); CHKERRQ(ierr);
 
      paramElastic[0] = 0;
      paramElastic[1] = 1;
      int nmax = 2;
      ierr = PetscOptionsGetScalarArray(
        PETSC_NULL,"-param_elastic",paramElastic,&nmax,&flg
      ); CHKERRQ(ierr);
 
      paramDensity[0] = 1;
      paramDensity[1] = 0;
      int n_max = 2;
      ierr = PetscOptionsGetScalarArray(PETSC_NULL,"-param_density",paramDensity,&n_max,&flg
  ); CHKERRQ(ierr);
 
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
 
  PetscFunctionReturn(0);
}

◆ getAverage()

double BoneRemodeling::DataFromMetaIO::getAverage ( const vector< double > & vec_data )

inline

Averages data within given vector.

Parameters

vector data from metaFile

Returns: averaged data

Definition at line 485 of file DensityMaps.hpp.

                                                      {
      PetscFunctionBegin;
      double data = 0;
      int i=0;
      for(vector<double>::const_iterator vec_it = vec_data.begin();vec_it!=vec_data.end();vec_it++) {
        if (*vec_it < tHreshold) continue;
        data += *vec_it;
        i++;
      }
      if (i == 0) i = 1;    // to avoid division by zero
      data /= i;
      //cout <<"data "<< data << " SIZE " << vec_data.size() << endl;
      return data;
    }

◆ getDataForGiveIndex()

PetscErrorCode BoneRemodeling::DataFromMetaIO::getDataForGiveIndex	(	const vector< int > &	vec_idx,
		vector< double > &	vec_data
	)

inline

returns vertors of data for given index of the point

Parameters

vec_idx	reference vector of indices
vec_data	reference vector of data (colors)

Returns: error code

Definition at line 302 of file DensityMaps.hpp.

      {
      PetscFunctionBegin;
      int size_of_vec_idx = vec_idx.size();
      vec_data.resize(size_of_vec_idx);
 
      int ii = 0;
      for(vector<int>::const_iterator vit = vec_idx.begin();vit!=vec_idx.end();vit++) {
        vec_data[ii] = metaFile.ElementData(*vit);
        // if (vec_data[ii] < 500) SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong value (lower than 1000)");
        // cout << "INDEX "<< *vit << " Data " << vec_data[ii] << endl;
        ii++;
      } //cout << "-----------------------" << endl;
      PetscFunctionReturn(0);
    }

◆ getElasticMod() [1/2]

double BoneRemodeling::DataFromMetaIO::getElasticMod	(	double *	param,
		double	rho
	)

inline

Transform densities to Youngs modulus.

Parameters

a,b	parameter for relation of density - density, b * density ^ a = E
rho	density

Returns: youngs modulus

Parameters

param	parameters a,b for density-elasticity relation
rho	density

Definition at line 538 of file DensityMaps.hpp.

       {
       return param[1] * pow(rho , param[0]); //         ///< power density mapping: param_a and param_b
     }

◆ getElasticMod() [2/2]

double BoneRemodeling::DataFromMetaIO::getElasticMod ( double rho )

inline

Definition at line 545 of file DensityMaps.hpp.

                                      {
       return getElasticMod(paramElastic,rho); //
     }

◆ getInidcesForGivenCoordinate()

PetscErrorCode BoneRemodeling::DataFromMetaIO::getInidcesForGivenCoordinate	(	const double	x,
		const double	y,
		const double	z,
		const double	dx,
		const double	dy,
		const double	dz,
		vector< int > &	vec_ix,
		vector< int > &	vec_iy,
		vector< int > &	vec_iz,
		vector< int > &	vec_idx,
		vector< double > &	vec_dist
	)

inline

returns vertors of coordinates and indices for given coordinate

Parameters

x	coordinate of point x
y	coordinate of point y
z	coordinate of point z
dx	dimension of the block in x direction
dy	dimension of the block in y direction
dz	dimension of the block in z direction
vec_ix	reference vector of x coordinates
vec_iy	reference vector of y coordinates
vec_iz	reference vector of z coordinates
vec_idx	reference vector of indices

Returns: error code

Parameters

z	center of the cube
dz	size of the cube
vec_ix	vectors of coord points inside cube
vec_idx	indices of inside cube with dimensions dx dy dz
vec_dist	vector of squared distances between center of the cube and points

Definition at line 334 of file DensityMaps.hpp.

      {
      PetscFunctionBegin;
 
      const double *elemet_size = metaFile.ElementSize();
      const int *dim_size = metaFile.DimSize();
 
      int nx = ceil(dx/ (elemet_size[0] * sCale));                   // number of points within given size of block - dx
      int ny = ceil(dy/ (elemet_size[1] * sCale));                   // number of points within given size of block - dy
      int nz = ceil(dz/ (elemet_size[2] * sCale));                   // number of points within given size of block - dz
 
      //const double *offset = metaFile.Offset();           // coords before translation (offset)
      double x_befor_offset = x;
      double y_befor_offset = y;
      double z_befor_offset = z;
      // double x_befor_offset = x - metaFile.Offset()[0];  // probably its not necessary
      // double y_befor_offset = y - metaFile.Offset()[1];
      // double z_befor_offset = z - metaFile.Offset()[2];
      // center of the cube
      int ix  = ceil(x_befor_offset/ (elemet_size[0] * sCale));     // original coords/indices (before multiplying by elem size)
      int iy  = ceil(y_befor_offset/ (elemet_size[1] * sCale));
      int iz  = ceil(z_befor_offset/ (elemet_size[2] * sCale));
 
      vec_ix.resize(2*nx);                                // setting sizes of the vectors equal to number of points
      vec_iy.resize(2*ny);
      vec_iz.resize(2*nz);
 
      int ii = 0;
      for(int i = 0;i<2*nx;i++) {                         // vectors of coordinates of points within given cube
        int idx = ix-nx+i;
        if(idx >= dim_size[0] || idx < 0) continue;       // erasing indices beyond the border
        vec_ix[ii++] = idx;
      }
      vec_ix.resize(ii);
      ii = 0;
      for(int i = 0;i<2*ny;i++) {
        int idx = iy-ny+i;
        if(idx >= dim_size[1] || idx < 0) continue;       // erasing indices beyond the border
        vec_iy[ii++] = idx;
      }
      vec_iy.resize(ii);
      ii = 0;
      for(int i = 0;i<2*nz;i++) {
        int idx = iz-nz+i;
        if(idx >= dim_size[2] || idx < 0) continue;       // erasing indices beyond the border
        vec_iz[ii++] = idx;
      }
      vec_iz.resize(ii);
 
 
      if (dx <= 0) {     //
        vec_ix.resize(1);
        vec_ix.push_back(round(x_befor_offset/ (elemet_size[0] * sCale)));
      }
      if (dy <= 0) {
        vec_iy.resize(1);
        vec_iy.push_back(round(y_befor_offset/ (elemet_size[1] * sCale)));
      }
      if (dz <= 0) {
        vec_iz.resize(1);
        vec_iz.push_back(round(z_befor_offset/ (elemet_size[2] * sCale)));
      }
 
      // vec_idx.resize(8*nx*ny*nz);
      // vec_dist.resize(8*nx*ny*nz);
      vec_idx.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
      vec_dist.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
      int dm_size0_dim_size1 = dim_size[0] * dim_size[1];
      // calculating the vector of global indices and vector of distances from center of the cube
 
      if (fIlter == 0) {
 
        // vector<double> dy2(vec_iy.size());
        // {
        //   int ii = 0;
        //   for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
        //     double dyy = (*it_iy*elemet_size[1] - y_befor_offset);
        //     dy2[ii++] =  dyy * dyy;
        //   }
        // }
        // vector<double> dz2(vec_iz.size());
        // {
        //   int ii = 0;
        //   for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
        //     double dzz = (*it_iz*elemet_size[2] - z_befor_offset);
        //     dz2[ii++] = dzz * dzz;
        //   }
        // }
 
        ii = 0;
        for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
          double a = (*it_iz) * dm_size0_dim_size1 ;
          double dz2 = (*it_iz*(elemet_size[2]*sCale) - z_befor_offset) * (*it_iz*(elemet_size[2]*sCale) - z_befor_offset);
          for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
            double b = (*it_iy) * dim_size[0];
            double dy2 = (*it_iy*(elemet_size[1]*sCale) - y_befor_offset) * (*it_iy*(elemet_size[1]*sCale) - y_befor_offset);
            for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
              vec_idx[ii] =  a + b + (*it_ix);
              vec_dist[ii] =  ( dz2 + dy2 + (*it_ix*(elemet_size[0]*sCale) - x_befor_offset)
               * (*it_ix*(elemet_size[0]*sCale) - x_befor_offset) );
              ii++;
              // cout << "Ix Iy Iz" << vec_idx[ii] << " " << (*it_ix) << " " << (*it_iy) << " " << (*it_iz) << endl;
              // cout << "Z " << *it_iz << "-" <<  vec_idx[ii] / (dim_size[0] * dim_size[1]) << " ";
              // cout << "Y " << *it_iy << "-" << (vec_idx[ii] - (*it_iz) * (dim_size[0] * dim_size[1])) / dim_size[0] << " ";
              // cout << "X " << *it_ix << "-" << vec_idx[ii] - (*it_iz) * dim_size[0] * dim_size[1] - (*it_iy) * dim_size[0] << " " << endl;
              // cout << "POINT "<< ii << " GLOBAL  " << vec_idx[ii] <<
              // " DATA "<< metaFile.ElementData(vec_idx[ii]) << " DISTANCE  " << vec_dist[ii] << endl;
 
 
              // if we require more indices than mhd file has ----> ERROR
              // if(vec_idx[ii] >= metaFile.Quantity()) {
              //   SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong index of voxels (> Quantity)");
              // }
 
 
            }
          }
        }
 
      }
      else {
 
        ii = 0;
        for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
          double a = (*it_iz) * dm_size0_dim_size1 ;
          for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
            double b = (*it_iy) * dim_size[0];
            for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
              vec_idx[ii++] =  a + b + (*it_ix);
            }
          }
        }
      }
      // cout << ix << " " << " " << iy << " " << iz <<endl;
      //cout << "--------------------------------------------------------" << endl;
      PetscFunctionReturn(0);
    }

◆ mapDensityLinear() [1/2]

PetscErrorCode BoneRemodeling::DataFromMetaIO::mapDensityLinear	(	const vector< double > &	vec_data,
		vector< double > &	vec_density
	)

inline

Definition at line 523 of file DensityMaps.hpp.

      {
      PetscFunctionBegin;
      ierr = mapDensityLinear(paramDensity,vec_data,vec_density); CHKERRQ(ierr);
      PetscFunctionReturn(0);
    }

◆ mapDensityLinear() [2/2]

PetscErrorCode BoneRemodeling::DataFromMetaIO::mapDensityLinear	(	double *	paramDensity,
		const vector< double > &	vec_data,
		vector< double > &	vec_density
	)

inline

Transform data to densities.

Parameters

a	array of parameters a,b for linear relation of grayscale - density, a * HU + b = density
vec_data	data from metaFile
vec_density	output vector with density

Returns: vector of densities

Definition at line 508 of file DensityMaps.hpp.

      {
      PetscFunctionBegin;
      vec_density.resize(vec_data.size());
      int ii = 0;
      for(vector<double>::const_iterator vec_it = vec_data.begin();vec_it!=vec_data.end();vec_it++) {
        vec_density[ii++] = fmax(paramDensity[0] * (*vec_it) + paramDensity[1],0); //       // linear density mapping: param_a and param_b
        //cout << "density " << fmax(a * (*vec_it) + b,0) << endl;
        //cerr << "tHreshold :" << tHreshold << " fIlter: " << fIlter << " dist toler" << dIst_tol << endl;
      }
      PetscFunctionReturn(0);
    }

◆ mapGaussSmooth() [1/2]

double BoneRemodeling::DataFromMetaIO::mapGaussSmooth	(	const double	sigma,
		const vector< double > &	vec_density,
		const vector< double > &	vec_dist
	)

inline

function smooths values in a given vector depending on distance from center (vec_dist)

< normalisation of kernel and smooting values

< avoid division by zero (ad exception) this should not happen /FIXME

Parameters

sigma	sigma parameter for gausssian smoothing
vec_density	vector with non-smoothed densities
vec_dist	vector with squared distances form center of the cubes

Definition at line 581 of file DensityMaps.hpp.

      {
      //vec_density_smooth.resize(vec_density.size());
      kErnel.resize(vec_density.size());
      //calculating kernel
      int i = 0;
      double sum = 0;
      const double k = 2 * sigma*sigma;
      const double m = (1 / sqrt( M_PI * k));
      for(int ii = 0; ii < vec_dist.size(); ii++) {
        if (vec_density[ii] < tHreshold) continue;
        kErnel[i] = m * exp(- (vec_dist[ii]) / k); // distance is already squared
        sum += kErnel[i++];
        //cout <<"KERNEL"<< kernel[i] << endl;;
      }
      kErnel.resize(i);
      ///< normalisation of kernel and smooting values
      int ii = 0;
      double smooth_data = 0;
      for(vector<double>::const_iterator vec_itr = vec_density.begin();vec_itr!=vec_density.end();vec_itr++) {
        if (*vec_itr < tHreshold) continue;
        kErnel[ii] /= sum;
        //cout << kernel[ii] << endl;
        smooth_data += (*vec_itr) * kErnel[ii++];
      }
      if (sum == 0) return 0;    ///< avoid division by zero  (ad exception) this should not happen /FIXME
      return smooth_data;
    }

◆ mapGaussSmooth() [2/2]

double BoneRemodeling::DataFromMetaIO::mapGaussSmooth	(	const vector< double > &	vec_density,
		const vector< double > &	vec_dist
	)

inline

Parameters

vec_density	vector with non-smoothed densities
vec_dist	vector with distances form center of the cubes

Definition at line 614 of file DensityMaps.hpp.

      {
      return mapGaussSmooth(sIgma,vec_density,vec_dist);
    }

◆ medianFilter()

double BoneRemodeling::DataFromMetaIO::medianFilter ( const vector< double > & vec_density )

inline

function returns median of data from given vector

Parameters

vec_density vector with densities

Returns: median value

< copy of the vector to allow sorting data

< applying threshold

< sorting vector

Definition at line 553 of file DensityMaps.hpp.

                                                           {
      vector<double> vec_copy;                                  ///< copy of the vector to allow sorting data
      vec_copy.resize(vec_density.size());
      int i = 0;
      for(vector<double>::const_iterator vec_it = vec_density.begin();vec_it!=vec_density.end();vec_it++){
        if (*vec_it < tHreshold) continue;                     ///< applying threshold
        vec_copy[i++] = *vec_it;
        //cout << *vec_it << "-";
      }
      vec_copy.resize(i);
 
      //cout << " " << endl;
      sort(vec_copy.begin(), vec_copy.end());                  ///< sorting vector
      //cout << " " << endl;
      if (vec_copy.size() % 2 == 0) {
        return (vec_copy[vec_copy.size() / 2 - 1] + vec_copy[vec_copy.size() / 2]) / 2.0;
      } else {
        return vec_copy[vec_copy.size() / 2];
      }
      return 0;
    }

Member Data Documentation

◆ cUbe_size

double BoneRemodeling::DataFromMetaIO::cUbe_size

Definition at line 214 of file DensityMaps.hpp.

◆ dIst_tol

double BoneRemodeling::DataFromMetaIO::dIst_tol

Definition at line 217 of file DensityMaps.hpp.

◆ fIlter

int BoneRemodeling::DataFromMetaIO::fIlter

Definition at line 215 of file DensityMaps.hpp.

◆ kErnel

vector<double> BoneRemodeling::DataFromMetaIO::kErnel

Definition at line 575 of file DensityMaps.hpp.

◆ lAmbda

double BoneRemodeling::DataFromMetaIO::lAmbda

Definition at line 221 of file DensityMaps.hpp.

◆ metaFile

MetaImage& BoneRemodeling::DataFromMetaIO::metaFile

Definition at line 207 of file DensityMaps.hpp.

◆ paramDensity

double BoneRemodeling::DataFromMetaIO::paramDensity[2]

Definition at line 220 of file DensityMaps.hpp.

◆ paramElastic

double BoneRemodeling::DataFromMetaIO::paramElastic[2]

Definition at line 219 of file DensityMaps.hpp.

◆ sCale

double BoneRemodeling::DataFromMetaIO::sCale

Definition at line 216 of file DensityMaps.hpp.

◆ sIgma

double BoneRemodeling::DataFromMetaIO::sIgma

Definition at line 213 of file DensityMaps.hpp.

◆ tHreshold

double BoneRemodeling::DataFromMetaIO::tHreshold

Definition at line 218 of file DensityMaps.hpp.

The documentation for this struct was generated from the following file:

users_modules/bone_remodelling/src/DensityMaps.hpp

Public Member Functions

Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ DataFromMetaIO()

Member Function Documentation

◆ fromOptions()

◆ getAverage()

◆ getDataForGiveIndex()

◆ getElasticMod() [1/2]

◆ getElasticMod() [2/2]

◆ getInidcesForGivenCoordinate()

◆ mapDensityLinear() [1/2]

◆ mapDensityLinear() [2/2]

◆ mapGaussSmooth() [1/2]

◆ mapGaussSmooth() [2/2]

◆ medianFilter()

Member Data Documentation

◆ cUbe_size

◆ dIst_tol

◆ fIlter

◆ kErnel

◆ lAmbda

◆ metaFile

◆ paramDensity

◆ paramElastic

◆ sCale

◆ sIgma

◆ tHreshold