v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
BoneRemodeling::DataFromMetaIO Struct Reference

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

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

Collaboration diagram for BoneRemodeling::DataFromMetaIO:
[legend]

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< doublekErnel
 

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.

208 :
210 }

Member Function Documentation

◆ fromOptions()

PetscErrorCode BoneRemodeling::DataFromMetaIO::fromOptions ( )
inline

Definition at line 223 of file DensityMaps.hpp.

223 {
224 PetscFunctionBegin;
225 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Get parameters for DataFromMetaIO","none"); CHKERRQ(ierr);
226
227 PetscBool flg = PETSC_TRUE;
228
229 sIgma = 1;
230 ierr = PetscOptionsScalar(
231 "-density_map_sigma",
232 "get sigma for gauss weighted average","",
233 sIgma,&sIgma,PETSC_NULL
234 ); CHKERRQ(ierr);
235
236 cUbe_size = 0;
237 ierr = PetscOptionsScalar(
238 "-cube_size",
239 "get cube size for the averaging","",
240 0,&cUbe_size,PETSC_NULL
241 ); CHKERRQ(ierr);
242
243 dIst_tol = 0;
244 ierr = PetscOptionsScalar(
245 "-dist_tol",
246 "get distance tolerance from mesh surface","",
247 0,&dIst_tol,PETSC_NULL
248 ); CHKERRQ(ierr);
249
250 sCale = 1;
251 ierr = PetscOptionsScalar(
252 "-scale",
253 "scale the distance between voxels (eg. from mm to m)","",
254 sCale,&sCale,PETSC_NULL
255 ); CHKERRQ(ierr);
256
257 tHreshold = -1024;
258 ierr = PetscOptionsScalar(
259 "-threshold",
260 "minimum value of ct threshold","",
261 tHreshold,&tHreshold,PETSC_NULL
262 ); CHKERRQ(ierr);
263
264 fIlter = 0;
265 ierr = PetscOptionsInt(
266 "-filter",
267 "number of the filter","",
268 fIlter,&fIlter,PETSC_NULL
269 ); CHKERRQ(ierr);
270
271 lAmbda = 1;
272 ierr = PetscOptionsScalar(
273 "-lambda",
274 "get lambda coefficient for smoothing","",
275 lAmbda,&lAmbda,PETSC_NULL
276 ); CHKERRQ(ierr);
277
278 paramElastic[0] = 0;
279 paramElastic[1] = 1;
280 int nmax = 2;
282 PETSC_NULL,"-param_elastic",paramElastic,&nmax,&flg
283 ); CHKERRQ(ierr);
284
285 paramDensity[0] = 1;
286 paramDensity[1] = 0;
287 int n_max = 2;
288 ierr = PetscOptionsGetScalarArray(PETSC_NULL,"-param_density",paramDensity,&n_max,&flg
289 ); CHKERRQ(ierr);
290
291 ierr = PetscOptionsEnd(); CHKERRQ(ierr);
292
293 PetscFunctionReturn(0);
294}
static PetscErrorCode ierr
PetscErrorCode PetscOptionsGetScalarArray(PetscOptions options, const char pre[], const char name[], PetscScalar dvalue[], PetscInt *nmax, PetscBool *set)

◆ getAverage()

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

Averages data within given vector.

Parameters
vectordata from metaFile
Returns
averaged data

Definition at line 485 of file DensityMaps.hpp.

485 {
486 PetscFunctionBegin;
487 double data = 0;
488 int i=0;
489 for(vector<double>::const_iterator vec_it = vec_data.begin();vec_it!=vec_data.end();vec_it++) {
490 if (*vec_it < tHreshold) continue;
491 data += *vec_it;
492 i++;
493 }
494 if (i == 0) i = 1; // to avoid division by zero
495 data /= i;
496 //cout <<"data "<< data << " SIZE " << vec_data.size() << endl;
497 return data;
498 }
FTensor::Index< 'i', SPACE_DIM > i

◆ 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_idxreference vector of indices
vec_datareference vector of data (colors)
Returns
error code

