mhd_merger.cpp File Reference

#include <MoFEM.hpp>
#include <metaImage.h>

Go to the source code of this file.

Functions
float	mapGaussSmooth (const float sigma, const vector< float > &vec_density, const vector< float > &vec_dist)
PetscErrorCode	getInidcesForGivenCoordinate (const float x, const float y, const float z, const float size, vector< int > &vec_ix, vector< int > &vec_iy, vector< int > &vec_iz, vector< int > &vec_idx, vector< float > &vec_dist, vector< float > &vec_density, MetaImage &meta_file)
int	main (int argc, char *argv[])

Variables
static char	help [] = "\n"

Function Documentation

getInidcesForGivenCoordinate()

PetscErrorCode getInidcesForGivenCoordinate	(	const float	x,
		const float	y,
		const float	z,
		const float	size,
		vector< int > &	vec_ix,
		vector< int > &	vec_iy,
		vector< int > &	vec_iz,
		vector< int > &	vec_idx,
		vector< float > &	vec_dist,
		vector< float > &	vec_density,
		MetaImage &	meta_file
	)

Parameters

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

Definition at line 245 of file mhd_merger.cpp.

  {
  PetscFunctionBegin;
 
  const float elemet_size[3] = {1,1,1};
  const int *dim_size = meta_file.DimSize();
 
  int nx = ceil(size/ (elemet_size[0]));                   // number of points within given size of block - dx
  int ny = ceil(size/ (elemet_size[1]));                   // number of points within given size of block - dy
  int nz = ceil(size/ (elemet_size[2]));                   // number of points within given size of block - dz
 
  //const float *offset = meta_file.Offset();           // coords before translation (offset)
  float x_befor_offset = x;
  float y_befor_offset = y;
  float z_befor_offset = z;
  // float x_befor_offset = x - meta_file.Offset()[0];  // probably its not necessary
  // float y_befor_offset = y - meta_file.Offset()[1];
  // float z_befor_offset = z - meta_file.Offset()[2];
  // center of the cube
  int ix  = ceil(x_befor_offset/ (elemet_size[0]));     // original coords/indices (before multiplying by elem size)
  int iy  = ceil(y_befor_offset/ (elemet_size[1]));
  int iz  = ceil(z_befor_offset/ (elemet_size[2]));
 
  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);
 
 
  // 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());
  vec_density.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
 
 
 
    // vector<float> dy2(vec_iy.size());
    // {
    //   int ii = 0;
    //   for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
    //     float dyy = (*it_iy*elemet_size[1] - y_befor_offset);
    //     dy2[ii++] =  dyy * dyy;
    //   }
    // }
    // vector<float> dz2(vec_iz.size());
    // {
    //   int ii = 0;
    //   for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
    //     float 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++) {
      float a = (*it_iz) * dm_size0_dim_size1 ;
      float dz2 = (*it_iz*(elemet_size[2]) - z_befor_offset) * (*it_iz*(elemet_size[2]) - z_befor_offset);
      for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
        float b = (*it_iy) * dim_size[0];
        float dy2 = (*it_iy*(elemet_size[1]) - y_befor_offset) * (*it_iy*(elemet_size[1]) - 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_density[ii] = meta_file.ElementData(vec_idx[ii]);
          vec_dist[ii] =  ( dz2 + dy2 + (*it_ix*(elemet_size[0]) - x_befor_offset)
           * (*it_ix*(elemet_size[0]) - 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)");
          // }
 
 
        }
      }
    }
 
 
 
  // cout << ix << " " << " " << iy << " " << iz <<endl;
  //cout << "--------------------------------------------------------" << endl;
  PetscFunctionReturn(0);
}

main()

int main	(	int	argc,
		char *	argv[]
	)

We create a simple 3D metaImage

Write the image in binary format

Write the object

Clear the object

