voxel_map.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/util/assertion.hpp>
#include <kernel/util/dist.hpp>
#include <kernel/util/dist_file_io.hpp>
#include <kernel/geometry/cgal.hpp>
#include <kernel/geometry/atlas/chart.hpp>
 
// includes, system
#include <array>
#include <vector>
#include <iostream>
#include <cstdint>
#include <cstring>
 
// includes, thirdparty
#ifdef FEAT_HAVE_FPARSER
#include <fparser.hh>
#endif // FEAT_HAVE_FPARSER
 
namespace FEAT
{
  namespace Geometry
  {
    class VoxelMapFileError :
      public FileError
    {
    public:
      explicit VoxelMapFileError(const String& filename, const String& msg) :
        FileError(filename.empty() ? msg : (filename + ": " + msg))
      {
      }
    }; // class VoxelMapFileError
 
    class VoxelMapFormulaParseError :
      public ParseError
    {
    public:
      explicit VoxelMapFormulaParseError(const String& msg) :
        ParseError(msg)
      {
      }
    }; // class VoxelMapFormulaParseError
 
    template<typename CoordType_, int dim_>
    class VoxelMasker
    {
    public:
      typedef CoordType_ CoordType;
      typedef Tiny::Vector<CoordType_, dim_> PointType;
 
      virtual ~VoxelMasker()
      {
      }
 
      static CoordType x_coord(const CoordType x_min, const CoordType x_max, const std::size_t i, const std::size_t n)
      {
        return x_min + (x_max - x_min) * CoordType(i) / CoordType(n - 1u);
      }
 
      virtual void mask_line(std::vector<int>& mask, const CoordType x_min, const CoordType x_max, const PointType& point) = 0;
    }; // class VoxelMasker<...>
 
#if defined(FEAT_HAVE_CGAL) || defined(DOXYGEN)
    template<typename CoordType_>
    class VoxelCGALMasker :
      public VoxelMasker<CoordType_, 3>
    {
    public:
      typedef VoxelMasker<CoordType_, 3> BaseClass;
      using typename BaseClass::CoordType;
      using typename BaseClass::PointType;
 
    private:
      const Geometry::CGALWrapper<CoordType_>& _cgal_wrapper;
      const bool _invert;
 
    public:
      explicit VoxelCGALMasker(const Geometry::CGALWrapper<CoordType_>& cgal_wrapper, bool invert) :
        _cgal_wrapper(cgal_wrapper),
        _invert(invert)
      {
      }
 
      virtual void mask_line(std::vector<int>& mask, const CoordType x_min, const CoordType x_max, const PointType& point) override
      {
        if(mask.empty())
          return;
 
        const std::size_t nv = mask.size();
        for(std::size_t i(0); i < nv; ++i)
        {
          mask[i] = (this->_cgal_wrapper.point_inside(BaseClass::x_coord(x_min, x_max, i, nv),
            point[1], point[2]) ? int(!_invert) : int(_invert));
        }
      }
    }; // class VoxelCGALMasker<...>
#endif // FEAT_HAVE_CGAL || DOXYGEN
 
    template<typename MeshType>
    class VoxelChartMasker :
      public VoxelMasker<typename MeshType::CoordType, MeshType::shape_dim>
    {
    public:
      typedef VoxelMasker<typename MeshType::CoordType, MeshType::shape_dim> BaseClass;
      using typename BaseClass::CoordType;
      using typename BaseClass::PointType;
 
    private:
      const Geometry::Atlas::ChartBase<MeshType>& _chart;
      const bool _invert;
 
    public:
      explicit VoxelChartMasker(const Geometry::Atlas::ChartBase<MeshType>& chart, bool invert) :
        _chart(chart),
        _invert(invert)
      {
      }
 
      virtual void mask_line(std::vector<int>& mask, const CoordType x_min, const CoordType x_max, const PointType& point) override
      {
        if(mask.empty())
          return;
 
        PointType pt(point);
 
        const std::size_t nv = mask.size();
        for(std::size_t i(0); i < nv; ++i)
        {
          pt[0] = BaseClass::x_coord(x_min, x_max, i, nv);
          mask[i] = (_chart.signed_dist(pt) <= CoordType(0) ? int(!_invert) : int(_invert));
        }
      }
    }; // class VoxelChartMasker<...>
 
