OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
main.cpp
Go to the documentation of this file.
1 //
2 // Created by jstark on 2020-05-18.
3 //
45 
47 // Include redistancing header files
48 #include "util/PathsAndFiles.hpp"
53 
54 
74 int main(int argc, char* argv[])
76 {
77  typedef double phi_type;
78  // initialize library
79  openfpm_init(&argc, &argv);
81  // Set current working directory, define output paths and create folders where output will be saved
82  std::string cwd = get_cwd();
83  const std::string path_output = cwd + "/output_3D_sphere";
84  create_directory_if_not_exist(path_output);
85 
87  // Now we set the input paths. We need a binary file with the pixel values and a csv file with the
88  // size of the stack (in #pixels / dimension)
89  const std::string path_input ="input";
90  const std::string path_to_zstack = path_input + "/sphere.bin";
91  const std::string path_to_size = path_input + "/size_sphere.csv";
92  /*
93  * in case of 3D (image volume):
94  */
95  const unsigned int grid_dim = 3;
96  const size_t x = 0;
97  const size_t y = 1;
98  const size_t z = 2;
99 
100  const size_t Phi_0_grid = 0;
101  const size_t Phi_SDF_grid = 1;
102 
103  const size_t Phi_SDF_vd = 0;
104  const size_t Phi_grad_vd = 1;
105  const size_t Phi_magnOfGrad_vd = 2;
107 
108 
123  const phi_type refinement [] = {1.0, 1.0, 1.0}; // without refinement
125 // const phi_type refinement [] = {0.8, 1.5, 2.0}; // factors by which grid should be finer as underlying image stack in each dimension (e.g. to get isotropic grid from anisotropic stack resolution)
128  // read the stack size (number of pixel values per dimension) from a binary file
129  // alternatively, you can directly define the stack size as: std::vector<size_t> stack_size {#pixels in x, #pixels in y, #pixels in z}
144  std::vector<int> stack_size = get_size(path_to_size);
146  auto & v_cl = create_vcluster();
147  if (v_cl.rank() == 0)
148  {
149  for(std::vector<int>::size_type i = 0; i != stack_size.size(); i++)
150  {
151  std::cout << "#Pixel in dimension " << i << " = " << stack_size[i] << std::endl;
152  std::cout << "original refinement in dimension " << i << " = " << refinement[i] << std::endl;
153  }
154  }
155 
156 
157  // Array with grid dimensions after refinement. This size-array will be used for the grid creation.
158  const size_t sz[grid_dim] = {(size_t)std::round(stack_size[x] * refinement[x]), (size_t)std::round(stack_size[y] * refinement[y]), (size_t)std::round(stack_size[z] * refinement[z])}; // 3D
160  // Here we create a 3D grid that stores 2 properties:
161  // Prop1: store the initial Phi;
162  // Prop2: here the re-initialized Phi (signed distance function) will be written to in the re-distancing step
176  Box<grid_dim, phi_type> box({0.0, 0.0, 0.0}, {5.0, 5.0, 5.0}); // 3D
178 
179  Ghost<grid_dim, long int> ghost(0);
180  typedef aggregate<phi_type, phi_type> props;
181  typedef grid_dist_id<grid_dim, phi_type, props > grid_in_type;
182  grid_in_type g_dist(sz, box, ghost);
183  g_dist.setPropNames({"Phi_0", "Phi_SDF"});
184 
185  // initialize complete grid including ghost layer with -1
186  init_grid_and_ghost<Phi_0_grid>(g_dist, -1.0);
187 
188  // Now we can initialize the grid with the pixel values from the image stack
189  load_pixel_onto_grid<Phi_0_grid>(g_dist, path_to_zstack, stack_size);
190  g_dist.write(path_output + "/grid_from_images_initial", FORMAT_BINARY); // Save the initial grid as vtk file
191 
192 
194  // Now we want to convert the initial Phi into a signed distance function (SDF) with magnitude of gradient = 1
195  // For the initial re-distancing we use the Sussman method
196  // 1.) Set some redistancing options (for details see example sussman disk or sphere)
197  Redist_options<phi_type> redist_options;
198  redist_options.min_iter = 1e3;
199  redist_options.max_iter = 1e4;
200 
201  redist_options.convTolChange.value = 1e-7;
202  redist_options.convTolChange.check = true;
203  redist_options.convTolResidual.value = 1e-6; // is ignored if convTolResidual.check = false;
204  redist_options.convTolResidual.check = false;
205 
206  redist_options.interval_check_convergence = 1e3;
207  redist_options.width_NB_in_grid_points = 10;
208  redist_options.print_current_iterChangeResidual = true;
209  redist_options.print_steadyState_iter = true;
210  redist_options.save_temp_grid = true;
211  RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options); // Instantiation of Sussman-redistancing
212  // class
213 // std::cout << "dt = " << redist_obj.get_time_step() << std::endl;
214  // Run the redistancing. in the <> brackets provide property-index where 1.) your initial Phi is stored and 2.) where the resulting SDF should be written to.
215  redist_obj.run_redistancing<Phi_0_grid, Phi_SDF_grid>();
216 
217  g_dist.write(path_output + "/grid_images_post_redistancing", FORMAT_BINARY); // Save the result as vtk file
218 
220  // Get narrow band: Place particles on interface (narrow band width e.g. 2 grid points on each side of the interface)
221  size_t bc[grid_dim] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
222  // Create an empty vector to which narrow-band particles will be added. You can choose, how many properties you want.
223  // Minimum is 1 property, to which the Phi_SDF can be written
224  // In this example we chose 3 properties. The 1st for the Phi_SDF, the 2nd for the gradient of phi and the 3rd for
225  // the magnitude of the gradient
226  typedef aggregate<phi_type, Point<grid_dim, phi_type>, phi_type> props_nb;
228  Ghost<grid_dim, phi_type> ghost_vd(0);
229  vd_type vd_narrow_band(0, box, bc, ghost_vd);
230  vd_narrow_band.setPropNames({"Phi_SDF", "Phi_grad", "Phi_magnOfGrad"});
231 
232  NarrowBand<grid_in_type, phi_type> narrowBand(g_dist, redist_options.width_NB_in_grid_points); // Instantiation of
233  // NarrowBand class
234 
235  // Get the narrow band. You can decide, if you only want the Phi_SDF saved to your particles or
236  // if you also want the gradients or gradients and magnitude of gradient.
237  // The function will know depending on how many property-indices you provide in the <> brackets.
238  // First property-index must always be the index of the SDF on the grid!
239  // E.g.: The following line would only write only the Phi_SDF from the grid to your narrow-band particles
240  // narrowBand.get_narrow_band<Phi_SDF_grid, Phi_SDF_vd>(g_dist, vd_narrow_band);
241  // Whereas this will give you the gradients and magnOfGrad of Phi as well:
242  narrowBand.get_narrow_band<Phi_SDF_grid, Phi_SDF_vd, Phi_grad_vd, Phi_magnOfGrad_vd>(g_dist, vd_narrow_band);
243 
244  vd_narrow_band.write(path_output + "/vd_narrow_band_images", FORMAT_BINARY); // Save particles as vtk file
245 
246  openfpm_finalize(); // Finalize openFPM library
247  return 0;
248 }
250 
static void create_directory_if_not_exist(std::string path)
Creates a directory if not already existent.
Class for getting the narrow band around the interface.
Definition: NarrowBand.hpp:42
Class for reinitializing a level-set function into a signed distance function using Sussman redistanc...
std::vector< int > get_size(const std::string &path_to_file)
Read the number of pixels per dimension from a csv-file in order to create a grid with the same size.
Structure to bundle options for redistancing.
static std::string get_cwd()
Gets the current working directory and returns path as string.
Definition: Ghost.hpp:39
Class for reinitializing a level-set function into a signed distance function using Sussman redistanc...
This is a distributed grid.
This class represent an N-dimensional box.
Definition: Box.hpp:60
Header file containing functions for loading pixel values from 2D image or 3D image stack (volume) st...
Distributed vector.
Header file containing functions for creating files and folders.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:214
Class for getting the narrow band around the interface.