Definition at line 302 of file DensityMaps.hpp.

305 {
306 PetscFunctionBegin;
307 int size_of_vec_idx = vec_idx.size();
308 vec_data.resize(size_of_vec_idx);
309
310 int ii = 0;
311 for(vector<int>::const_iterator vit = vec_idx.begin();vit!=vec_idx.end();vit++) {
312 vec_data[ii] = metaFile.ElementData(*vit);
313 // if (vec_data[ii] < 500) SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong value (lower than 1000)");
314 // cout << "INDEX "<< *vit << " Data " << vec_data[ii] << endl;
315 ii++;
316 } //cout << "-----------------------" << endl;
317 PetscFunctionReturn(0);
318 }

◆ getElasticMod() [1/2]

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

Transform densities to Youngs modulus.

Parameters
a,bparameter for relation of density - density, b * density ^ a = E
rhodensity
Returns
youngs modulus
Parameters
paramparameters a,b for density-elasticity relation
rhodensity

Definition at line 538 of file DensityMaps.hpp.

541 {
542 return param[1] * pow(rho , param[0]); // ///< power density mapping: param_a and param_b
543 }
double rho
Definition: plastic.cpp:191

◆ getElasticMod() [2/2]

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

Definition at line 545 of file DensityMaps.hpp.

545 {
546 return getElasticMod(paramElastic,rho); //
547 }
double getElasticMod(double *param, double rho)
Transform densities to Youngs modulus.

◆ 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
xcoordinate of point x
ycoordinate of point y
zcoordinate of point z
dxdimension of the block in x direction
dydimension of the block in y direction
dzdimension of the block in z direction
vec_ixreference vector of x coordinates
vec_iyreference vector of y coordinates
vec_izreference vector of z coordinates
vec_idxreference vector of indices
Returns
error code
Parameters
zcenter of the cube
dzsize of the cube
vec_ixvectors of coord points inside cube
vec_idxindices of inside cube with dimensions dx dy dz
vec_distvector of squared distances between center of the cube and points

Definition at line 334 of file DensityMaps.hpp.