    template<typename Lambda_, typename CoordType_, int dim_>
    class VoxelLambdaMasker :
      public VoxelMasker<CoordType_, dim_>
    {
    public:
      typedef VoxelMasker<CoordType_, dim_> BaseClass;
      using typename BaseClass::CoordType;
      using typename BaseClass::PointType;
 
    public:
      Lambda_ _lambda;
 
      explicit VoxelLambdaMasker(Lambda_&& lambda) :
        _lambda(std::forward<Lambda_>(lambda))
      {
      }
 
      virtual void mask_line(std::vector<int>& mask, const CoordType x_min, const CoordType x_max, const PointType& point) override
      {
        if(mask.empty())
          return;
 
        PointType pt(point);
        // define a const reference to ensure that the lambda expression does not modify the point
        const PointType&  cpt(pt);
 
        const std::size_t nv = mask.size();
        for(std::size_t i(0); i < nv; ++i)
        {
          pt[0] = BaseClass::x_coord(x_min, x_max, i, nv);
          mask[i] = (_lambda(cpt) ? 1 : 0);
        }
      }
    }; // class VoxelLambdaMasker<...>
 
#if defined(FEAT_HAVE_FPARSER) || defined(DOXYGEN)
    template<int dim_>
    class VoxelFormulaMasker :
      public VoxelMasker<double, dim_>
    {
    public:
      typedef VoxelMasker<double, dim_> BaseClass;
      using typename BaseClass::CoordType;
      using typename BaseClass::PointType;
 
    public:
      String _formula;
      ::FunctionParser _parser;
 
      explicit VoxelFormulaMasker(const String& formula) :
        _formula(formula),
        _parser()
      {
        _parser.AddConstant("pi", Math::pi<double>());
        _parser.AddConstant("eps", Math::eps<double>());
 
        String vars("x");
        if(dim_ > 1) vars += ",y";
        if(dim_ > 2) vars += ",z";
 
        // try to parse the function
        const int ret = _parser.Parse(_formula.c_str(), vars.c_str());
        if(ret >= 0)
        {
          String msg(_parser.ErrorMsg());
          msg.append("\n>>> '");
          msg.append(_formula);
          msg.append("'");
          throw VoxelMapFormulaParseError(msg);
        }
 
        // optimize the parsed function
        _parser.Optimize();
      }
 
      virtual void mask_line(std::vector<int>& mask, const CoordType x_min, const CoordType x_max, const PointType& point) override
      {
        if(mask.empty())
          return;
 
        PointType pt(point);
 
        const std::size_t nv = mask.size();
        for(std::size_t i(0); i < nv; ++i)
        {
          pt[0] = BaseClass::x_coord(x_min, x_max, i, nv);
          mask[i] = (_parser.Eval(pt.v) > 0.0 ? 1 : 0);
        }
      }
    }; // class VoxelFormulaMasker<...>
#endif // FEAT_HAVE_FPARSER || DOXYGEN
 
    class VoxelMap
    {
    public:
      typedef std::int64_t i64;
 
      typedef std::uint64_t u64;
 
      static constexpr u64 magic_number = 0x70614D6C65786F56ull; // "VoxelMap"
 
      static constexpr u64 header_size = 136ull;
 
      static constexpr i64 unit_size = 1'000'000'000;
 
