v0.14.0
mhd_merger.cpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2 * \ingroup bone_remodelling
3 
4 * MoFEM is free software: you can redistribute it and/or modify it under
5 * the terms of the GNU Lesser General Public License as published by the
6 * Free Software Foundation, either version 3 of the License, or (at your
7 * option) any later version.
8 *
9 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12 * License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
16 
17 #include <MoFEM.hpp>
18 using namespace MoFEM;
19 
20 //include <DensityMaps.hpp>
21 //using namespace DensityMaps;
22 
23 #ifdef WITH_METAIO
24 #include <metaImage.h>
25 #endif
26 
27 static char help[] = "\n";
28 #ifdef WITH_METAIO
29 
30 float mapGaussSmooth(
31  const float sigma, ///< sigma parameter for gausssian smoothing
32  const vector<float> &vec_density, ///< vector with non-smoothed densities
33  const vector<float> &vec_dist ///< vector with squared distances form center of the cubes
34  // vector<float> &vec_density_smooth ///< vector with the smoothed densities
35 );
36 
37 
38 PetscErrorCode getInidcesForGivenCoordinate(
39  const float x,const float y,const float z, ///< center of the cube
40  const float size, ///< size of the cube
41  vector<int> &vec_ix, ///< vectors of coord points inside cube
42  vector<int> &vec_iy, ///<
43  vector<int> &vec_iz, ///<
44  vector<int> &vec_idx, ///< indices of inside cube with dimensions dx dy dz
45  vector<float> &vec_dist,
46  vector<float> &vec_density, ///< vector of squared distances between center of the cube and points
47  MetaImage &meta_file
48 
49 );
50 
51 int main(int argc, char *argv[]) {
52  PetscInitialize(&argc,&argv,(char *)0,help);
53 
54  PetscBool flg = PETSC_TRUE;
55  char meta_file_name_1[255];
56  PetscBool meta_flg_file_1 = PETSC_FALSE;
57  char meta_file_name_2[255];
58  PetscBool meta_flg_file_2 = PETSC_FALSE;
59 
60  char output_name[255] = "newImage.mhd";
61  PetscBool gauss_flag = PETSC_FALSE;
62 
63 
64  // Parameters from command line
65  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Meta merge","none"); CHKERRQ(ierr);
66 
67  // file name of first mhd file
68  ierr = PetscOptionsString(
69  "-meta_image1",
70  "meta image name","",
71  "file.mhd",meta_file_name_1,255,&meta_flg_file_1
72  ); CHKERRQ(ierr);
73 
74  ierr = PetscOptionsString(
75  "-o",
76  "output file","",
77  "newImage.mhd",output_name,255,PETSC_NULL
78  ); CHKERRQ(ierr);
79 
80 
81  //file name of second mhd file
82  ierr = PetscOptionsString(
83  "-meta_image2",
84  "meta image name","",
85  "file.mhd",meta_file_name_2,255,&meta_flg_file_2
86  ); CHKERRQ(ierr);
87  //translation vector in case when images are shifted to each other
88  ierr = PetscOptionsEnd(); CHKERRQ(ierr);
89  int my_vector[3] = {0, 0, 0};
90  int my_vector_gauss[2] = {3, 2};
91  int nmax = 3;
92  ierr = PetscOptionsGetIntArray(PETSC_NULL,"-my_vector",my_vector,&nmax,&flg); CHKERRQ(ierr);
93  nmax=2;
94  ierr = PetscOptionsGetIntArray(PETSC_NULL,"-gauss",my_vector_gauss,&nmax,&gauss_flag); CHKERRQ(ierr);
95 
96  //Load data for both files
97  MetaImage meta_file1(meta_file_name_1);
98  cout << "Quantity of file 1: " << meta_file1.Quantity() << endl;
99  const int *dim_size1 = meta_file1.DimSize();
100  cout << "Image dimension: " << dim_size1[0] << " " << dim_size1[1] << " " << dim_size1[2] << endl;
101  const double *element_size1 = meta_file1.ElementSize();
102  cout << "Image element size: " << element_size1[0] << " " << element_size1[1] << " " << element_size1[2] << endl;
103 
104  MetaImage meta_file2(meta_file_name_2);
105  cout << "\nQuantity of file 2: " << meta_file2.Quantity() << endl;
106  const int *dim_size2 = meta_file2.DimSize();
107  cout << "Image dimension: " << dim_size2[0] << " " << dim_size2[1] << " " << dim_size2[2] << endl;
108  const double *element_size2 = meta_file2.ElementSize();
109  cout << "Image element size: " << element_size2[0] << " " << element_size2[1] << " " << element_size2[2] << endl;
110  //show translation vector
111  cout<< "\nTranslation vector: " << my_vector[0] << " " << my_vector[1] << " " << my_vector[2] << endl;
112 
113 
114  //check if images have the same resolution/size (not needed)
115  // if(dim_size2[0] != dim_size1[0] || dim_size2[1] != dim_size2[1]) {
116  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"images have different size! Rescale them.");
117  // }
118  // if(element_size1[0] != element_size2[0] || element_size1[1] != element_size2[1] || element_size1[2] != element_size2[2]) {
119  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"images have different resolution! Rescale them.");
120  // }
121 
122  int quantity = meta_file1.Quantity();//number of voxels in new file
123  if(meta_flg_file_2) quantity += meta_file2.Quantity();
124  std::vector<float> data_vec(quantity); //allocating data for entire new file
125  float *data = &data_vec[0];
126 
127  //for gaussian smoothing
128 
129 
130 
131 
132  //loop over all elements in file 1
133  float sum1 = 0;
134  float sum2= 0;
135  int i = 0;
136  for(int i0 = 0; i0 != dim_size1[2]; i0++) {
137  for(int i1 = 0; i1 != dim_size1[1]; i1++) {
138  for(int i2 = 0; i2 != dim_size1[0]; i2++) {
139  // 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;
140  vector<float> vec_density, vec_dist;
141  vector<int> vec_ix,vec_iy,vec_iz, vec_idx;
142 
143  data[i] = meta_file1.ElementData(i);
144  sum1 += data[i];
145 
147  i2,i1,i0,
148  my_vector_gauss[0],
149  vec_ix,
150  vec_iy,
151  vec_iz,
152  vec_idx,
153  vec_dist,
154  vec_density,
155  meta_file1
156  );
157 
158  if (gauss_flag) {
159  data[i] = mapGaussSmooth(my_vector_gauss[1],vec_density,vec_dist);
160  }
161  sum2 += data[i];
162  // cout << "i " << i << " " << i0 * dim_size[0] * dim_size[1] + i1 * dim_size[0] + i2 << endl;
163  // cout << "Z " << i0 << " " << i / (dim_size[0] * dim_size[1]) << endl;
164  // cout << "Y " << i1 << " " << (i - i0 * (dim_size[0] * dim_size[1])) / dim_size[0] << endl;
165  // cout << "X " << i2 << " " << i - i0 * dim_size[0] * dim_size[1] - i1 * dim_size[2] << endl;
166  i++;
167  }
168  }
169  }
170 
171  if(meta_flg_file_2) {
172  //loop for second file
173  for(int i0 = 0; i0 != dim_size2[2]; i0++) {
174  for(int i1 = 0; i1 != dim_size2[1]; i1++) {
175  for(int i2 = 0; i2 != dim_size2[0]; i2++) {
176  //throw out voxels beyond the border after shifting by vector from command line
177  int z = dim_size2[0]*dim_size2[1] * ((i0 - my_vector[2])%dim_size2[2]);
178  int y = ((i1 - my_vector[1])%dim_size2[1])*dim_size2[1];
179  int x = (i2 - my_vector[0])%dim_size2[0];
180  if (i0 - my_vector[2] >= dim_size2[2] || i1 - my_vector[1] >= dim_size2[1] || i2 - my_vector[0] >= dim_size2[0] ||
181  i0 - my_vector[2] < 0 || i1 - my_vector[1] < 0 || i2 - my_vector[0] < 0){
182  data[i++] = -1024; continue;
183  }
184  data[i++] = meta_file2.ElementData(x + y + z);
185  }
186  }
187  }
188  }
189  int Z_dimensions = dim_size1[2]; // Z_dimension of the new file
190 
191  if(meta_flg_file_2) Z_dimensions += dim_size2[2];
192  int global_dimsize[3] = {dim_size1[0], dim_size1[1], Z_dimensions};
193  cout << "MHD converter. Writing file." << endl;
194  /** We create a simple 3D metaImage */
195  MetaImage newImage(3, global_dimsize, element_size1, MET_CHAR, 1, data); // Create a new 3D metaImage
196  newImage.FileName(output_name); // Define the name of the file
197  newImage.ObjectTypeName("Image"); // Define the type of the object
198  /** Write the image in binary format */
199  /** Write the object */
200  newImage.BinaryData(true);
201  newImage.ElementType(MET_FLOAT);
202  newImage.Write(output_name);
203  if (gauss_flag) printf("Sum without smoothing: %g Sum with smoothing: %g \n",sum1,sum2 );
204  cout << "File saved." << endl;
205  /** Clear the object */
206  newImage.Clear();
207 
208 
209  PetscFinalize();
210  PetscFunctionReturn(0);
211 }
212 
214  const float sigma, ///< sigma parameter for gausssian smoothing
215  const vector<float> &vec_density, ///< vector with non-smoothed densities
216  const vector<float> &vec_dist ///< vector with squared distances form center of the cubes
217  // vector<float> &vec_density_smooth ///< vector with the smoothed densities
218 ) {
219  //vec_density_smooth.resize(vec_density.size());
220  vector<float> kErnel;
221  kErnel.resize(vec_density.size());
222  //calculating kernel
223  int i = 0;
224  float sum = 0;
225  const float k = 2 * sigma*sigma;
226  const float m = (1 / sqrt( M_PI * k));
227  for(int ii = 0; ii < vec_dist.size(); ii++) {
228  kErnel[i] = m * exp(- (vec_dist[ii]) / k); // distance is already squared
229  sum += kErnel[i++];
230  //cout <<"KERNEL"<< kernel[i] << endl;;
231  }
232  kErnel.resize(i);
233  ///< normalisation of kernel and smooting values
234  int ii = 0;
235  float smooth_data = 0;
236  for(vector<float>::const_iterator vec_itr = vec_density.begin();vec_itr!=vec_density.end();vec_itr++) {
237  kErnel[ii] /= sum;
238  smooth_data += (*vec_itr) * kErnel[ii++];
239  }
240  if (sum == 0) return 0; ///< avoid division by zero (ad exception) this should not happen /FIXME
241  return smooth_data;
242 }
243 
244 
246  const float x,const float y,const float z, ///< center of the cube
247  const float size, ///< size of the cube
248  vector<int> &vec_ix, ///< vectors of coord points inside cube
249  vector<int> &vec_iy, ///<
250  vector<int> &vec_iz, ///<
251  vector<int> &vec_idx, ///< indices of inside cube with dimensions dx dy dz
252  vector<float> &vec_dist,
253  vector<float> &vec_density, ///< vector of squared distances between center of the cube and points
254  MetaImage &meta_file
255 
256 ) {
257  PetscFunctionBegin;
258 
259  const float elemet_size[3] = {1,1,1};
260  const int *dim_size = meta_file.DimSize();
261 
262  int nx = ceil(size/ (elemet_size[0])); // number of points within given size of block - dx
263  int ny = ceil(size/ (elemet_size[1])); // number of points within given size of block - dy
264  int nz = ceil(size/ (elemet_size[2])); // number of points within given size of block - dz
265 
266  //const float *offset = meta_file.Offset(); // coords before translation (offset)
267  float x_befor_offset = x;
268  float y_befor_offset = y;
269  float z_befor_offset = z;
270  // float x_befor_offset = x - meta_file.Offset()[0]; // probably its not necessary
271  // float y_befor_offset = y - meta_file.Offset()[1];
272  // float z_befor_offset = z - meta_file.Offset()[2];
273  // center of the cube
274  int ix = ceil(x_befor_offset/ (elemet_size[0])); // original coords/indices (before multiplying by elem size)
275  int iy = ceil(y_befor_offset/ (elemet_size[1]));
276  int iz = ceil(z_befor_offset/ (elemet_size[2]));
277 
278  vec_ix.resize(2*nx); // setting sizes of the vectors equal to number of points
279  vec_iy.resize(2*ny);
280  vec_iz.resize(2*nz);
281 
282  int ii = 0;
283  for(int i = 0;i<2*nx;i++) { // vectors of coordinates of points within given cube
284  int idx = ix-nx+i;
285  if(idx >= dim_size[0] || idx < 0) continue; // erasing indices beyond the border
286  vec_ix[ii++] = idx;
287  }
288  vec_ix.resize(ii);
289  ii = 0;
290  for(int i = 0;i<2*ny;i++) {
291  int idx = iy-ny+i;
292  if(idx >= dim_size[1] || idx < 0) continue; // erasing indices beyond the border
293  vec_iy[ii++] = idx;
294  }
295  vec_iy.resize(ii);
296  ii = 0;
297  for(int i = 0;i<2*nz;i++) {
298  int idx = iz-nz+i;
299  if(idx >= dim_size[2] || idx < 0) continue; // erasing indices beyond the border
300  vec_iz[ii++] = idx;
301  }
302  vec_iz.resize(ii);
303 
304 
305  // vec_idx.resize(8*nx*ny*nz);
306  // vec_dist.resize(8*nx*ny*nz);
307  vec_idx.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
308  vec_dist.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
309  vec_density.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
310  int dm_size0_dim_size1 = dim_size[0] * dim_size[1];
311  // calculating the vector of global indices and vector of distances from center of the cube
312 
313 
314 
315  // vector<float> dy2(vec_iy.size());
316  // {
317  // int ii = 0;
318  // for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
319  // float dyy = (*it_iy*elemet_size[1] - y_befor_offset);
320  // dy2[ii++] = dyy * dyy;
321  // }
322  // }
323  // vector<float> dz2(vec_iz.size());
324  // {
325  // int ii = 0;
326  // for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
327  // float dzz = (*it_iz*elemet_size[2] - z_befor_offset);
328  // dz2[ii++] = dzz * dzz;
329  // }
330  // }
331 
332  ii = 0;
333  for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
334  float a = (*it_iz) * dm_size0_dim_size1 ;
335  float dz2 = (*it_iz*(elemet_size[2]) - z_befor_offset) * (*it_iz*(elemet_size[2]) - z_befor_offset);
336  for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
337  float b = (*it_iy) * dim_size[0];
338  float dy2 = (*it_iy*(elemet_size[1]) - y_befor_offset) * (*it_iy*(elemet_size[1]) - y_befor_offset);
339  for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
340  vec_idx[ii] = a + b + (*it_ix);
341  vec_density[ii] = meta_file.ElementData(vec_idx[ii]);
342  vec_dist[ii] = ( dz2 + dy2 + (*it_ix*(elemet_size[0]) - x_befor_offset)
343  * (*it_ix*(elemet_size[0]) - x_befor_offset) );
344  ii++;
345  // cout << "Ix Iy Iz" << vec_idx[ii] << " " << (*it_ix) << " " << (*it_iy) << " " << (*it_iz) << endl;
346  // cout << "Z " << *it_iz << "-" << vec_idx[ii] / (dim_size[0] * dim_size[1]) << " ";
347  // cout << "Y " << *it_iy << "-" << (vec_idx[ii] - (*it_iz) * (dim_size[0] * dim_size[1])) / dim_size[0] << " ";
348  // cout << "X " << *it_ix << "-" << vec_idx[ii] - (*it_iz) * dim_size[0] * dim_size[1] - (*it_iy) * dim_size[0] << " " << endl;
349  // cout << "POINT "<< ii << " GLOBAL " << vec_idx[ii] <<
350  // " DATA "<< metaFile.ElementData(vec_idx[ii]) << " DISTANCE " << vec_dist[ii] << endl;
351 
352 
353  // if we require more indices than mhd file has ----> ERROR
354  // if(vec_idx[ii] >= metaFile.Quantity()) {
355  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"wrong index of voxels (> Quantity)");
356  // }
357 
358 
359  }
360  }
361  }
362 
363 
364 
365  // cout << ix << " " << " " << iy << " " << iz <<endl;
366  //cout << "--------------------------------------------------------" << endl;
367  PetscFunctionReturn(0);
368 }
369 
370 
371 #endif // WITH_METAIO
main
int main(int argc, char *argv[])
Definition: mhd_merger.cpp:51
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)
Definition: mhd_merger.cpp:245
MoFEM.hpp
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::PetscOptionsGetIntArray
PetscErrorCode PetscOptionsGetIntArray(PetscOptions options, const char pre[], const char name[], PetscInt dvalue[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:215
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
mapGaussSmooth
float mapGaussSmooth(const float sigma, const vector< float > &vec_density, const vector< float > &vec_dist)
Definition: mhd_merger.cpp:213
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
help
static char help[]
Definition: mhd_merger.cpp:27
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20