342 {
343 PetscFunctionBegin;
344
345 const double *elemet_size = metaFile.ElementSize();
346 const int *dim_size = metaFile.DimSize();
347
348 int nx = ceil(dx/ (elemet_size[0] * sCale)); // number of points within given size of block - dx
349 int ny = ceil(dy/ (elemet_size[1] * sCale)); // number of points within given size of block - dy
350 int nz = ceil(dz/ (elemet_size[2] * sCale)); // number of points within given size of block - dz
351
352 //const double *offset = metaFile.Offset(); // coords before translation (offset)
353 double x_befor_offset = x;
354 double y_befor_offset = y;
355 double z_befor_offset = z;
356 // double x_befor_offset = x - metaFile.Offset()[0]; // probably its not necessary
357 // double y_befor_offset = y - metaFile.Offset()[1];
358 // double z_befor_offset = z - metaFile.Offset()[2];
359 // center of the cube
360 int ix = ceil(x_befor_offset/ (elemet_size[0] * sCale)); // original coords/indices (before multiplying by elem size)
361 int iy = ceil(y_befor_offset/ (elemet_size[1] * sCale));
362 int iz = ceil(z_befor_offset/ (elemet_size[2] * sCale));
363
364 vec_ix.resize(2*nx); // setting sizes of the vectors equal to number of points
365 vec_iy.resize(2*ny);
366 vec_iz.resize(2*nz);
367
368 int ii = 0;
369 for(int i = 0;i<2*nx;i++) { // vectors of coordinates of points within given cube
370 int idx = ix-nx+i;
371 if(idx >= dim_size[0] || idx < 0) continue; // erasing indices beyond the border
372 vec_ix[ii++] = idx;
373 }
374 vec_ix.resize(ii);
375 ii = 0;
376 for(int i = 0;i<2*ny;i++) {
377 int idx = iy-ny+i;
378 if(idx >= dim_size[1] || idx < 0) continue; // erasing indices beyond the border
379 vec_iy[ii++] = idx;
380 }
381 vec_iy.resize(ii);
382 ii = 0;
383 for(int i = 0;i<2*nz;i++) {
384 int idx = iz-nz+i;
385 if(idx >= dim_size[2] || idx < 0) continue; // erasing indices beyond the border
386 vec_iz[ii++] = idx;
387 }
388 vec_iz.resize(ii);
389
390
391 if (dx <= 0) { //
392 vec_ix.resize(1);
393 vec_ix.push_back(round(x_befor_offset/ (elemet_size[0] * sCale)));
394 }
395 if (dy <= 0) {
396 vec_iy.resize(1);
397 vec_iy.push_back(round(y_befor_offset/ (elemet_size[1] * sCale)));
398 }
399 if (dz <= 0) {
400 vec_iz.resize(1);
401 vec_iz.push_back(round(z_befor_offset/ (elemet_size[2] * sCale)));
402 }
403
404 // vec_idx.resize(8*nx*ny*nz);
405 // vec_dist.resize(8*nx*ny*nz);
406 vec_idx.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
407 vec_dist.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
408 int dm_size0_dim_size1 = dim_size[0] * dim_size[1];
409 // calculating the vector of global indices and vector of distances from center of the cube
410
411 if (fIlter == 0) {
412
413 // vector<double> dy2(vec_iy.size());
414 // {
415 // int ii = 0;
416 // for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
417 // double dyy = (*it_iy*elemet_size[1] - y_befor_offset);
418 // dy2[ii++] = dyy * dyy;
419 // }
420 // }
421 // vector<double> dz2(vec_iz.size());
422 // {
423 // int ii = 0;
424 // for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
425 // double dzz = (*it_iz*elemet_size[2] - z_befor_offset);
426 // dz2[ii++] = dzz * dzz;
427 // }
428 // }
429
430 ii = 0;
431 for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
432 double a = (*it_iz) * dm_size0_dim_size1 ;
433 double dz2 = (*it_iz*(elemet_size[2]*sCale) - z_befor_offset) * (*it_iz*(elemet_size[2]*sCale) - z_befor_offset);
434 for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
435 double b = (*it_iy) * dim_size[0];
436 double dy2 = (*it_iy*(elemet_size[1]*sCale) - y_befor_offset) * (*it_iy*(elemet_size[1]*sCale) - y_befor_offset);
437 for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
438 vec_idx[ii] = a + b + (*it_ix);
439 vec_dist[ii] = ( dz2 + dy2 + (*it_ix*(elemet_size[0]*sCale) - x_befor_offset)
440 * (*it_ix*(elemet_size[0]*sCale) - x_befor_offset) );
441 ii++;
442 // cout << "Ix Iy Iz" << vec_idx[ii] << " " << (*it_ix) << " " << (*it_iy) << " " << (*it_iz) << endl;
443 // cout << "Z " << *it_iz << "-" << vec_idx[ii] / (dim_size[0] * dim_size[1]) << " ";
444 // cout << "Y " << *it_iy << "-" << (vec_idx[ii] - (*it_iz) * (dim_size[0] * dim_size[1])) / dim_size[0] << " ";
445 // cout << "X " << *it_ix << "-" << vec_idx[ii] - (*it_iz) * dim_size[0] * dim_size[1] - (*it_iy) * dim_size[0] << " " << endl;
446 // cout << "POINT "<< ii << " GLOBAL " << vec_idx[ii] <<
447 // " DATA "<< metaFile.ElementData(vec_idx[ii]) << " DISTANCE " << vec_dist[ii] << endl;
448
449
450 // if we require more indices than mhd file has ----> ERROR
451 // if(vec_idx[ii] >= metaFile.Quantity()) {
452 // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong index of voxels (> Quantity)");
453 // }
454
455
456 }
457 }
458 }
459
460 }
461 else {
462
463 ii = 0;
464 for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
465 double a = (*it_iz) * dm_size0_dim_size1 ;
466 for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
467 double b = (*it_iy) * dim_size[0];
468 for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
469 vec_idx[ii++] = a + b + (*it_ix);
470 }
471 }
472 }
473 }
474 // cout << ix << " " << " " << iy << " " << iz <<endl;
475 //cout << "--------------------------------------------------------" << endl;
476 PetscFunctionReturn(0);
477 }
constexpr double a