      struct FileHeader
      {
        u64 magic_number;     //< "VoxelMap" = 0x70614D6C65786F56ll
        u64 header_size;      //< 136 bytes
        i64 min_x;            //< minimum X-coordinate in nanometer
        i64 max_x;            //< maximum X-coordinate in nanometer
        u64 num_x;            //< number of points in X dimension
        u64 stride_line;      //< number of bytes for a single uncompressed X-line; a multiple of 16; usually stride_line = (((num_x + 7ull)/8ull + 15ull) & ~15ull)
        i64 min_y;            //< minimum Y-coordinate in nanometer
        i64 max_y;            //< maximum Y-coordinate in nanometer
        u64 num_y;            //< number of points in Y dimension
        u64 stride_plane;     //< number of bytes for a single uncompressed XY-plane; a multiple of 16; usually stride_plane = stride_line * num_y
        i64 min_z;            //< minimum Z-coordinate in nanometer
        i64 max_z;            //< maximum Z-coordinate in nanometer
        u64 num_z;            //< number of points in Z dimension
        u64 stride_volume;    //< number of bytes for a single uncompressed XYZ-volume; a multiple of 16; usually stride_volume = stride_plane * num_z
        u64 coverage;         //< coverage of domain in voxel map, relative to unit size
        u64 planes_per_block; //< number of planes per compression block (0 if uncompressed)
        u64 num_blocks;       //< total number of compression blocks (0 if uncompressed)
      }; // struct FileHeader: 136 bytes
 
      static_assert(sizeof(FileHeader) == header_size, "invalid FileHeader size");
 
      static u64 calc_line_stride(u64 num_x)
      {
        // compute number of bytes by rounding to next multiple of 8 then dividing by 8,
        // afterwards round up to next multiple of 16
        return (((num_x + 7ull)/8ull + 15ull) & ~15ull);
      }
 
      class WriteResult
      {
      public:
        u64 filesize;            //< total size of written file in bytes
        u64 compress_block_size; //< chosen maximum compression block size in MB
        u64 compress_level;      //< chosen compression level
        u64 num_compress;        //< number of compression blocks
        u64 compressed_size;     //< total compressed buffer size in bytes
 
        WriteResult() :
          filesize(0u),
          compress_block_size(0u),
          compress_level(0u),
          num_compress(0u),
          compressed_size(0u)
        {
        }
 
        explicit WriteResult(u64 fs, u64 cbs = 0u, u64 cl = 0u, u64 nc = 0u, u64 cs = 0u) :
          filesize(fs),
          compress_block_size(cbs),
          compress_level(cl),
          num_compress(nc),
          compressed_size(cs)
        {
        }
 
        operator bool() const
        {
          return filesize > u64(0);
        }
      }; // class WriteResult
 
      class ReadResult
      {
      public:
        u64 filesize;            //< total size of read file in bytes
        u64 num_compress;        //< number of compression blocks
        u64 compressed_size;     //< total compressed buffer size in bytes
 
        ReadResult() :
          filesize(0u),
          num_compress(0u),
          compressed_size(0u)
        {
        }
 
        explicit ReadResult(u64 fs, u64 nc = 0u, u64 cs = 0u) :
          filesize(fs),
          num_compress(nc),
          compressed_size(cs)
        {
        }
 
        operator bool() const
        {
          return filesize > u64(0);
        }
      }; // class ReadResult
 
    protected:
      int _create_stage;
      std::array<i64, 3u> _bbox_min, _bbox_max;
      std::array<u64, 3u> _num_points;
      u64 _stride_line;
      u64 _stride_plane;
      u64 _num_lines;
      u64 _num_planes;
      std::vector<char> _voxel_map;
      bool _out_of_bounds_value;
      u64 _coverage;
 
    public:
      VoxelMap();
 
      VoxelMap(const VoxelMap&) = delete;
      VoxelMap& operator=(const VoxelMap&) = delete;
 
      VoxelMap(VoxelMap&& other);
 
      VoxelMap& operator=(VoxelMap&& other);
 
      virtual ~VoxelMap();
 
      void set_bounding_box_2d(Real x_min, Real x_max, Real y_min, Real y_max);
 
