DensityMap.h

Go to the documentation of this file.

/**
 *  \file IMP/em/DensityMap.h
 *  \brief Class for handling density maps.
 *
 *  Copyright 2007-2022 IMP Inventors. All rights reserved.
 *
 */

#ifndef IMPEM_DENSITY_MAP_H
#define IMPEM_DENSITY_MAP_H
#include <IMP/Pointer.h>
#include <IMP/em/em_config.h>
#include "DensityHeader.h"
#include "MapReaderWriter.h"
#include <IMP/Object.h>
#include <IMP/algebra/Vector3D.h>
#include <IMP/algebra/BoundingBoxD.h>
#include <IMP/algebra/Transformation3D.h>
#include <IMP/Object.h>
#include <boost/scoped_array.hpp>
#include <iostream>
#include <iomanip>
#include <IMP/algebra/standard_grids.h>
#include <IMP/base_types.h>
//#include <IMP/statistics/Histogram.h>

IMPEM_BEGIN_NAMESPACE

class DensityMap;

//! Create a copy of another map
IMPEMEXPORT DensityMap *create_density_map(const DensityMap *other);

//! Create an empty density map from a bounding box
IMPEMEXPORT DensityMap *create_density_map(const algebra::BoundingBox3D &bb,
                                           double spacing);

//! Create an empty density map
IMPEMEXPORT DensityMap *create_density_map(int nx, int ny, int nz,
                                           double spacing);

//! Read a density map from a file (using the given reader) and return it.
/** \relates DensityMap
 */
IMPEMEXPORT DensityMap *read_map(std::string filename, MapReaderWriter *reader);

//! Read a density map from a file and return it.
/** Guess the file type from the file name. The file formats supported are:
    - .mrc/.map
    - .em
    - .vol
    - .xplor
    \relates DensityMap
*/
IMPEMEXPORT DensityMap *read_map(std::string filename);

//! Write a density map to a file using the given writer.
/** \relates DensityMap
*/
IMPEMEXPORT void write_map(DensityMap *m, std::string filename,
                           MapReaderWriter *writer);

//! Write a density map to a file.
/** Guess the file type from the
    file name. The file formats supported are:
    - .mrc/.map
    - .em
    - .vol
    - .xplor
    \relates DensityMap
*/
IMPEMEXPORT void write_map(DensityMap *m, std::string filename);

//! Get the bounding box for a map.
/** \param[in] m a density map
    \param[in] threshold find the bounding box for voxels
                         with value above the threshold
 */
IMPEMEXPORT algebra::BoundingBoxD<3> get_bounding_box(const DensityMap *m,
                                                      Float threshold);
//! Estimate the molecular mass from a map.
/** \param[in] m a density map
    \param[in] threshold consider volume of only voxels above this threshold
    \note The method assumes 1.21 cubic Å per Dalton (Harpaz 1994).
 */
IMPEMEXPORT Float approximate_molecular_mass(DensityMap *m, Float threshold);

//! Class for handling density maps.
/** \note The location of a voxel is its center. That is important
     for sampling function as well as for functions
     like get_location_in_dim_by_voxel.
 */
class IMPEMEXPORT DensityMap : public IMP::Object {
  typedef IMP::algebra::DenseGrid3D<double> DGrid;
  IMP_NO_SWIG(friend IMPEMEXPORT DensityMap *read_map(std::string filename,
                                                      MapReaderWriter *reader));
  IMP_NO_SWIG(friend IMPEMEXPORT void write_map(DensityMap *m,
                                                std::string filename,
                                                MapReaderWriter *writer));

 public:
  DensityMap(std::string name = "DensityMap%1%");
  //! Construct a density map as instructed in the input header
  DensityMap(const DensityHeader &header, std::string name = "DensityMap%1%");

  //! Set the density voxels to some value and reset the management flags.
  /** \param[in] value all of the density voxels will have this value
   */
  void reset_data(float value = 0.0);

  //! Calculates RMSD and mean of a map values as stored in the header.
  /** The header stores whether the map is normalized.
   */
  emreal calcRMS();

  //! Normalize the density voxels according to standard deviation (stdv).
  /** The mean is subtracted from the map, which is then divided by the stdv.
      The normalization flag is set to avoid repeated computation. */
  void std_normalize();

  inline bool is_normalized() const { return normalized_; }

