v0.14.0
Loading...
Searching...
No Matches
Functions | Variables
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
zcenter of the cube
sizesize of the cube
vec_ixvectors of coord points inside cube
vec_idxindices of inside cube with dimensions dx dy dz
vec_densityvector of squared distances between center of the cube and points

Definition at line 245 of file mhd_merger.cpp.

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}
constexpr double a
FTensor::Index< 'i', SPACE_DIM > i

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

51 {
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}
static char help[]
Definition: mhd_merger.cpp:27
float mapGaussSmooth(const float sigma, const vector< float > &vec_density, const vector< float > &vec_dist)
Definition: mhd_merger.cpp:213
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode PetscOptionsGetIntArray(PetscOptions options, const char pre[], const char name[], PetscInt dvalue[], PetscInt *nmax, PetscBool *set)

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

Definition at line 213 of file mhd_merger.cpp.

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}
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'k', 3 > k

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 27 of file mhd_merger.cpp.