      void set_bounding_box_3d(Real x_min, Real x_max, Real y_min, Real y_max, Real z_min, Real z_max);
 
      Real get_bounding_box_min(int dim) const
      {
        XASSERT((dim >= 0) && (dim < 3));
        return Real(this->_bbox_min[std::size_t(dim)]) / Real(unit_size);
      }
 
      Real get_bounding_box_max(int dim) const
      {
        XASSERT((dim >= 0) && (dim < 3));
        return Real(this->_bbox_max[std::size_t(dim)]) / Real(unit_size);
      }
 
      Real get_bounding_box_volume() const;
 
      void set_out_of_bounds_value(bool value)
      {
        this->_out_of_bounds_value = value;
      }
 
      bool get_out_of_bounds_value() const
      {
        return this->_out_of_bounds_value;
      }
 
      void set_num_points(Index num_x, Index num_y, Index num_z = 0u);
 
      Index get_num_points(int dim) const
      {
        return Index(this->_num_points.at(Index(dim)));
      }
 
      void set_resolution(Real max_res);
 
      Index get_num_voxels() const
      {
        return Index(this->_num_points[0] * this->_num_points[1] * this->_num_points[2]);
      }
 
      u64 get_stride_line() const
      {
        return this->_stride_line;
      }
 
      u64 get_stride_plane() const
      {
        return this->_stride_plane;
      }
 
      u64 get_stride_volume() const
      {
        return u64(this->_voxel_map.size());
      }
 
      const std::vector<char>& get_map() const
      {
        return this->_voxel_map;
      }
 
      Real get_domain_coverage() const
      {
        // we store the coverage relative to unit size
        return Real(this->_coverage) / Real(unit_size);
      }
 
      Real get_domain_volume() const
      {
        return get_domain_coverage() * get_bounding_box_volume();
      }
 
      template<typename Coord_, int dim_>
      void compute_map(const Dist::Comm& comm, VoxelMasker<Coord_, dim_>& masker, bool gather_to_all = true)
      {
        XASSERT(!_voxel_map.empty());
        this->_compute_voxel_map(comm, masker, gather_to_all);
      }
 
      template<typename MeshType_>
      void compute_map_from_chart(const Dist::Comm& comm, const Geometry::Atlas::ChartBase<MeshType_>& chart, bool invert, bool gather_to_all = true)
      {
        XASSERT(!_voxel_map.empty());
        VoxelChartMasker<MeshType_> masker(chart, invert);
        this->_compute_voxel_map(comm, masker, gather_to_all);
      }
 
      template<typename Lambda_>
      void compute_map_from_lambda_2d(const Dist::Comm& comm, Lambda_&& lambda, bool gather_to_all = true)
      {
        VoxelLambdaMasker<Lambda_, Real, 2> masker(std::forward<Lambda_>(lambda));
        this->_compute_voxel_map(comm, masker, gather_to_all);
      }
 
      template<typename Lambda_>
      void compute_map_from_lambda_3d(const Dist::Comm& comm, Lambda_&& lambda, bool gather_to_all = true)
      {
        VoxelLambdaMasker<Lambda_, Real, 3> masker(std::forward<Lambda_>(lambda));
        this->_compute_voxel_map(comm, masker, gather_to_all);
      }
 
      void compute_map_from_formula_2d(const Dist::Comm& comm, const String& formula, bool gather_to_all = true);
 
      void compute_map_from_formula_3d(const Dist::Comm& comm, const String& formula, bool gather_to_all = true);
 
      void compute_map_from_off_3d(const Dist::Comm& comm, const String& filename, bool invert, bool gather_to_all = true);
 
      void compute_map_from_off_3d(const Dist::Comm& comm, std::istream& is, bool invert, bool gather_to_all = true);
 
      WriteResult write(const String& filename, const u64 compress_block_size = 128u, const int compress_level = 9) const;
 
      WriteResult write(std::ostream& os, const u64 compress_block_size = 128u, const int compress_level = 9) const;
 