  //! Calculate the location of a given voxel in a given dimension
  /** \param[in] index The voxel index
      \param[in] dim   The dimension of interest (x:=0,y:=1,z:=2)
      \return the location (x,y,z) (in angstroms) of a given voxel. If the
              index is not part of the map, the function returns -1.
   */
  float get_location_in_dim_by_voxel(long index, int dim) const;

  //! Calculate the voxel of a given xyz index
  /** \param[in] x The voxel index on the x axis of the grid
      \param[in] y The voxel index on the y axis of the grid
      \param[in] z The voxel index on the z axis of the grid
      \return the voxel index.
   */
  inline long xyz_ind2voxel(int x, int y, int z) const {
    return z * header_.get_nx() * header_.get_ny() + y * header_.get_nx() + x;
  }

  //! Calculate the voxel of a given location
  /** \param[in] x The position (in angstroms) of the x coordinate
      \param[in] y The position (in angstroms) of the y coordinate
      \param[in] z The position (in angstroms) of the z coordinate
      \return the voxel index of a given position. If the position is out of
              the boundaries of the map, the function returns -1.
   */
  long get_voxel_by_location(float x, float y, float z) const;

  //! Calculate the voxel of a given location
  /** \param[in] v The position (in angstroms)
      \return the voxel index of a given position. If the position is out of
              the boundaries of the map, the function returns -1.
   */
  long get_voxel_by_location(const algebra::Vector3D &v) const {
    return get_voxel_by_location(v[0], v[1], v[2]);
  }
  //! Calculate dimension index of a given location
  /** \param[in] v The position (in angstroms)
      \param[in] ind dimension index (X:0,Y:1 or Z:2)
   */
  int get_dim_index_by_location(const algebra::Vector3D &v, int ind) const;

  //! Calculate the location of a given voxel.
  /** \param[in] index The voxel index
      \return the location (x,y,z) (in angstroms) of a given voxel.
  */
  algebra::Vector3D get_location_by_voxel(long index) const {
    IMP_USAGE_CHECK(index >= 0 && index < get_number_of_voxels(),
                    "invalid map index");
    IMP_USAGE_CHECK(
        loc_calculated_,
        "locations should be calculated prior to calling this function");
    return algebra::Vector3D(x_loc_[index], y_loc_[index], z_loc_[index]);
  }

  bool is_xyz_ind_part_of_volume(int ix, int iy, int iz) const;

  //! Checks whether a given point is in the grid
  /** \param[in] x The position (in angstroms) of the x coordinate
      \param[in] y The position (in angstroms) of the y coordinate
      \param[in] z The position (in angstroms) of the z coordinate
      \return true if the point is part of the grid, false otherwise.
   */
  bool is_part_of_volume(float x, float y, float z) const;
  //! Checks whether a given point is in the grid
  /** \param[in] v The position (in angstroms)
      \return true if the point is part of the grid, false otherwise.
   */
  bool is_part_of_volume(const algebra::Vector3D &v) const {
    return is_part_of_volume(v[0], v[1], v[2]);
  }

  //! Gets the value of the voxel located at (x,y,z)
  /** \param[in] x The position (in angstroms) of the x coordinate
      \param[in] y The position (in angstroms) of the y coordinate
      \param[in] z The position (in angstroms) of the z coordinate
      \return the value of the voxel located at (x,y,z)
      \exception IndexException The point is not covered by the grid.
      \note the value is not interpolated between this and neighboring
            voxels. For that, see get_density().
   */
  emreal get_value(float x, float y, float z) const;
  emreal get_value(const algebra::Vector3D &point) const {
    return get_value(point[0], point[1], point[2]);
  }

  //! Gets the value of the voxel at a given index
  /** \param[in] index voxel number in physical sense, NOT logical
   */
  emreal get_value(long index) const;

  //! Set the value of the voxel at a given index
  /** \param[in] index voxel number in physical sense, NOT logical
      \param[in] value value
   */
  void set_value(long index, emreal value);

  //! Set the value of the voxel at given coordinates
  void set_value(float x, float y, float z, emreal value);

  //! Sets the origin of the header
  /** \param x the new x (angstroms)
      \param y the new y (angstroms)
      \param z the new z (angstroms)
   */
  void set_origin(float x, float y, float z);
  void set_origin(const IMP::algebra::Vector3D &v) {
    set_origin(v[0], v[1], v[2]);
  }

  algebra::Vector3D get_origin() const {
    return algebra::Vector3D(header_.get_origin(0), header_.get_origin(1),
                             header_.get_origin(2));
  }

  algebra::Vector3D get_top() const {
    return algebra::Vector3D(header_.get_top(0), header_.get_top(1),
                             header_.get_top(2));
  }

  //! Returns a read-only pointer to the header of the map
  const DensityHeader *get_header() const { return &header_; }

  //! Returns a pointer to the header of the map in a writable version
  DensityHeader *get_header_writable() { return &header_; }

#if !defined(IMP_DOXYGEN) && !defined(SWIG)
  //! Returns the x-location of the map
  /**
  \exception InvalidStateException The locations have not been calculated.
  */
  float *get_x_loc() const {
    IMP_USAGE_CHECK(loc_calculated_,
                    "x location requested before being calculated");
    return x_loc_.get();
  }
  //! Returns the y-location of the map
  /**
  \exception InvalidStateException The locations have not been calculated.
  */
  float *get_y_loc() const {
    IMP_USAGE_CHECK(loc_calculated_,
                    "y location requested before being calculated");
    return y_loc_.get();
  }
  //! Returns the z-location of the map
  /**
  \exception InvalidStateException The locations have not been calculated.
  */
  float *get_z_loc() const {
    IMP_USAGE_CHECK(loc_calculated_,
                    "z location requested before being calculated");
    return z_loc_.get();
  }

  emreal *get_data() const { return data_.get(); }
#endif

  //! Checks if two maps have the same origin
  /** \param[in] other the map to compare with
      \return true if the two maps have the same origin
   */
  bool same_origin(const DensityMap *other) const;

  //! Checks if two maps have the same dimensions
  /** \param[in] other the map to compare with
      \return true if the two maps have the same dimensions
   */
  bool same_dimensions(const DensityMap *other) const;

  //! Checks if two maps have the same voxel size
  /** \param[in] other the map to compare with
      \return true if the two maps have the same voxel size
   */
  bool same_voxel_size(const DensityMap *other) const;
  //! Get the centroid of all voxels with density above a given threshold
  /** \param[in] threshold the input threshold
  */
  algebra::Vector3D get_centroid(emreal threshold = 0.0) const;
  //! Returns the value of the voxel with the highest density.
  emreal get_max_value() const;
  //! Returns the value of the voxel with the lowest density.
  emreal get_min_value() const;

  //! Sums two grids; the result is kept in the map.
  /** \param[in] other the other map
      \note The shared extend is summed
      \note The two maps should have the same voxelsize and other
            should be contained within this map
   */
  void add(const DensityMap *other);

  //! Add a number to every voxel in the map
  /** \param[in] d Value to add

   */
  void add(Float d);

  //! Pick the max value between two corresponding voxels between two maps
  /** The result is kept in the map.
      \param[in] other the other map
      \note The two maps should have the same voxelsize and the same dimensions
   */
  void pick_max(const DensityMap *other);

  //! Get the number of map voxels
  long get_number_of_voxels() const {
    return header_.get_number_of_voxels();
  }

  //! Set the map dimension and reset all voxels to 0
  /** \param[in] nx x-dimension (voxels)
      \param[in] ny y-dimension (voxels)
      \param[in] nz z-dimension (voxels)
      \note the origin and spacing remain unchanged
   */
  void set_void_map(int nx, int ny, int nz);

  //! Increase the dimension of the map
  /** The function pads zeroes to the right-upper section on the map.
      The original content of the map will be in the lower XYZ part of the map.
      \param[in] nx the number of voxels on the X axis
      \param[in] ny the number of voxels on the Y axis
      \param[in] nz the number of voxels on the Z axis
      \param[in] val all additional voxels will have this value
      \exception UsageException if the input numbers of x/y/z voxels are
                 smaller than those currently in the map
   */
  void pad(int nx, int ny, int nz, float val = 0.0);

  //! Create a new padded map
  /** Given this map of size [nx,ny,nz],
      the new map is of size [2*mrg_x+nx,2*mrg_y+ny,2*mrg_z+nz].
      The new map will consist of the values of the old map,
      with a padded margin on all sides.

      \param[in] mrg_x number of margin voxels to add on both right and
                       left on the X axis
      \param[in] mrg_y number of margin voxels to add on both right and
                       left on the Y axis
      \param[in] mrg_z number of margin voxels to add on both right and
                       left on the Z axis
      \param[in] val   ignored
      \return the new padded map.
      \exception UsageException if the input numbers of x/y/z voxels are
                 smaller than those currently in the map
   */
  DensityMap *pad_margin(int mrg_x, int mrg_y, int mrg_z, float val = 0.0);

  //! Create and return a new cropped map
  /** The margins are determined to be the bounding box
      with density values below the input
      \param[in] threshold used for cropping
      \return the new cropped map.
   */
  DensityMap *get_cropped(float threshold);
  
  //! Create and return a new cropped map based on particles
  /** The map is cropped based on voxel proximity to a set of particles.
      All voxel centers farther than [distance] away from any
      particle center in [ps] will have their density value set to 0.0.
      
      setting inverse to true will set all voxel centers within
      [distance] of a particle to 0.0.

      \param[in] ps List of particles used to crop map
      \param[in] distance Distance from particles at which to crop map
      \param[in] inverse Set to true to crop the particle volume of the map
      \return the new cropped map.
   */
  DensityMap *get_cropped(Particles ps, double distance, bool inverse = false, bool keep_map_dimensions = false);

  //! Create and return a new cropped map with the given bounding box extent
  /** \param[in] bb the bounding box
      \return the new cropped map.
      \note If the input bounding box is larger than the density box,
            it is snapped to the right size.
   */
  DensityMap *get_cropped(const algebra::BoundingBox3D &bb);

  //! Get the maximum value in a XY plane indicated by a Z index
  float get_maximum_value_in_xy_plane(int z_ind);
  //! Get the maximum value in a XZ plane indicated by a Y index
  float get_maximum_value_in_xz_plane(int y_ind);
  //! Get the maximum value in a YZ plane indicated by a X index
  float get_maximum_value_in_yz_plane(int x_ind);

  //! Multiply in place each voxel in the map by the input factor
  /** The result is kept in the map.
      \param[in] factor the multiplication factor
   */
  void multiply(float factor);

  //! Print the locations of all voxels with value above a given threshold
  std::string get_locations_string(float t);

  //! Update the voxel size of the map
  void update_voxel_size(float new_apix);

  //! Get the voxel size of the map
  /** \note Use update_voxel_size() to set the spacing value.
  */
  Float get_spacing() const { return header_.get_spacing(); }
  //! Calculates the coordinates that correspond to all voxels.
  void calc_all_voxel2loc();
  //! Copy map into this map
  void copy_map(const DensityMap *other);
  IMP_OBJECT_METHODS(DensityMap);
#if !defined(IMP_DOXYGEN) && !defined(SWIG)
  //! Convolute a kernel with a map and write results into current map
  /**
     \param[in] other the map to convolute
     \param[in] kernel an array of kernel values. The data is in ZYX
     order, Z is the slowest.
     \param[in] dim_len the kernel array
   */
  void convolute_kernel(DensityMap *other, double *kernel, int dim_len);
  //! Convolute a kernel with this map and write results to this map
  /**
     \param[in] kernel an array of kernel values. The data is in ZYX
     order, Z is the slowest.
     \param[in] dim_len the kernel array
   */
  void convolute_kernel(double *kernel, int dim_len) {
    Pointer<DensityMap> cmap = em::create_density_map(this);
    cmap->set_was_used(true);
    convolute_kernel(cmap, kernel, dim_len);
    cmap = static_cast<DensityMap *>(nullptr);
  }
#endif
  int lower_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
  int upper_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
  inline bool get_rms_calculated() const { return rms_calculated_; }
  int get_dim_index_by_location(float loc_val, int ind) const;

 protected:
  //!update the header values  -- still in work
  void update_header();
  void reset_all_voxel2loc();

  void allocated_data();
  void float2real(float *f_data, boost::scoped_array<emreal> &r_data);
  void real2float(emreal *r_data, boost::scoped_array<float> &f_data);

  DensityHeader header_;              // holds all the info about the map
  boost::scoped_array<emreal> data_;  // the order is ZYX (Z-slowest)
  bool data_allocated_;

  //! Locations (centers) for each of the voxels of the map (they are
  //! precomputed and each one is of size nvox, where nvox is the
  //! size of the map)
  boost::scoped_array<float> x_loc_, y_loc_, z_loc_;
  //! true if the locations have already been computed
  bool loc_calculated_;

  bool normalized_;
  bool rms_calculated_;
};

inline algebra::BoundingBoxD<3> get_bounding_box(const DensityMap *m) {
  const DensityHeader *h = m->get_header();
  float hspace = m->get_spacing() / 2.0;
  algebra::Vector3D lb = m->get_origin()
                         - algebra::Vector3D(hspace, hspace, hspace);
  return algebra::BoundingBoxD<3>(
      lb,
      lb + algebra::Vector3D(m->get_spacing() * h->get_nx(),
                             m->get_spacing() * h->get_ny(),
                             m->get_spacing() * h->get_nz()));
}

#if !defined(IMP_DOXYGEN) && !defined(SWIG)
//! Calculate a bounding box around a 3D point within the EM grid
/**
\param[in] d_map the density map
\param[in] x coordinates of the point to sample around
\param[in] y coordinates of the point to sample around
\param[in] z coordinates of the point to sample around
\param[in] kdist the length of the box
\param[out] iminx the minimum index on the X axis of the output bounding box
\param[out] iminy the minimum index on the Y axis of the output bounding box
\param[out] iminz the minimum index on the Z axis of the output bounding box
\param[out] imaxx the maximum index on the X axis of the output bounding box
\param[out] imaxy the maximum index on the Y axis of the output bounding box
\param[out] imaxz the maximum index on the Z axis of the output bounding box
 */
inline void calc_local_bounding_box(const em::DensityMap *d_map, double x,
                                    double y, double z, float kdist, int &iminx,
                                    int &iminy, int &iminz, int &imaxx,
                                    int &imaxy, int &imaxz) {
  const DensityHeader *h = d_map->get_header();
  iminx = d_map->lower_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
  iminy = d_map->lower_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
  iminz = d_map->lower_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
  imaxx = d_map->upper_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
  imaxy = d_map->upper_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
  imaxz = d_map->upper_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
}
#endif

IMP_OBJECTS(DensityMap, DensityMaps);

//! Get density value at point v, interpolating linearly from the sample values.
/** The resulting function is C0 over R3.
    \relates DensityMap
*/
IMPEMEXPORT double get_density(const DensityMap *m, const algebra::Vector3D &v);

//! Return a new density map containing a rotated version of the old one.
/** Only voxels whose value is above threshold are considered when
    computing the bounding box of the new map (set IMP::em::get_bounding_box()).
    \relates DensityMap
*/
IMPEMEXPORT DensityMap *get_transformed(const DensityMap *input,
                                        const algebra::Transformation3D &tr,
                                        double threshold);

//! Return a new density map containing a rotated version of the old one.
/** The dimension of the new map is the same as the old one.
    \relates DensityMap
*/
IMPEMEXPORT DensityMap *get_transformed(DensityMap *input,
                                        const algebra::Transformation3D &tr);

//! Get a resampled version of the map.
/** The spacing is multiplied by scaling, i.e. scaling values greater
    than 1 increase the voxel size.
    \relates DensityMap
*/
IMPEMEXPORT DensityMap *get_resampled(DensityMap *input, double scaling);

//! Rotate a density map into another map
/** \param[in] source the map to transform
    \param[in] tr transform the from density map by this transformation
    \param[out] into the map to transform into
    \param[in] calc_rms if true RMS is calculated on the transformed map
    \relates DensityMap
*/
IMPEMEXPORT void get_transformed_into(const DensityMap *source,
                                      const algebra::Transformation3D &tr,
                                      DensityMap *into, bool calc_rms = true);
IMPEMEXPORT void get_transformed_into2(const DensityMap *source,
                                       const algebra::Transformation3D &tr,
                                       DensityMap *into);

inline bool get_interiors_intersect(const DensityMap *d1,
                                    const DensityMap *d2) {
  return get_interiors_intersect(get_bounding_box(d1), get_bounding_box(d2));
}
#if 0
//! Get a histogram of density values
/**
\param[in] dmap the density map to analyze
\param[in] threshold only add voxels with value above this threshold
                     to the histogram
\param[in] num_bins the number of bins to have in the histogram
 \see DensityMap
*/
// IMPEMEXPORT statistics::Histogram
// get_density_histogram(const DensityMap *dmap, float threshold,int num_bins);
#endif

//! Get a segment of the map according to xyz indexes
/** \note the output map will cover the region
    [[nx_start,nx_end],[]ny_start,ny_end,[nz_start,nz_end]]
 */
IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment, int nx_start,
                                    int nx_end, int ny_start, int ny_end,
                                    int nz_start, int nz_end);

//! Get a segment of the map covered by the input points
IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment,
                                    algebra::Vector3Ds vecs, float dist);
//! Get a segment of the map covered by another map
IMPEMEXPORT DensityMap *get_segment_by_masking(DensityMap *map_to_segment,
                                               DensityMap *mask,
                                               float mas_threshold);

//! Return a map with 0 for all voxels below the threshold and 1 for those above
/** \param[in] orig_map the map to binarize
    \param[in] threshold values below the threshold are set to 0 and 1 otherwise
    \param[in] reverse if true values above the threshold
                       are set to 0 and 1 otherwise
 */
IMPEMEXPORT DensityMap *binarize(DensityMap *orig_map, float threshold,
                                 bool reverse = false);

//! Return a map with 0 for all voxels below the threshold
/** \param[in] orig_map the map to threshold
    \param[in] threshold values below the threshold are set to 0 and
                         left unchanged otherwise
 */
IMPEMEXPORT DensityMap *get_threshold_map(const DensityMap *orig_map,
                                          float threshold);

//! Return a density map for which voxel i contains the result of m1[i]*m2[i].
/** The function assumes m1 and m2 are of the same dimensions.
 */
IMPEMEXPORT DensityMap *multiply(const DensityMap *m1, const DensityMap *m2);

//! Return a convolution between density maps m1 and m2.
/** The function assumes m1 and m2 are of the same dimensions.
 */
IMPEMEXPORT double convolute(const DensityMap *m1, const DensityMap *m2);

//! Return the sum of all voxels
IMPEMEXPORT double get_sum(const DensityMap *m1);

//! Return a map where each voxel is the maximum value from the input maps.
/** \note all input maps should have the same extent
 */
IMPEMEXPORT DensityMap *get_max_map(DensityMaps maps);

//! Return a new map with an updated spacing
/** Input data is interpolated accordingly.
 */
IMPEMEXPORT
DensityMap *interpolate_map(DensityMap *in_map, double new_spacing);

//! Return a dense grid containing the voxels of the passed density map
/** The returned grid has the same bounding box as the map.
 */
IMPEMEXPORT
algebra::GridD<3, algebra::DenseGridStorageD<3, float>, float> get_grid(
    DensityMap *in_map);

//! Create a density map from an arbitrary IMP::algebra::GridD
template <class S, class V, class E>
inline DensityMap *create_density_map(
    const IMP::algebra::GridD<3, S, V, E> &arg) {
  IMP_FUNCTION_LOG;
  typedef IMP::algebra::GridD<3, S, V, E> Grid;
  IMP_USAGE_CHECK(
      std::abs(arg.get_unit_cell()[0] - arg.get_unit_cell()[1]) < .01,
      "The passed grid does not seem to have cubic voxels");
  Pointer<DensityMap> ret = create_density_map(
      algebra::get_bounding_box(arg), arg.get_unit_cell()[0]);
  IMP_USAGE_CHECK(arg.get_number_of_voxels(0) ==
                      static_cast<unsigned int>(ret->get_header()->get_nx()),
                  "X voxels don't match");
  IMP_USAGE_CHECK(arg.get_number_of_voxels(1) ==
                      static_cast<unsigned int>(ret->get_header()->get_ny()),
                  "Y voxels don't match");
  IMP_USAGE_CHECK(arg.get_number_of_voxels(2) ==
                      static_cast<unsigned int>(ret->get_header()->get_nz()),
                  "Z voxels don't match");
  for(typename Grid::Index i : arg.get_all_indexes()) {
    long vi = ret->xyz_ind2voxel(i[0], i[1], i[2]);
    ret->set_value(vi, arg[vi]);
  }
  return ret.release();
}

//! Return a density map with the values taken from the grid.
/** \relates DensityMap
 */
IMPEMEXPORT
DensityMap *create_density_map(algebra::DenseGrid3D<float> &grid);

//! Return a density map with the values taken from the grid.
/** \relates DensityMap
 */
IMPEMEXPORT
DensityMap *create_density_map(algebra::DenseGrid3D<double> &grid);


//! Return a binarized map with 1 for voxels that are internal in the input map
/** \relates DensityMap
 */
IMPEMEXPORT
DensityMap *get_binarized_interior(DensityMap *dmap);

//! Rasterize the particles into an existing density map.
/** \relates DensityMap
 */
IMPEMEXPORT
void add_to_map(DensityMap *dm, const Particles &pis);  // defined in
                                                        // SampledDensityMap.cpp

IMPEM_END_NAMESPACE

#endif /* IMPEM_DENSITY_MAP_H */