mhd_merger.cpp

Go to the documentation of this file.

/* This file is part of MoFEM.
* \ingroup bone_remodelling
 
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
//include <DensityMaps.hpp>
//using namespace DensityMaps;
 
#ifdef WITH_METAIO
#include <metaImage.h>
#endif
 
static char help[] = "\n";
#ifdef WITH_METAIO
 
float mapGaussSmooth(
  const float sigma,                                         ///< sigma parameter for gausssian smoothing
  const vector<float> &vec_density,                          ///< vector with non-smoothed densities
  const vector<float> &vec_dist                              ///< vector with squared distances form center of the cubes
  //  vector<float> &vec_density_smooth                        ///< vector with the smoothed densities
);
 
 
PetscErrorCode getInidcesForGivenCoordinate(
  const float x,const float y,const float z,     ///< center of the cube
  const float size,  ///< size of the cube
  vector<int> &vec_ix,                    ///<   vectors of coord points inside cube
  vector<int> &vec_iy,                    ///<
  vector<int> &vec_iz,                    ///<
  vector<int> &vec_idx,                   ///<   indices of inside cube with dimensions dx dy dz
  vector<float> &vec_dist,
  vector<float> &vec_density,                ///<   vector of squared distances between center of the cube and points
  MetaImage &meta_file
 
);
 
int main(int argc, char *argv[]) {
    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
    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
    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);
}
 
float mapGaussSmooth(
  const float sigma,                                         ///< sigma parameter for gausssian smoothing
  const vector<float> &vec_density,                          ///< vector with non-smoothed densities
  const vector<float> &vec_dist                              ///< vector with squared distances form center of the cubes
  //  vector<float> &vec_density_smooth                        ///< vector with the smoothed densities
) {
  //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;
}
 
 
PetscErrorCode getInidcesForGivenCoordinate(
  const float x,const float y,const float z,     ///< center of the cube
  const float size,  ///< size of the cube
  vector<int> &vec_ix,                    ///<   vectors of coord points inside cube
  vector<int> &vec_iy,                    ///<
  vector<int> &vec_iz,                    ///<
  vector<int> &vec_idx,                   ///<   indices of inside cube with dimensions dx dy dz
  vector<float> &vec_dist,
  vector<float> &vec_density,                ///<   vector of squared distances between center of the cube and points
  MetaImage &meta_file
 
) {
  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);
}
 
 
#endif // WITH_METAIO