      ReadResult read(const Dist::Comm& comm, const String& filename);
 
      ReadResult read(std::istream& is, const String& filename = "");
 
      bool check_point(u64 ix, u64 iy = 0u, u64 iz = 0u) const
      {
        ASSERT((ix < this->_num_points[0]) || ((this->_num_points[0] == 0u) && (ix == 0u)));
        ASSERT((iy < this->_num_points[1]) || ((this->_num_points[1] == 0u) && (iy == 0u)));
        ASSERT((iz < this->_num_points[2]) || ((this->_num_points[2] == 0u) && (iz == 0u)));
        return this->_check_point(i64(ix), i64(iy), i64(iz));
      }
 
      template<typename Coord_, int dim_>
      bool check_point_nearest(const Tiny::Vector<Coord_, dim_>& point) const
      {
        // convert point to internal units
        std::array<i64, dim_> ipt;
        for(int i(0); i < dim_; ++i)
          ipt[std::size_t(i)] = i64(point[i] * Coord_(unit_size));
        return this->_check_point_nearest(ipt);
      }
 
      template<typename Coord_, int dim_>
      Real sample_point(const Tiny::Vector<Coord_, dim_>& point) const
      {
        // convert point to internal units
        std::array<i64, dim_> ipt;
        for(int i(0); i < dim_; ++i)
          ipt[std::size_t(i)] = i64(point[i] * Coord_(unit_size));
        return this->_sample_point(ipt);
      }
 
      template<typename Coord_, int dim_>
      bool check_box(const Tiny::Vector<Coord_, dim_>& box_min, const Tiny::Vector<Coord_, dim_>& box_max) const
      {
        // convert coordinates to map indices
        std::array<u64, dim_> ibox_min, ibox_max;
        for(std::size_t i(0); i < std::size_t(dim_); ++i)
        {
          // make sure the map indices do not exceed the bounding box of the voxel map
          i64 imin = _map_coord_idx_lower(i64(box_min[int(i)] * Coord_(unit_size)), i);
          if(imin < i64(0))
          {
            if(this->_out_of_bounds_value)
              return true;
            ibox_min[i] = u64(0);
          }
          else
            ibox_min[i] = u64(imin);
 
          i64 imax = _map_coord_idx_upper(i64(box_max[int(i)] * Coord_(unit_size)), i);
          if(i64(this->_num_points[i]) < imax + 1)
          {
            if(this->_out_of_bounds_value)
              return true;
            ibox_max[i] = i64(this->_num_points[i]) - 1;
          }
          else
            ibox_max[i] = u64(imax);
        }
 
        return this->_check_box(ibox_min, ibox_max);
      }
 
      template<typename Coord_, int dim_>
      Real sample_box(const Tiny::Vector<Coord_, dim_>& box_min, const Tiny::Vector<Coord_, dim_>& box_max) const
      {
        std::array<i64, dim_> ibox_min, ibox_max;
        for(std::size_t i(0); i < std::size_t(dim_); ++i)
        {
          ibox_min[i] = i64(box_min[int(i)] * Coord_(unit_size));
          ibox_max[i] = i64(box_max[int(i)] * Coord_(unit_size));
        }
        return this->_sample_box(ibox_min, ibox_max);
      }
 
      void export_to_bmp(const String& filename_prefix) const;
 
      void export_plane_to_bmp(const String& filename, Real z_coord) const;
 
      void export_plane_to_bmp(std::ostream& os, Real z_coord) const;
 
      void render_to_bmp(const String& filename_prefix, Index width, Index height, Index depth) const;
 
      void render_plane_to_bmp(const String& filename, Index width, Index height, Real z_min, Real z_max) const;
 
      void render_plane_to_bmp(std::ostream& os, Index width, Index height, Real z_min, Real z_max) const;
 
    protected:
      u64 _map_coord_idx_nearest(i64 xyz, std::size_t sdim) const;
 
      i64 _map_coord_idx_lower(i64 xyz, std::size_t sdim) const;
 
      i64 _map_coord_idx_upper(i64 xyz, std::size_t sdim) const;
 
      Real _calc_sample_rate(i64 xyz, i64 idx, std::size_t sdim) const;
 
      bool _check_point(i64 xidx, i64 yidx, i64 zidx) const;
 
      bool _check_point_nearest(const std::array<i64, 1>& p) const;
 
      bool _check_point_nearest(const std::array<i64, 2>& p) const;
 
      bool _check_point_nearest(const std::array<i64, 3>& p) const;
 
      bool _check_box(const std::array<u64, 1>& box_min, const std::array<u64, 1>& box_max) const;
 
      bool _check_box(const std::array<u64, 2>& box_min, const std::array<u64, 2>& box_max) const;
 
      bool _check_box(const std::array<u64, 3>& box_min, const std::array<u64, 3>& box_max) const;
 
      Real _sample_point_1d_x(i64 px, i64 xidx, i64 yidx = i64(0), i64 zidx = i64(0)) const;
 
      Real _sample_point_1d_y(i64 py, i64 xidx, i64 yidx, i64 zidx = i64(0)) const;
 
      Real _sample_point_1d_z(i64 pz, i64 xidx, i64 yidx, i64 zidx) const;
 
      Real _sample_point_2d(i64 px, i64 py, i64 xidx, i64 yidx, i64 zidx = i64(0)) const;
 
      Real _sample_point_3d(i64 px, i64 py, i64 pz, i64 xidx, i64 yidx, i64 zidx) const;
 
      Real _sample_point(const std::array<i64, 1>& p) const;
 
      Real _sample_point(const std::array<i64, 2>& p) const;
 
      Real _sample_point(const std::array<i64, 3>& p) const;
 
      Real _sample_box_1d(i64 px0, i64 px1, i64 xidx0, i64 xidx1, i64 yidx = i64(0), i64 zidx = i64(0)) const;
 
      Real _sample_box_2d(i64 px0, i64 px1, i64 py0, i64 py1, i64 xidx0, i64 xidx1, i64 yidx0, i64 yidx1, i64 zidx = i64(0)) const;
 
      Real _sample_box_3d(i64 px0, i64 px1, i64 py0, i64 py1, i64 pz0, i64 pz1, i64 xidx0, i64 xidx1, i64 yidx0, i64 yidx1, i64 zidx0, i64 zidx1) const;
 
      Real _sample_box(const std::array<i64, 1>& box_min, const std::array<i64, 1>& box_max) const;
 
      Real _sample_box(const std::array<i64, 2>& box_min, const std::array<i64, 2>& box_max) const;
 
      Real _sample_box(const std::array<i64, 3>& box_min, const std::array<i64, 3>& box_max) const;
 
      u64 _compute_domain_coverage() const;
 
      void _compress_voxel_map_line(const std::vector<int>& mask, const u64 line);
 
      template<typename CoordType_>
      void _compute_voxel_map_lines(VoxelMasker<CoordType_, 1>& masker, u64 DOXY(beg), u64 DOXY(end), u64 DOXY(offset))
      {
        XASSERTM(_num_points[1] == 1ull, "cannot use a 1D VoxelMasker to create a 2D/3D voxel mask");
        XASSERTM(_num_points[2] == 1ull, "cannot use a 1D VoxelMasker to create a 2D/3D voxel mask");
        const CoordType_ x_min = CoordType_(_bbox_min[0]) / CoordType_(unit_size);
        const CoordType_ x_max = CoordType_(_bbox_max[0]) / CoordType_(unit_size);
        std::vector<int> line_mask(this->_num_points[0], 0);
        Tiny::Vector<CoordType_, 1> coords;
        masker.mask_line(line_mask, x_min, x_max, coords);
        this->_compress_voxel_map_line(line_mask, 0);
      }
 