◆ 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.

526 {
527 PetscFunctionBegin;
528 ierr = mapDensityLinear(paramDensity,vec_data,vec_density); CHKERRQ(ierr);
529 PetscFunctionReturn(0);
530 }
PetscErrorCode mapDensityLinear(double *paramDensity, const vector< double > &vec_data, vector< double > &vec_density)
Transform data to densities.

◆ mapDensityLinear() [2/2]

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

Transform data to densities.

Parameters
aarray of parameters a,b for linear relation of grayscale - density, a * HU + b = density
vec_datadata from metaFile
vec_densityoutput vector with density
Returns
vector of densities

Definition at line 508 of file DensityMaps.hpp.

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

◆ 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
sigmasigma parameter for gausssian smoothing
vec_densityvector with non-smoothed densities
vec_distvector with squared distances form center of the cubes

Definition at line 581 of file DensityMaps.hpp.

586 {
587 //vec_density_smooth.resize(vec_density.size());
588 kErnel.resize(vec_density.size());
589 //calculating kernel
590 int i = 0;
591 double sum = 0;
592 const double k = 2 * sigma*sigma;
593 const double m = (1 / sqrt( M_PI * k));
594 for(int ii = 0; ii < vec_dist.size(); ii++) {
595 if (vec_density[ii] < tHreshold) continue;
596 kErnel[i] = m * exp(- (vec_dist[ii]) / k); // distance is already squared
597 sum += kErnel[i++];
598 //cout <<"KERNEL"<< kernel[i] << endl;;
599 }
600 kErnel.resize(i);
601 ///< normalisation of kernel and smooting values
602 int ii = 0;
603 double smooth_data = 0;
604 for(vector<double>::const_iterator vec_itr = vec_density.begin();vec_itr!=vec_density.end();vec_itr++) {
605 if (*vec_itr < tHreshold) continue;
606 kErnel[ii] /= sum;
607 //cout << kernel[ii] << endl;
608 smooth_data += (*vec_itr) * kErnel[ii++];
609 }
610 if (sum == 0) return 0; ///< avoid division by zero (ad exception) this should not happen /FIXME
611 return smooth_data;
612 }
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'k', 3 > k

◆ mapGaussSmooth() [2/2]

double BoneRemodeling::DataFromMetaIO::mapGaussSmooth ( const vector< double > &  vec_density,
const vector< double > &  vec_dist 
)
inline
Parameters
vec_densityvector with non-smoothed densities
vec_distvector with distances form center of the cubes

Definition at line 614 of file DensityMaps.hpp.

617 {
618 return mapGaussSmooth(sIgma,vec_density,vec_dist);
619 }
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)

◆ medianFilter()

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

function returns median of data from given vector

Parameters
vec_densityvector 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.

553 {
554 vector<double> vec_copy; ///< copy of the vector to allow sorting data
555 vec_copy.resize(vec_density.size());
556 int i = 0;
557 for(vector<double>::const_iterator vec_it = vec_density.begin();vec_it!=vec_density.end();vec_it++){
558 if (*vec_it < tHreshold) continue; ///< applying threshold
559 vec_copy[i++] = *vec_it;
560 //cout << *vec_it << "-";
561 }
562 vec_copy.resize(i);
563
564 //cout << " " << endl;
565 sort(vec_copy.begin(), vec_copy.end()); ///< sorting vector
566 //cout << " " << endl;
567 if (vec_copy.size() % 2 == 0) {
568 return (vec_copy[vec_copy.size() / 2 - 1] + vec_copy[vec_copy.size() / 2]) / 2.0;
569 } else {
570 return vec_copy[vec_copy.size() / 2];
571 }
572 return 0;
573 }

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: