main.cpp

Go to the documentation of this file.

//
// Created by jstark on 2020-05-18.
//

// Include redistancing header files
#include "util/PathsAndFiles.hpp"
#include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
#include "level_set/redistancing_Sussman/NarrowBand.hpp"
#include "RawReader/InitGridWithPixel.hpp"


int main(int argc, char* argv[])
{
    typedef double phi_type;
    //  initialize library
    openfpm_init(&argc, &argv);
    // Set current working directory, define output paths and create folders where output will be saved
    std::string cwd                     = get_cwd();
    const std::string path_output       = cwd + "/output_3D_sphere";
    create_directory_if_not_exist(path_output);
    
    // Now we set the input paths. We need a binary file with the pixel values and a csv file with the
    // size of the stack (in #pixels / dimension)
    const std::string path_input        ="input";
    const std::string path_to_zstack    = path_input + "/sphere.bin";
    const std::string path_to_size      = path_input + "/size_sphere.csv";
    /*
     * in case of 3D (image volume):
     */
    const unsigned int grid_dim         = 3;
    const size_t x                      = 0;
    const size_t y                      = 1;
    const size_t z                      = 2;

    const size_t Phi_0_grid             = 0;
    const size_t Phi_SDF_grid           = 1;

    const size_t Phi_SDF_vd             = 0;
    const size_t Phi_grad_vd            = 1;
    const size_t Phi_magnOfGrad_vd      = 2;
    
    
    const phi_type refinement []          = {1.0, 1.0, 1.0};      // without refinement
//  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)
    //  read the stack size (number of pixel values per dimension) from a binary file
    // alternatively, you can directly define the stack size as: std::vector<size_t> stack_size {#pixels in x, #pixels in y, #pixels in z}
    std::vector<int> stack_size = get_size(path_to_size);
    auto & v_cl = create_vcluster();
    if (v_cl.rank() == 0)
    {
        for(std::vector<int>::size_type i = 0; i != stack_size.size(); i++)
        {
            std::cout << "#Pixel in dimension " << i << " = " << stack_size[i] << std::endl;
            std::cout << "original refinement in dimension " << i << " = " << refinement[i] << std::endl;
        }
    }
    
    
    //  Array with grid dimensions after refinement. This size-array will be used for the grid creation.
    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
    // Here we create a 3D grid that stores 2 properties:
    // Prop1: store the initial Phi;
    // Prop2: here the re-initialized Phi (signed distance function) will be written to in the re-distancing step
    Box<grid_dim, phi_type> box({0.0, 0.0, 0.0}, {5.0, 5.0, 5.0}); // 3D

    Ghost<grid_dim, long int> ghost(0);
    typedef aggregate<phi_type, phi_type> props;
    typedef grid_dist_id<grid_dim, phi_type, props > grid_in_type;
    grid_in_type g_dist(sz, box, ghost);
    g_dist.setPropNames({"Phi_0", "Phi_SDF"});
    
    // initialize complete grid including ghost layer with -1
    init_grid_and_ghost<Phi_0_grid>(g_dist, -1.0);
    
    // Now we can initialize the grid with the pixel values from the image stack
    load_pixel_onto_grid<Phi_0_grid>(g_dist, path_to_zstack, stack_size);
    g_dist.write(path_output + "/grid_from_images_initial", FORMAT_BINARY); // Save the initial grid as vtk file


    // Now we want to convert the initial Phi into a signed distance function (SDF) with magnitude of gradient = 1
    // For the initial re-distancing we use the Sussman method
    // 1.) Set some redistancing options (for details see example sussman disk or sphere)
    Redist_options<phi_type> redist_options;
    redist_options.min_iter                             = 1e3;
    redist_options.max_iter                             = 1e4;
    
    redist_options.convTolChange.value                  = 1e-7;
    redist_options.convTolChange.check                  = true;
    redist_options.convTolResidual.value                = 1e-6; // is ignored if convTolResidual.check = false;
    redist_options.convTolResidual.check                = false;
    
    redist_options.interval_check_convergence           = 1e3;
    redist_options.width_NB_in_grid_points              = 10;
    redist_options.print_current_iterChangeResidual     = true;
    redist_options.print_steadyState_iter               = true;
    redist_options.save_temp_grid                       = true;
    RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options);   // Instantiation of Sussman-redistancing
    // class
//  std::cout << "dt = " << redist_obj.get_time_step() << std::endl;
    // 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.
    redist_obj.run_redistancing<Phi_0_grid, Phi_SDF_grid>();
    
    g_dist.write(path_output + "/grid_images_post_redistancing", FORMAT_BINARY); // Save the result as vtk file
    
    //  Get narrow band: Place particles on interface (narrow band width e.g. 2 grid points on each side of the interface)
    size_t bc[grid_dim] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
    // Create an empty vector to which narrow-band particles will be added. You can choose, how many properties you want.
    // Minimum is 1 property, to which the Phi_SDF can be written
    // In this example we chose 3 properties. The 1st for the Phi_SDF, the 2nd for the gradient of phi and the 3rd for
    // the magnitude of the gradient
    typedef aggregate<phi_type, Point<grid_dim, phi_type>, phi_type> props_nb;
    typedef vector_dist<grid_dim, phi_type, props_nb> vd_type;
        Ghost<grid_dim, phi_type> ghost_vd(0);
        vd_type vd_narrow_band(0, box, bc, ghost_vd);   
        vd_narrow_band.setPropNames({"Phi_SDF", "Phi_grad", "Phi_magnOfGrad"});
    
    NarrowBand<grid_in_type, phi_type> narrowBand(g_dist, redist_options.width_NB_in_grid_points); // Instantiation of
    // NarrowBand class
    
    // Get the narrow band. You can decide, if you only want the Phi_SDF saved to your particles or
    // if you also want the gradients or gradients and magnitude of gradient.
    // The function will know depending on how many property-indices you provide in the <> brackets.
    // First property-index must always be the index of the SDF on the grid!
    // E.g.: The following line would only write only the Phi_SDF from the grid to your narrow-band particles
    // narrowBand.get_narrow_band<Phi_SDF_grid, Phi_SDF_vd>(g_dist, vd_narrow_band);
    // Whereas this will give you the gradients and magnOfGrad of Phi as well:
    narrowBand.get_narrow_band<Phi_SDF_grid, Phi_SDF_vd, Phi_grad_vd, Phi_magnOfGrad_vd>(g_dist, vd_narrow_band);
    
    vd_narrow_band.write(path_output + "/vd_narrow_band_images", FORMAT_BINARY); // Save particles as vtk file

    openfpm_finalize(); // Finalize openFPM library
    return 0;
}