Definition at line 51 of file mhd_merger.cpp.

                                 {
    PetscInitialize(&argc,&argv,(char *)0,help);
 
    PetscBool flg = PETSC_TRUE;
    char meta_file_name_1[255];
    PetscBool meta_flg_file_1 = PETSC_FALSE;
    char meta_file_name_2[255];
    PetscBool meta_flg_file_2 = PETSC_FALSE;
 
    char output_name[255] = "newImage.mhd";
    PetscBool gauss_flag = PETSC_FALSE;
 
 
    // Parameters from command line
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Meta merge","none"); CHKERRQ(ierr);
 
    // file name of first mhd file
    ierr = PetscOptionsString(
      "-meta_image1",
      "meta image name","",
      "file.mhd",meta_file_name_1,255,&meta_flg_file_1
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsString(
      "-o",
      "output file","",
      "newImage.mhd",output_name,255,PETSC_NULL
    ); CHKERRQ(ierr);
 
 
    //file name of second mhd file
    ierr = PetscOptionsString(
      "-meta_image2",
      "meta image name","",
      "file.mhd",meta_file_name_2,255,&meta_flg_file_2
    ); CHKERRQ(ierr);
    //translation vector in case when images are shifted to each other
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);
    int my_vector[3] = {0, 0, 0};
    int my_vector_gauss[2] = {3, 2};
    int nmax = 3;
    ierr = PetscOptionsGetIntArray(PETSC_NULL,"-my_vector",my_vector,&nmax,&flg); CHKERRQ(ierr);
    nmax=2;
    ierr = PetscOptionsGetIntArray(PETSC_NULL,"-gauss",my_vector_gauss,&nmax,&gauss_flag); CHKERRQ(ierr);
 
    //Load data for both files
    MetaImage meta_file1(meta_file_name_1);
    cout << "Quantity of file 1: " << meta_file1.Quantity() << endl;
    const int *dim_size1 = meta_file1.DimSize();
    cout << "Image dimension: " << dim_size1[0] << " " << dim_size1[1] << " " << dim_size1[2] << endl;
    const double *element_size1 = meta_file1.ElementSize();
    cout << "Image element size: " << element_size1[0] << " " << element_size1[1] << " " << element_size1[2] << endl;
 
    MetaImage meta_file2(meta_file_name_2);
    cout << "\nQuantity of file 2: " << meta_file2.Quantity() << endl;
    const int *dim_size2 = meta_file2.DimSize();
    cout << "Image dimension: " << dim_size2[0] << " " << dim_size2[1] << " " << dim_size2[2] << endl;
    const double *element_size2 = meta_file2.ElementSize();
    cout << "Image element size: " << element_size2[0] << " " << element_size2[1] << " " << element_size2[2] << endl;
    //show translation vector
    cout<< "\nTranslation vector: " << my_vector[0] << " " << my_vector[1] << " " << my_vector[2] << endl;
 
 
    //check if images have the same resolution/size (not needed)
    // if(dim_size2[0] != dim_size1[0] || dim_size2[1] != dim_size2[1]) {
    //   SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"images have different size! Rescale them.");
    // }
    // if(element_size1[0] != element_size2[0] || element_size1[1] != element_size2[1] || element_size1[2] != element_size2[2]) {
    //   SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"images have different resolution! Rescale them.");
    // }
 
    int quantity = meta_file1.Quantity();//number of voxels in new file
    if(meta_flg_file_2) quantity += meta_file2.Quantity();
    std::vector<float> data_vec(quantity); //allocating data for entire new file
    float *data = &data_vec[0];
 
    //for gaussian smoothing
 
 
 
 
    //loop over all elements in file 1
    float sum1 = 0;
    float sum2= 0;
    int i = 0;
    for(int i0 = 0; i0 != dim_size1[2]; i0++) {
      for(int i1 = 0; i1 != dim_size1[1]; i1++) {
        for(int i2 = 0; i2 != dim_size1[0]; i2++) {
          //    int i = i0*dim_size[0]*dim_size[1]+i1*dim_size[0]+i2; // int i = i2*dim_size[0]*dim_size[1]+i1*dim_size[0]+i0;
          vector<float> vec_density, vec_dist;
          vector<int> vec_ix,vec_iy,vec_iz, vec_idx;
 
          data[i] = meta_file1.ElementData(i);
          sum1 += data[i];
 
          getInidcesForGivenCoordinate(
            i2,i1,i0,
            my_vector_gauss[0],
            vec_ix,
            vec_iy,
            vec_iz,
            vec_idx,
            vec_dist,
            vec_density,
            meta_file1
          );
 
          if (gauss_flag) {
             data[i] = mapGaussSmooth(my_vector_gauss[1],vec_density,vec_dist);
          }
          sum2 += data[i];
          // cout << "i " << i << " " <<  i0 * dim_size[0] *  dim_size[1] + i1 * dim_size[0] + i2 << endl;
          // cout << "Z " << i0 << " " <<  i / (dim_size[0] * dim_size[1]) << endl;
          // cout << "Y " << i1 << " " << (i - i0 * (dim_size[0] * dim_size[1])) / dim_size[0] << endl;
          // cout << "X " << i2 << " " << i - i0 * dim_size[0] * dim_size[1] - i1 * dim_size[2] << endl;
          i++;
        }
      }
    }
 
    if(meta_flg_file_2) {
    //loop for second file
    for(int i0 = 0; i0 != dim_size2[2]; i0++) {
      for(int i1 = 0; i1 != dim_size2[1]; i1++) {
        for(int i2 = 0; i2 != dim_size2[0]; i2++) {
          //throw out voxels beyond the border after shifting by vector from command line
          int z = dim_size2[0]*dim_size2[1] * ((i0 - my_vector[2])%dim_size2[2]);
          int y = ((i1 - my_vector[1])%dim_size2[1])*dim_size2[1];
          int x = (i2 - my_vector[0])%dim_size2[0];
          if (i0 - my_vector[2] >= dim_size2[2] || i1 - my_vector[1] >= dim_size2[1] || i2 - my_vector[0] >= dim_size2[0] ||
          i0 - my_vector[2] < 0 || i1 - my_vector[1] < 0 || i2 - my_vector[0] < 0){
            data[i++] = -1024; continue;
          }
          data[i++] = meta_file2.ElementData(x + y + z);
        }
      }
    }
  }
    int Z_dimensions = dim_size1[2]; // Z_dimension of the new file
 
    if(meta_flg_file_2) Z_dimensions += dim_size2[2];
    int global_dimsize[3] = {dim_size1[0], dim_size1[1], Z_dimensions};
    cout << "MHD converter. Writing file." << endl;
    /** We create a simple 3D metaImage */
    MetaImage newImage(3, global_dimsize, element_size1, MET_CHAR, 1, data); // Create a new 3D metaImage
    newImage.FileName(output_name); // Define the name of the file
    newImage.ObjectTypeName("Image"); // Define the type of the object
    /** Write the image in binary format */
    /** Write the object */
    newImage.BinaryData(true);
    newImage.ElementType(MET_FLOAT);
    newImage.Write(output_name);
    if (gauss_flag) printf("Sum without smoothing: %g Sum with smoothing: %g \n",sum1,sum2 );
    cout << "File saved." << endl;
    /** Clear the object */
    newImage.Clear();
 
 
    PetscFinalize();
    PetscFunctionReturn(0);
}

mapGaussSmooth()

float mapGaussSmooth	(	const float	sigma,
		const vector< float > &	vec_density,
		const vector< float > &	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 213 of file mhd_merger.cpp.

  {
  //vec_density_smooth.resize(vec_density.size());
  vector<float> kErnel;
  kErnel.resize(vec_density.size());
  //calculating kernel
  int i = 0;
  float sum = 0;
  const float k = 2 * sigma*sigma;
  const float m = (1 / sqrt( M_PI * k));
  for(int ii = 0; ii < vec_dist.size(); ii++) {
    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;
  float smooth_data = 0;
  for(vector<float>::const_iterator vec_itr = vec_density.begin();vec_itr!=vec_density.end();vec_itr++) {
    kErnel[ii] /= sum;
    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;
}

Variable Documentation

help

char help[] = "\n"

static

Definition at line 27 of file mhd_merger.cpp.