      template<typename CoordType_>
      void _compute_voxel_map_lines(VoxelMasker<CoordType_, 2>& masker, u64 beg, u64 end, u64 offset)
      {
        XASSERTM(_num_points[2] == 1ull, "cannot use a 2D VoxelMasker to create a 3D voxel mask");
        const CoordType_ x_min = CoordType_(_bbox_min[0]) / CoordType_(unit_size);
        const CoordType_ x_max = CoordType_(_bbox_max[0]) / CoordType_(unit_size);
        const CoordType_ y_min = CoordType_(_bbox_min[1]) / CoordType_(unit_size);
        const CoordType_ y_max = CoordType_(_bbox_max[1]) / CoordType_(unit_size);
        std::vector<int> line_mask(this->_num_points[0], 0);
        Tiny::Vector<CoordType_, 2> coords;
 
        for(u64 i = beg; i < end; ++i)
        {
          coords[1] = y_min + (y_max - y_min) * CoordType_(i % this->_num_points[1]) / CoordType_(this->_num_points[1] - 1u);
          masker.mask_line(line_mask, x_min, x_max, coords);
          this->_compress_voxel_map_line(line_mask, i - offset);
        }
      }
 
      template<typename CoordType_>
      void _compute_voxel_map_lines(VoxelMasker<CoordType_, 3>& masker, u64 beg, u64 end, u64 offset)
      {
        const CoordType_ x_min = CoordType_(_bbox_min[0]) / CoordType_(unit_size);
        const CoordType_ x_max = CoordType_(_bbox_max[0]) / CoordType_(unit_size);
        const CoordType_ y_min = CoordType_(_bbox_min[1]) / CoordType_(unit_size);
        const CoordType_ y_max = CoordType_(_bbox_max[1]) / CoordType_(unit_size);
        const CoordType_ z_min = CoordType_(_bbox_min[2]) / CoordType_(unit_size);
        const CoordType_ z_max = CoordType_(_bbox_max[2]) / CoordType_(unit_size);
 
        FEAT_PRAGMA_OMP(parallel)
        {
          std::vector<int> line_mask(this->_num_points[0], 0);
          Tiny::Vector<CoordType_, 3> coords;
 
          FEAT_PRAGMA_OMP(for schedule(dynamic, 16))
          for(u64 i = beg; i < end; ++i)
          {
            // line = iz * this->_num_points[1] + iy
            coords[1] = y_min + (y_max - y_min) * CoordType_(i % this->_num_points[1]) / CoordType_(this->_num_points[1] - 1u);
            coords[2] = z_min + (z_max - z_min) * CoordType_(i / this->_num_points[1]) / CoordType_(this->_num_points[2] - 1u);
            masker.mask_line(line_mask, x_min, x_max, coords);
            this->_compress_voxel_map_line(line_mask, i - offset);
          }
        }
      }
 
