SPHinXsys  alpha version
image_mhd.hpp
Go to the documentation of this file.
1 /* -------------------------------------------------------------------------*
2 * SPHinXsys *
3 * --------------------------------------------------------------------------*
4 * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
5 * Hydrodynamics for industrial compleX systems. It provides C++ APIs for *
6 * physical accurate simulation and aims to model coupled industrial dynamic *
7 * systems including fluid, solid, multi-body dynamics and beyond with SPH *
8 * (smoothed particle hydrodynamics), a meshless computational method using *
9 * particle discretization. *
10 * *
11 * SPHinXsys is partially funded by German Research Foundation *
12 * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1 *
13 * and HU1527/12-1. *
14 * *
15 * Portions copyright (c) 2017-2020 Technical University of Munich and *
16 * the authors' affiliations. *
17 * *
18 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
19 * not use this file except in compliance with the License. You may obtain a *
20 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
21 * *
22 * --------------------------------------------------------------------------*/
30 #ifndef IMAGE_MHD_HPP
31 #define IMAGE_MHD_HPP
32 
33 #ifndef __EMSCRIPTEN__
34 
35 #include "image_mhd.h"
36 #include "boost/algorithm/string.hpp"
37 
38 namespace SPH {
39 
40  template<typename T, int nDims>
41  ImageMHD<T, nDims>::ImageMHD(std::string full_path_to_file):
42  objectType_("Image"),
43  ndims_(nDims),
44  binaryData_(true),
45  binaryDataByteOrderMSB_(false),
46  compressedData_(false),
47  transformMatrix_(Mat3d(1.0)),
48  offset_(Vec3d(0.0,0.0,0.0)),
49  centerOfRotation_(Vec3d(0.0,0.0,0.0)),
50  elementSpacing_(Vec3d(1.0,1.0,1.0)),
51  dimSize_(Vec3i(1,1,1)),
52  width_(dimSize_[0]),
53  height_(dimSize_[1]),
54  depth_(dimSize_[2]),
55  size_(dimSize_[0] * dimSize_[1] * dimSize_[2]),
56  anatomicalOrientation_("???"),
57  elementType_(MET_FLOAT),
58  elementDataFile_(""),
59  min_value_(Infinity),
60  max_value_(-Infinity),
61  data_(nullptr)
62  {
63  //- read mhd file
64  std::ifstream dataFile(full_path_to_file, std::ifstream::in);
65  std::string file_path_to_raw_file;
66  if (dataFile.is_open())
67  {
68  std::string line;
69  std::vector<std::string> values;
70  while (std::getline(dataFile, line))
71  {
72  std::stringstream ss(line);
73  std::vector<std::string> elements;
74  split(line, '=', elements);
75  if (elements.size() == 2)
76  {
77  boost::trim(elements[0]);
78  boost::trim(elements[1]);
79  if (elements[0].compare("TransformMatrix") == 0)
80  {
81  std::vector<std::string> values;
82  split(elements[1], ' ', values);
83  transformMatrix_[0][0] = std::stof(values[0]);
84  transformMatrix_[0][1] = std::stof(values[1]);
85  transformMatrix_[0][2] = std::stof(values[2]);
86  transformMatrix_[1][0] = std::stof(values[3]);
87  transformMatrix_[1][1] = std::stof(values[4]);
88  transformMatrix_[1][2] = std::stof(values[5]);
89  transformMatrix_[2][0] = std::stof(values[6]);
90  transformMatrix_[2][1] = std::stof(values[7]);
91  transformMatrix_[2][2] = std::stof(values[8]);
92  }
93  else if (elements[0].compare("Offset") == 0)
94  {
95  std::vector<std::string> values;
96  split(elements[1], ' ', values);
97  offset_[0] = std::stof(values[0]);
98  offset_[1] = std::stof(values[1]);
99  offset_[2] = std::stof(values[2]);
100  }
101  else if (elements[0].compare("ElementSpacing") == 0)
102  {
103  std::vector<std::string> values;
104  split(elements[1], ' ', values);
105  elementSpacing_[0] = std::stof(values[0]);
106  elementSpacing_[1] = std::stof(values[1]);
107  elementSpacing_[2] = std::stof(values[2]);
108  }
109  else if (elements[0].compare("DimSize") == 0)
110  {
111  std::vector<std::string> values;
112  split(elements[1], ' ', values);
113  dimSize_[0] = std::stoi(values[0]);
114  dimSize_[1] = std::stoi(values[1]);
115  dimSize_[2] = std::stoi(values[2]);
116  width_ = dimSize_[0];
117  height_ = dimSize_[1];
118  depth_ = dimSize_[2];
119  size_ = width_ * height_*depth_;
120  if(data_ == nullptr)
121  data_ = new T[dimSize_[0] * dimSize_[1] * dimSize_[2]];
122  }
123  else if (elements[0].compare("ElementDataFile") == 0)
124  {
125  full_path_to_file = full_path_to_file.substr(0, full_path_to_file.find_last_of("\\/"));
126  file_path_to_raw_file = full_path_to_file + '/' + elements[1];
127  }
128  }
129  }
130  }
131 
132  dataFile.close();
133  std::cout << "dimensions: " << dimSize_ << std::endl;
134  std::cout << "spacing: " << elementSpacing_ << std::endl;
135  std::cout << "offset: " << offset_ << std::endl;
136  std::cout << "transformMatrix: " << transformMatrix_ << std::endl;
137 
138  //- read raw file
139  std::ifstream dataFileRaw(file_path_to_raw_file, std::ios::in | std::ios::binary);
140 
141  if (dataFileRaw.is_open())
142  {
143  dataFileRaw.read((char*)data_, sizeof(T)*size_);
144  T distance = 0;
145  for(int index = 0; index<size_; index++)
146  {
147  distance = data_[index];
148  data_[index] = distance;
149  //std::cout <<index <<" "<< distance << '\n';
150  if (distance < min_value_) min_value_ = distance;
151  if (distance > max_value_) max_value_ = distance;
152  }
153  }
154  dataFileRaw.close();
155 
156  //write(std::string("sphere-binary"),ASCII);
157  }
158  template<typename T, int nDims>
159  ImageMHD<T, nDims>::ImageMHD(Real radius, Vec3i NxNyNz, Vec3d spacings):
160  objectType_("Image"),
161  ndims_(nDims),
162  binaryData_(true),
163  binaryDataByteOrderMSB_(false),
164  compressedData_(false),
165  transformMatrix_(Mat3d(1.0)),
166  offset_(Vec3d(-0.5*NxNyNz[0]*spacings[0], -0.5*NxNyNz[1] * spacings[1], -0.5*NxNyNz[2] * spacings[2])),
167  centerOfRotation_(Vec3d(0.0, 0.0, 0.0)),
168  elementSpacing_(spacings),
169  dimSize_(NxNyNz),
170  width_(dimSize_[0]),
171  height_(dimSize_[1]),
172  depth_(dimSize_[2]),
173  size_(width_*height_*depth_),
174  anatomicalOrientation_("???"),
175  elementType_(MET_FLOAT),
176  elementDataFile_(""),
177  min_value_(Infinity),
178  max_value_(-Infinity),
179  data_(nullptr)
180  {
181  if(data_ == nullptr)
182  data_ = new float[size_];
183 
184  Vec3d center(0.5*width_, 0.5*height_, 0.5*depth_);
185 
186  for (int z = 0; z < depth_; z++)
187  {
188  for (int y = 0; y < height_; y++)
189  {
190  for (int x = 0; x < width_; x++)
191  {
192  int index = z * width_*height_ + y * width_ + x;
193  double distance = (Vec3d(x, y, z) - center).norm() - radius;
194  if (distance < min_value_) min_value_ = distance;
195  if (distance > max_value_) max_value_ = distance;
196  data_[index] = float(distance);
197  }
198  }
199  }
200  write(std::string("sphere"), BINARY);
201  }
202 
203  template<typename T, int nDims>
204  ImageMHD<T, nDims>::~ImageMHD()
205  {
206  if (data_)
207  {
208  delete data_;
209  data_ = nullptr;
210  }
211  }
212 
213  //=================================================================================================//
214  template<typename T, int nDims>
215  std::vector<int> ImageMHD<T, nDims>::findNeighbors(const Vec3d& input_pnt, Vec3i& this_cell)
216  {
217  std::vector<int> neighbors;
218 
219  Vec3d image_coord = transformMatrix_.invert()*(input_pnt - offset_);
220  // std::cout <<"findNeighbor of " << input_pnt << " ........... " << image_coord << std::endl;
221 
222  int z = int(floor(image_coord[2]));
223  int y = int(floor(image_coord[1]));
224  int x = int(floor(image_coord[0]));
225 
226  //- cannot count cells in buffer zone
227  if (x < 0 || x > width_ - 1 || y < 0 || y > height_ - 1 || z < 0 || z > depth_ - 1) return neighbors;
228 
229  for (int k = z - 1; k < z + 2; k = k + 2)
230  {
231  for (int j = y - 1; j < y + 2; j = j + 2)
232  {
233  for (int i = x - 1; i < x + 2; i = i + 2)
234  {
235  if (i<0 || i >width_ - 1 || j <0 || j >height_ - 1 || k<0 || k >depth_)
236  continue;
237  int index = z * width_*height_ + y * width_ + x;
238  neighbors.push_back(index);
239  }
240  }
241 
242  }
243 
244  return neighbors;
245 
246  }
247  //=================================================================================================//
248  template<typename T, int nDims>
249  Vec3d ImageMHD<T, nDims>::computeGradientAtCell(int i)
250  {
251  //- translate 1D index to 3D index
252  int width = width_;
253  int height = height_;
254  int depth = depth_;
255  int sliceSize = width * height;
256  int z = i / sliceSize;
257  int y = (i%sliceSize) / width;
258  int x = (i%sliceSize) % width;
259 
260  double gradx = 0.0; double grady = 0.0; double gradz = 0.0;
261  //- cds (if inner cell)
262  //- otherwise back/forward scheme
263  if (x == 0)
264  {
265  int indexHigh = z * sliceSize + y * width + (x + 1);
266  gradx = (getValueAtCell(indexHigh) - getValueAtCell(i));
267  }
268  else if (x == width - 1)
269  {
270  int indexLow = z * sliceSize + y * width + (x - 1);
271  gradx = -(getValueAtCell(indexLow) - getValueAtCell(i));
272  }
273  else if (x > 0 && x < width_ - 1)
274  {
275  int indexHigh = z * sliceSize + y * width + (x + 1);
276  int indexLow = z * sliceSize + y * width + (x - 1);
277  gradx = (getValueAtCell(indexHigh) - getValueAtCell(indexLow)) / 2.0;
278  }
279 
280  if (y == 0)
281  {
282  int indexHigh = z * sliceSize + (y + 1)*width + x;
283  grady = (getValueAtCell(indexHigh) - getValueAtCell(i));
284  }
285  else if (y == height - 1)
286  {
287  int indexLow = z * sliceSize + (y - 1)*width + x;
288  grady = -(getValueAtCell(indexLow) - getValueAtCell(i));
289  }
290  else if (y > 0 && y < height_ - 1)
291  {
292  int indexHigh = z * sliceSize + (y + 1)*width + x;
293  int indexLow = z * sliceSize + (y - 1)*width + x;
294  grady = (getValueAtCell(indexHigh) - getValueAtCell(indexLow)) / 2.0;
295  }
296 
297  if (z == 0)
298  {
299  int indexHigh = (z + 1)*sliceSize + y * width + x;
300  gradz = (getValueAtCell(indexHigh) - getValueAtCell(i));
301  }
302  else if (z == depth - 1)
303  {
304  int indexLow = (z - 1)*sliceSize + y * width + x;
305  gradz = -(getValueAtCell(indexLow) - getValueAtCell(i));
306  }
307  else if (z > 0 && z < depth_ - 1)
308  {
309  int indexHigh = (z + 1)*sliceSize + y * width + x;
310  int indexLow = (z - 1)*sliceSize + y * width + x;
311  gradz = (getValueAtCell(indexHigh) - getValueAtCell(indexLow)) / 2.0;
312  }
313  gradx = gradx / elementSpacing_[0];
314  grady = grady / elementSpacing_[1];
315  gradz = gradz / elementSpacing_[2];
316  return Vec3d(gradx, grady, gradz);
317 
318  }
319  //=================================================================================================//
320  template<typename T, int nDims>
321  Vec3d ImageMHD<T, nDims>::computeNormalAtCell(int i)
322  {
323  Vec3d grad_phi = computeGradientAtCell(i);
324  Vec3d n = grad_phi.normalize();
325  return n;
326  }
327 
328  template<typename T, int nDims>
329  T ImageMHD<T, nDims>::getValueAtCell(int i)
330  {
331  if (i < 0 || i > size_)
332  {
333  return float(max_value_);
334  }
335  else
336  {
337  return data_[i];
338  }
339 
340  }
341  //=================================================================================================//
342  template<typename T, int nDims>
343  Vec3d ImageMHD<T, nDims>::convertToPhysicalSpace(Vec3d p)
344  {
345  Vec3d position = transformMatrix_ * p + offset_;
346  for (int i = 0; i < position.size(); i++)
347  {
348  position[i] = position[i] * elementSpacing_[i];
349  }
350  return position;
351 
352  }
353  //=================================================================================================//
354  template<typename T, int nDims>
355  void ImageMHD<T, nDims>::split(const std::string &s, char delim, \
356  std::vector<std::string> &elems)
357  {
358  std::stringstream ss(s);
359  std::string item;
360  while (std::getline(ss, item, delim)) {
361  if (item.length() > 0) elems.push_back(item);
362  }
363  }
364  //=================================================================================================//
365  template<typename T, int nDims>
366  Vec3d ImageMHD<T, nDims>::findClosestPoint(const Vec3d& input_pnt)
367  {
368  Vec3i this_cell;
369  std::vector<int> neighbors = findNeighbors(input_pnt, this_cell);
370  Vec3d n_sum(0.0, 0.0, 0.0);
371  double weight_sum = 0.0;
372  double d_sum = 0.0;
373  for (const int& i : neighbors)
374  {
375  // checkIndexBound(i);
376  Vec3d nCj = computeNormalAtCell(i);
377  double dCj = float(getValueAtCell(i));
378  double weight_Cj = 1.0 / (fabs(dCj) + Eps);
379  n_sum = n_sum + weight_Cj * nCj;
380  weight_sum = weight_sum + weight_Cj;
381  d_sum = d_sum + dCj;
382  }
383  Vec3d n = n_sum / (weight_sum + Eps);
384  double d = d_sum / (weight_sum + Eps);
385 
386  Vec3d p_image = Vec3d(this_cell[0], this_cell[1], this_cell[2]) + n.normalize()*d;
387  Vec3d p = convertToPhysicalSpace(p_image);
388  return p;
389 
390  }
391 
392  template<typename T, int nDims>
393  BoundingBox ImageMHD<T, nDims>::findBounds()
394  {
395  //initial reference values
396  Vec3d lower_bound = Vec3d(Infinity);
397  Vec3d upper_bound = Vec3d(-Infinity);
398 
399  for (int z = 0; z < depth_ + 1; z++)
400  {
401  for (int y = 0; y < height_ + 1; y++)
402  {
403  for (int x = 0; x < width_ + 1; x++)
404  {
405  Vec3d p_image = Vec3d(x, y, z);
406  Vec3d vertex_position = convertToPhysicalSpace(p_image);
407  for (int j = 0; j != 3; ++j) {
408  lower_bound[j] = SMIN(lower_bound[j], vertex_position[j]);
409  upper_bound[j] = SMAX(upper_bound[j], vertex_position[j]);
410  }
411  }
412  }
413  }
414  return BoundingBox(lower_bound, upper_bound);
415  }
416 
417  template<typename T, int nDims>
418  Real ImageMHD<T, nDims>::findValueAtPoint(const Vec3d& input_pnt)
419  {
420  Vec3i this_cell;
421  std::vector<int> neighbors = findNeighbors(input_pnt, this_cell);
422  double weight_sum = 0.0;
423  double d_sum = 0.0;
424  if (neighbors.size() > 0)
425  {
426  for (const int& i : neighbors)
427  {
428  // checkIndexBound(i);
429  double dCj = float(getValueAtCell(i));
430  double weight_Cj = 1.0 / (fabs(dCj) + Eps);
431  weight_sum = weight_sum + weight_Cj;
432  d_sum = d_sum + dCj;
433  }
434  return d_sum / (weight_sum + Eps);
435  }
436  else
437  {
438  return max_value_;
439  }
440 
441  }
442  //=================================================================================================//
443  template<typename T, int nDims>
444  Vec3d ImageMHD<T, nDims>::findNormalAtPoint(const Vec3d & input_pnt)
445  {
446  Vec3i this_cell;
447  std::vector<int> neighbors = findNeighbors(input_pnt, this_cell);
448  Vec3d n_sum(0.0, 0.0, 0.0);
449  double weight_sum = 0.0;
450  double d_sum = 0.0;
451  if (neighbors.size() > 0)
452  {
453  for (const int& i : neighbors)
454  {
455  // checkIndexBound(i);
456  Vec3d nCj = computeNormalAtCell(i);
457  double dCj = float(getValueAtCell(i));
458  double weight_Cj = 1.0 / (fabs(dCj) + Eps);
459  n_sum = n_sum + weight_Cj * nCj;
460  weight_sum = weight_sum + weight_Cj;
461  d_sum = d_sum + dCj;
462  }
463  Vec3d n = n_sum / (weight_sum + Eps);
464  return n.normalize();
465  }
466  else
467  {
468  return Vec3d(1.0, 1.0, 1.0).normalize();
469  }
470  }
471 
472  //=================================================================================================//
473  template<typename T, int nDims>
474  void ImageMHD<T, nDims>::write(std::string filename, Output_Mode mode)
475  {
476  std::ofstream output_file(filename+".mhd", std::ofstream::out);
477  output_file << "ObjectType = " << objectType_ << "\n";
478  output_file << "NDims = " << ndims_ << "\n";
479  if (mode == BINARY)
480  output_file << "BinaryData = True" << "\n";
481  else
482  output_file << "BinaryData = False" << "\n";
483  output_file << "BinaryDataByteOrderMSB = " << binaryDataByteOrderMSB_ << "\n";
484  output_file << "CompressedData = " << compressedData_ << "\n";
485  output_file << "TransformMatrix = "
486  << transformMatrix_[0][0] << " " << transformMatrix_[0][1] << " " << transformMatrix_[0][2] << " "
487  << transformMatrix_[1][0] << " " << transformMatrix_[1][1] << " " << transformMatrix_[1][2] << " "
488  << transformMatrix_[2][0] << " " << transformMatrix_[2][1] << " " << transformMatrix_[2][2] << "\n";
489  output_file << "Offset = "
490  << offset_[0] << " " << offset_[1] << " " << offset_[2] << "\n";
491  output_file << "CenterOfRotation = "
492  << centerOfRotation_[0] << " " << centerOfRotation_[1] << " " << centerOfRotation_[2] << "\n";
493  output_file << "ElementSpacing = "
494  << elementSpacing_[0] << " " << elementSpacing_[1] << " " << elementSpacing_[2] << "\n";
495  output_file << "DimSize = "
496  << dimSize_[0] << " " << dimSize_[1] << " " << dimSize_[2] << "\n";
497  output_file << "AnatomicalOrientation = " << anatomicalOrientation_ << "\n";
498  if (elementType_ == MET_FLOAT)
499  output_file << "ElementType = MET_FLOAT" << "\n";
500  else if (elementType_ == MET_UCHAR)
501  output_file << "ElementType = MET_UCHAR" << "\n";
502  if (elementType_ == MET_LONG)
503  output_file << "ElementType = MET_LONG" << "\n";
504  output_file << "ElementDataFile = "<< filename+".raw" << "\n";
505 
506  output_file.close();
507 
508  if (mode == BINARY)
509  {
510  std::ofstream output_file_raw(filename + ".raw", std::ios::binary | std::ios::out);
511  output_file_raw.write((const char*)data_, sizeof(T)*size_);
512  output_file_raw.close();
513  }
514  else
515  {
516  std::ofstream output_file_raw(filename + ".raw");
517  for (int index = 0; index < size_; index++)
518  {
519  output_file_raw << data_[index] << std::endl;
520  }
521  output_file_raw.close();
522  }
523 
524 
525  }
526 }
527 #endif //__EMSCRIPTEN__
528 
529 #endif //IMAGE_MHD_HPP
SPH
Definition: solid_body_supplementary.cpp:9
Generated by   doxygen 1.8.14