      template<typename CoordType, int dim_>
      void _compute_voxel_map(const Dist::Comm& comm, VoxelMasker<CoordType, dim_>& masker, bool gather_to_all)
      {
        // single process or MPI parallel?
        if(comm.size() <= 1)
        {
          // compute all voxel map lines
          this->_compute_voxel_map_lines(masker, 0u, this->_num_lines, 0u);
        }
        else if(u64(comm.size()) <= this->_num_lines)
        {
          // there are at least as many lines in the voxel map as there are MPI processes; this is the usual case
          // we can use the entire communicator (usually MPI_COMM_WORLD) for the voxel map computation
          u64 comm_rank = u64(comm.rank());
          u64 comm_size = u64(comm.size());
 
          // compute line count for a single process and round up; this may be bigger than _num_lines,
          // so the last process may have less lines to chew through than all other processes
          u64 line_count = (this->_num_lines + comm_size - 1u) / comm_size;
          std::size_t block_len = this->_stride_line * line_count;
 
          // allocate voxel map and format to zero
          std::size_t vmap_size = block_len;
          // do we need the entire map on this process?
          if(gather_to_all || (comm_rank == 0))
            vmap_size *= u64(comm_size);
          this->_voxel_map.resize(vmap_size, 0u);
 
          // compute first and last+1 line of this process; make sure the last process doesn't go beyond the total count
          u64 line_beg = comm_rank * line_count;
          u64 line_end = Math::min((comm_rank + 1u) * line_count, this->_num_lines);
 
          // debug output
          //this->_comm.allprint(String(">>>") + stringify(line_beg).pad_front(6) + ":" + stringify(line_count) + ">" + stringify(line_end).pad_front(6));
 
          // compute lines for this process
          this->_compute_voxel_map_lines(masker, line_beg, line_end, gather_to_all ? u64(0) : line_beg);
 
          // perform an in-place allgather to synchronize the voxel map over all processes
          if(gather_to_all)
            comm.allgather(this->_voxel_map.data(), block_len, this->_voxel_map.data(), block_len);
          else
            comm.gather(this->_voxel_map.data(), block_len, this->_voxel_map.data(), block_len, 0);
        }
        else if(this->_num_lines < u64(comm.size()))
        {
          // more processes than slag lines; this is a borderline scenario, but we have to support it anyways
          // create a sub-communicator and then compute mask on that sub-communicator
          Dist::Comm sub_comm = comm.comm_create_range_incl(int(this->_num_lines));
 
          // does this process participate in the sub-communicator?
          if(!sub_comm.is_null())
          {
            // get sub-communicator rank and size
            u64 comm_rank = u64(sub_comm.rank());
            u64 comm_size = u64(sub_comm.size());
 
            // compute line count for a single process and round up; this may be bigger than _num_lines,
            // so the last process may have less lines to chew through than all other processes
            u64 line_count = (this->_num_lines + comm_size - 1u) / comm_size;
            std::size_t block_len = this->_stride_line * line_count;
 
            // allocate voxel map and format to zero
            std::size_t vmap_size = block_len;
            // do we need the entire map on this process?
            if(gather_to_all || (comm_rank == 0))
              vmap_size *= u64(comm_size);
            this->_voxel_map.resize(vmap_size, 0u);
 
            // compute first and last+1 line of this process; make sure the last process doesn't go beyond the total count
            u64 line_beg = comm_rank * line_count;
            u64 line_end = Math::min((comm_rank + 1u) * line_count, this->_num_lines);
 
            // debug output
            //this->_comm.allprint(String(">>>") + stringify(line_beg).pad_front(6) + ":" + stringify(line_count) + ">" + stringify(line_end).pad_front(6));
 
            // compute lines for this process
            this->_compute_voxel_map_lines(masker, line_beg, line_end, gather_to_all ? u64(0) : line_beg);
 
            // perform an in-place gather to synchronize the voxel map on rank 0
            sub_comm.gather(this->_voxel_map.data(), block_len, this->_voxel_map.data(), block_len, 0);
          }
          else
          {
            // just allocate the voxel map vector for the following broadcast
            if(gather_to_all)
              this->_voxel_map.resize(this->_stride_line * this->_num_lines, 0u);
          }
 
          // broadcast from rank 0 to all processes our entire domain communicator
          if(gather_to_all)
            comm.bcast(this->_voxel_map.data(), this->_stride_line * this->_num_lines, 0);
        }
 
        // compute domain coverage on rank 0
        if(comm.rank() == 0)
          this->_coverage = this->_compute_domain_coverage();
 
        // broadcast domain coverage to all ranks
        comm.bcast(&this->_coverage, std::size_t(1u), 0);
 
        // that's it
      }
 
      void _export_plane_to_bmp(std::ostream& os, u64 plane) const;
 
      void _render_plane_to_bmp(std::ostream& os, Index width, Index height, i64 box_min_z, i64 box_max_z) const;
    }; // class VoxelMap
  } // namespace Geometry
} // namespace FEAT