IMP  2.0.0
The Integrative Modeling Platform
DensityMap.h
Go to the documentation of this file.
1 /**
2  * \file IMP/em/DensityMap.h
3  * \brief Class for handling density maps.
4  *
5  * Copyright 2007-2013 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPEM_DENSITY_MAP_H
10 #define IMPEM_DENSITY_MAP_H
11 #include <IMP/Pointer.h>
12 #include <IMP/em/em_config.h>
13 #include "DensityHeader.h"
14 #include "MapReaderWriter.h"
15 #include <IMP/Object.h>
16 #include <IMP/algebra/Vector3D.h>
19 #include <IMP/Object.h>
20 #include <boost/scoped_array.hpp>
21 #include <iostream>
22 #include <iomanip>
24 //#include <IMP/statistics/Histogram.h>
25 
26 IMPEM_BEGIN_NAMESPACE
27 
28 class DensityMap;
29 
30 //! create a copy of another map
31 IMPEMEXPORT DensityMap *create_density_map(const DensityMap *other);
32 
33 //! Create an empty density map from a boudning box
34 IMPEMEXPORT DensityMap *create_density_map(
35  const algebra::BoundingBox3D &bb,
36  double spacing);
37 //! Create an empty density map
38 IMPEMEXPORT DensityMap *create_density_map(
39  int nx,int ny,int nz,
40  double spacing);
41 
42 /** Read a density map from a file and return it.
43  \relatesalso DensityMap
44 */
45 IMPEMEXPORT DensityMap* read_map(std::string filename, MapReaderWriter *reader);
46 
47 /** Read a density map from a file and return it. Guess the file type from the
48  file name. The file formats supported are:
49  - .mrc
50  - .em
51  - .vol
52  - .xplor
53  \relatesalso DensityMap
54 */
55 IMPEMEXPORT DensityMap* read_map(std::string filename);
56 
57 
58 /** Write a density map to a file.
59  \relatesalso DensityMap
60 */
61 IMPEMEXPORT void write_map(DensityMap* m, std::string filename,
62  MapReaderWriter *writer);
63 
64 /** Write a density map to a file.
65  Guess the file type from the
66  file name. The file formats supported are:
67  - .mrc
68  - .em
69  - .vol
70  - .xplor
71  \relatesalso DensityMap
72 */
73 IMPEMEXPORT void write_map(DensityMap* m, std::string filename);
74 
75 
76 //!
77 /**
78 \param[in] m a density map
79 \param[in] threshold find the boudning box for voxels
80  with value above the threshold
81  */
82 IMPEMEXPORT algebra::BoundingBoxD<3>
83  get_bounding_box(const DensityMap* m, Float threshold);
84 //!
85 /**
86 \param[in] m a density map
87 \param[in] threshold consider volume of only voxels above this threshold
88 \note The method assumes 1.21 cubic A per dalton (Harpaz 1994).
89  */
90 IMPEMEXPORT Float approximate_molecular_mass(DensityMap* m, Float threshold);
91 
92 
93 //! Class for handling density maps.
94 /** \note The location of a voxel is its center. That is important
95  for sampling function as well as for functions
96  like get_location_in_dim_by_voxel.
97  */
98 class IMPEMEXPORT DensityMap: public IMP::base::Object
99 {
101  IMP_NO_SWIG(friend IMPEMEXPORT DensityMap* read_map(std::string filename,
102  MapReaderWriter *reader));
103  IMP_NO_SWIG(friend IMPEMEXPORT void write_map(DensityMap* m,
104  std::string filename,
105  MapReaderWriter *writer));
106 
107 public:
108  DensityMap(std::string name= "DensityMap%1%");
109  //! Construct a density map as intructed in the input header
110  DensityMap(const DensityHeader &header, std::string name= "DensityMap%1%");
111  void release(){
112  /* if (data_ != nullptr) {
113  data_.reset(nullptr);}
114  x_loc_.reset();
115  if (y_loc_ != nullptr){
116  y_loc_.reset(nullptr);}
117  if (z_loc_ != nullptr){
118  z_loc_.reset(nullptr);}*/
119  }
120 
121  //! Set the density voxels to some calue and reset the managment flags.
122  /**
123  \param[in] value all of the density voxels will have this value
124  */
125  void reset_data(float value=0.0);
126 
127  //! Calculates RMSD and mean of a map values are stored in the header.
128  /** The header stores whether map is normalized.
129  */
130  emreal calcRMS();
131 
132 
133  //! Normailze the density voxles according to standard deviation (stdv).
134  /** The mean is subtracted from the map, which is then divided by the stdv.
135  The normalization flag is set to avoid repeated computation */
136  void std_normalize();
137 
138  inline bool is_normalized() const {return normalized_;}
139 
140  //! Calculate the location of a given voxel in a given dimension
141  /** \param[in] index The voxel index
142  \param[in] dim The dimesion of intereset (x:=0,y:=1,z:=2)
143  \return the location (x,y,z) (in angstroms) of a given voxel. If the
144  index is not part of the map, the function returns -1.
145  */
146  float get_location_in_dim_by_voxel(long index,int dim) const;
147 
148 
149  //! Calculate the voxel of a given xyz indexes
150  /** \param[in] x The voxel index on the x axis of the grid
151  \param[in] y The voxel index on the y axis of the grid
152  \param[in] z The voxel index on the z axis of the grid
153  \return the voxel index.
154  */
155  inline long xyz_ind2voxel(int x,int y,int z) const {
156  return z * header_.get_nx() * header_.get_ny() +
157  y * header_.get_nx() + x;
158 
159  }
160 
161  //! Calculate the voxel of a given location
162  /** \param[in] x The position (in angstroms) of the x coordinate
163  \param[in] y The position (in angstroms) of the y coordinate
164  \param[in] z The position (in angstroms) of the z coordinate
165  \return the voxel index of a given position. If the position is out of
166  the boundaries of the map, the function returns -1.
167  */
168  long get_voxel_by_location(float x, float y, float z) const;
169 
170  //! Calculate the voxel of a given location
171  /** \param[in] v The position (in angstroms)
172  \return the voxel index of a given position. If the position is out of
173  the boundaries of the map, the function returns -1.
174  */
176  return get_voxel_by_location(v[0],v[1],v[2]);
177  }
178  //! Calculate dimension index of a given location
179  /** \param[in] v The position (in angstroms)
180  \param[in] ind dimension index (X:0,Y:1 or Z:2)
181  */
182  int get_dim_index_by_location(const algebra::Vector3D &v,int ind) const;
183 
184  //! Calculate the location of a given voxel.
185  /** \param[in] index The voxel index
186  \return the location (x,y,z) (in angstroms) of a given voxel.
187  */
189  IMP_USAGE_CHECK(index >= 0 && index < get_number_of_voxels(),
190  "invalid map index");
191  IMP_USAGE_CHECK(loc_calculated_,
192  "locations should be calculated prior to calling this function");
193  return algebra::Vector3D(x_loc_[index], y_loc_[index], z_loc_[index]);
194  }
195 
196  bool is_xyz_ind_part_of_volume(int ix,int iy,int iz) const;
197 
198  //! Checks whether a given point is in the grid the voxel of a given location
199  /** \param[in] x The position ( in angstroms) of the x coordinate
200  \param[in] y The position ( in angstroms) of the y coordinate
201  \param[in] z The position ( in angstroms) of the z coordinate
202  \return true if the point is part of the grid, false otherwise.
203  */
204  bool is_part_of_volume(float x,float y,float z) const;
205  //! Checks whether a given point is in the grid the voxel of a given location
206  /** \param[in] v The position ( in angstroms)
207  \return true if the point is part of the grid, false otherwise.
208  */
209  bool is_part_of_volume(const algebra::Vector3D &v) const {
210  return is_part_of_volume(v[0],v[1],v[2]);
211  }
212 
213  //! Gets the value of the voxel located at (x,y,z)
214  /** \param[in] x The position ( in angstroms) of the x coordinate
215  \param[in] y The position ( in angstroms) of the y coordinate
216  \param[in] z The position ( in angstroms) of the z coordinate
217  \return the value of the voxel located at (x,y,z)
218  \exception IndexException The point is not covered by the grid.
219  */
220  emreal get_value(float x,float y,float z) const;
221  emreal get_value(const algebra::Vector3D &point) const {
222  return get_value(point[0],point[1],point[2]);
223  }
224 
225 
226  //! Gets the value of the voxel at a given index
227  /**
228  \param[in] index voxel number in physical sense, NOT logical
229  */
230  emreal get_value(long index) const;
231 
232 
233  //! Set the value of the voxel at a given index
234  /**
235  \param[in] index voxel number in physical sense, NOT logical
236  \param[in] value value
237  */
238  void set_value(long index,emreal value);
239 
240 
241  //! Set the value of the voxel at a given index
242  /**
243  index voxel number in physical sense, NOT logical
244  */
245  void set_value(float x, float y, float z,emreal value);
246 
247  //! Sets the origin of the header
248  /**
249  \param x the new x (angstroms)
250  \param y the new y (angstroms)
251  \param z the new z (angstroms)
252  */
253  void set_origin(float x,float y,float z);
254  void set_origin(const IMP::algebra::Vector3D &v) {
255  set_origin(v[0],v[1],v[2]);
256  }
257 
258  algebra::Vector3D get_origin() const{
259  return algebra::Vector3D(header_.get_origin(0),
260  header_.get_origin(1),
261  header_.get_origin(2));
262  }
263 
264  algebra::Vector3D get_top() const {
265  return algebra::Vector3D(header_.get_top(0),
266  header_.get_top(1),
267  header_.get_top(2));
268  }
269 
270  // inspection functions
271  const DensityHeader *get_header()const {return &header_;}
272  //! Returns a pointer to the header of the map in a writable version
273  DensityHeader *get_header_writable() {return &header_;}
274 
275 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
276  //! Returns the x-location of the map
277  /**
278  \exception InvalidStateException The locations have not been calculated.
279  */
280  float* get_x_loc() const {
281  IMP_USAGE_CHECK(loc_calculated_,
282  "x location requested before being calculated");
283  return x_loc_.get();
284  }
285  //! Returns the y-location of the map
286  /**
287  \exception InvalidStateException The locations have not been calculated.
288  */
289  float* get_y_loc() const {
290  IMP_USAGE_CHECK(loc_calculated_,
291  "y location requested before being calculated");
292  return y_loc_.get();
293  }
294  //! Returns the z-location of the map
295  /**
296  \exception InvalidStateException The locations have not been calculated.
297  */
298  float* get_z_loc() const {
299  IMP_USAGE_CHECK(loc_calculated_,
300  "z location requested before being calculated");
301  return z_loc_.get();
302  }
303 
304  emreal* get_data() const {return data_.get();}
305 #endif
306 
307  //! Checks if two maps have the same origin
308  /** \param[in] other the map to compare with
309  \return true if the two maps have the same origin
310  */
311  bool same_origin(const DensityMap *other) const;
312 
313  //! Checks if two maps have the same dimensions
314  /** \param[in] other the map to compare with
315  \return true if the two maps have the same dimensions
316  */
317  bool same_dimensions(const DensityMap *other) const;
318 
319  //! Checks if two maps have the same voxel size
320  /** \param[in] other the map to compare with
321  \return true if the two maps have the same voxel size
322  */
323  bool same_voxel_size(const DensityMap *other) const;
324  //! Calculates the centroid of all the voxels with
325  //! density above a given threshold
326  /** \param[in] threshold the input threshold
327  */
328  algebra::Vector3D get_centroid(emreal threshold=0.0) const;
329  //! Returns the the value of the voxel with the highest density.
330  emreal get_max_value() const;
331  //! Returns the the value of the voxel with the lowest density.
332  emreal get_min_value() const;
333  //! Sums two grids.
334  //! The result is kept in the map.
335  /** \param[in] other the other map
336  \note The shared extend is sumed
337  \note The two maps should have the same voxelsize and other
338  should be contained within this map
339  */
340  void add(const DensityMap *other);
341 
342  //! Pick the max value between two corresponding voxels between two maps
343  //! The result is kept in the map.
344  /** \param[in] other the other map
345  \note The two maps should have the same voxelsize and the same dimensions
346  */
347  void pick_max(const DensityMap *other);
348 
349  //! number of map voxels
350  long get_number_of_voxels() const;
351 
352  //! Set the map dimension and reset all voxels to 0
353  /**
354  \param[in] nx x-dimension (voxels)
355  \param[in] ny y-dimension (voxels)
356  \param[in] nz z-dimension (voxels)
357  \note the origin and spacing remain unchanged
358  */
359  void set_void_map(int nx,int ny,int nz);
360 
361 
362  //! Increase the dimension of the map
363  //! The function pads zeros to the right-upper section on the map.
364  //! The original content of the map will be in the lower XYZ part of the map
365  /** \param[in] nx the number of voxels on the X axis
366  \param[in] ny the number of voxels on the Y axis
367  \param[in] nz the number of voxels on the Z axis
368  \param[in] val all additional voxels will have this value
369  \exception if the input x/y/z voxels is smaller than the one
370  currently in the map
371  */
372  void pad(int nx, int ny, int nz,float val=0.0);
373 
374  //! Create a new padded map
375  /** \brief Given this map of size [nx,ny,nz],
376  the new map is of size [2*mrg_x+nx,2*mrg_y+ny,2*mrg_z+nz].
377  The new map will consist of the values of the old map,
378  padded margin on all sides.
379  \param[in] mrg_x
380  number of margin voxels to add on both right and left on the X axis
381  \param[in] mrg_y
382  number of margin voxels to add on both right and left on the Y axis
383  \param[in] mrg_z
384  number of margin voxels to add on both right and left on the Z axis
385  \param[in] val ignored
386  \exception if the input x/y/z voxels is smaller than the one
387  currently in the map
388  */
389  DensityMap* pad_margin(int mrg_x, int mrg_y, int mrg_z,float val=0.0);
390 
391 
392  //! Create a new cropped map
393  /** \brief The margins are determined to be the bounding box
394  with density values below the input
395  \param[in] threshold used for cropping
396  */
397  DensityMap* get_cropped(float threshold);
398  //! Create a new cropped map with the bounding box extent
399  /**
400  \param[in] bb the bounding box
401  \note If the input bounding box is larger than the density box,
402  it is snapped to the right size.
403  */
404  DensityMap* get_cropped(const algebra::BoundingBox3D &bb);
405 
406  //! Get the maximum value in a XY plane indicated by a Z index
407  float get_maximum_value_in_xy_plane(int z_ind);
408  //! Get the maximum value in a XZ plane indicated by a Y index
409  float get_maximum_value_in_xz_plane(int y_ind);
410  //! Get the maximum value in a YZ plane indicated by a X index
411  float get_maximum_value_in_yz_plane(int x_ind);
412 
413  //! Multiply each voxel in the map by the input factor
414  //! The result is kept in the map.
415  /** \param[in] factor the multiplication factor
416  */
417  void multiply(float factor);
418 
419  //! Prints the locations of all of the voxels with value above a given
420  //! threshold into the input stream.
421  std::string get_locations_string(float t);
422 
423  //! Updated the voxel size of the map
424  void update_voxel_size(float new_apix);
425 
426  //! Updated the voxel size of the map
427  /** \note Use update_voxel_size() to set the spacing value.
428  */
429  Float get_spacing() const {return header_.get_spacing();}
430  //! Calculates the coordinates that correspond to all voxels.
431  void calc_all_voxel2loc();
432  //! copy map into this map
433  void copy_map(const DensityMap *other);
435 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
436  //! Convolution a kernel with a map and write results into current map
437  /**
438  \param[in] other the map to convolute
439  \param[in] kernel an array of kernel values. The data is in ZYX
440  order, Z is the slowest.
441  \param[in] dim_len the kernel array
442  */
443  void convolute_kernel(DensityMap *other,
444  double *kernel, int dim_len);
445  //! Convolution a kernel with this map and write results to this map
446  /**
447  \param[in] kernel an array of kernel values. The data is in ZYX
448  order, Z is the slowest.
449  \param[in] dim_len the kernel array
450  */
451  void convolute_kernel(
452  double *kernel, int dim_len) {
453  Pointer<DensityMap> cmap = em::create_density_map(this);
454  cmap->set_was_used(true);
455  convolute_kernel(cmap,kernel,dim_len);
456  cmap=static_cast<DensityMap*>(nullptr);
457  }
458 #endif
459  int lower_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
460  int upper_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
461  inline bool get_rms_calculated() const {return rms_calculated_;}
462  int get_dim_index_by_location(float loc_val,
463  int ind) const;
464 protected:
465  //!update the header values -- still in work
466  void update_header();
467  void reset_all_voxel2loc();
468 
469  void allocated_data();
470  void float2real(float *f_data, boost::scoped_array<emreal> &r_data);
471  void real2float(emreal *r_data, boost::scoped_array<float> &f_data);
472 
473  DensityHeader header_; // holds all the info about the map
474  boost::scoped_array<emreal> data_; // the order is ZYX (Z-slowest)
475  bool data_allocated_;
476 
477  //! Locations for each of the voxels of the map (they are precomputed and
478  //! each one is of size nvox, being nvox the size of the map)
479  boost::scoped_array<float> x_loc_, y_loc_, z_loc_;
480  //! true if the locations have already been computed
482 
483  bool normalized_;
484  bool rms_calculated_;
485 
486 };
487 
489  const DensityHeader *h=m->get_header();
491  m->get_origin(),
492  m->get_origin()+algebra::Vector3D(m->get_spacing()*h->get_nx(),
493  m->get_spacing()*h->get_ny(),
494  m->get_spacing()*h->get_nz()));
495 }
496 
497 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
498  //! Calculate a bounding box around a 3D point within the EM grid
499  /**
500 \param[in] d_map the density map
501 \param[in] x coordinates of the point to sample around
502 \param[in] y coordinates of the point to sample around
503 \param[in] z coordinates of the point to sample around
504 \param[in] kdist the lenght of the box
505 \param[out] iminx the minimum index on the X axis of the output bounding box
506 \param[out] iminy the minimum index on the Y axis of the output bounding box
507 \param[out] iminz the minimum index on the Z axis of the output bounding box
508 \param[out] imaxx the maximum index on the X axis of the output bounding box
509 \param[out] imaxy the maximum index on the Y axis of the output bounding box
510 \param[out] imaxz the maximum index on the Z axis of the output bounding box
511  */
512 inline void calc_local_bounding_box(
513  const em::DensityMap *d_map,
514  double x,double y, double z,
515  float kdist,
516  int &iminx,int &iminy, int &iminz,
517  int &imaxx,int &imaxy, int &imaxz) {
518  const DensityHeader *h=d_map->get_header();
519  iminx = d_map->lower_voxel_shift(x, kdist,h->get_xorigin(),h->get_nx());
520  iminy = d_map->lower_voxel_shift(y, kdist,h->get_yorigin(),h->get_ny());
521  iminz = d_map->lower_voxel_shift(z, kdist,h->get_zorigin(),h->get_nz());
522  imaxx = d_map->upper_voxel_shift(x, kdist,h->get_xorigin(),h->get_nx());
523  imaxy = d_map->upper_voxel_shift(y, kdist,h->get_yorigin(),h->get_ny());
524  imaxz = d_map->upper_voxel_shift(z, kdist,h->get_zorigin(),h->get_nz());
525 }
526 #endif
527 
528 //! rotate a grid
529 /**
530 /param[in] orig_dens the density map to rotate
531 /param[in] trans the transformation
532 \note this is a low resolution operation.
533 IMPEMEXPORT DensityMap* rotate_grid(const DensityMap *orig_dens,
534  const algebra::Transformation3D &trans);
535 */
536 
538 
539 /** Return the value for the density map, m, at point v, interpolating linearly
540  from the sample values. The resulting function is C0 over R3.
541  \relatesalso DensityMap
542 */
543 IMPEMEXPORT double get_density(const DensityMap *m,
544  const algebra::Vector3D &v);
545 
546 /** Return a new density map containing a rotated version of the old
547  one. Only voxels whose value is above threshold are considered when
548  computing the bounding box of the new map (set IMP::em::get_bounding_box()).
549  \relatesalso DensityMap
550 */
551 IMPEMEXPORT DensityMap* get_transformed(const DensityMap *input,
552  const algebra::Transformation3D &tr,
553  double threshold);
554 
555 /** Return a new density map containing a rotated version of the old
556  one. The dimension of the new map is the same as the old one.
557  \relatesalso DensityMap
558 */
559 IMPEMEXPORT DensityMap* get_transformed(DensityMap *input,
560  const algebra::Transformation3D &tr);
561 
562 
563 //! Get a resampled version of the map.
564 /** The spacing is multiplied by scaling.
565  That means, scaling values greater than 1 increase the voxel size.
566  \relatesalso DensityMap
567 */
568 IMPEMEXPORT DensityMap* get_resampled(DensityMap *input, double scaling);
569 
570 
571 //! Rotate a density map into another map
572 /**
573 \param[in] source the map to transform
574 \param[in] tr transform the from density map by this transformation
575 \param[out] into the map to tranform into
576 \param[in] calc_rms if true RMS is calculated on the transformed map
577  \relatesalso DensityMap
578 */
579 IMPEMEXPORT void get_transformed_into(const DensityMap *source,
580  const algebra::Transformation3D &tr,
581  DensityMap *into,
582  bool calc_rms=true);
583 IMPEMEXPORT void get_transformed_into2(const DensityMap *source,
584  const algebra::Transformation3D &tr,
585  DensityMap *into);
586 
587 inline bool get_interiors_intersect(const DensityMap *d1,
588  const DensityMap *d2){
590  get_bounding_box(d2));
591 }
592 #if 0
593 //! Get a histrogram of density values
594 /**
595 \param[in] dmap the density map to analyse
596 \param[in] threshold only add voxels with value above this threshold
597  to the histogram
598 \param[in] num_bins the number of bins to have in the histogram
599  \relatesalso DensityMap
600 */
601 // IMPEMEXPORT statistics::Histogram
602 // get_density_histogram(const DensityMap *dmap, float threshold,int num_bins);
603 #endif
604 
605 //! Get a segment of the map according to xyz indexes
606 /**
607 \note the output map will be cover
608 the region [[nx_start,nx_end],[]ny_start,ny_end,[nz_start,nz_end]]
609  */
610 IMPEMEXPORT DensityMap* get_segment(DensityMap *map_to_segment,
611  int nx_start,int nx_end,
612  int ny_start,int ny_end,
613  int nz_start,int nz_end);
614 
615 //! Get a segment of the map covered by the input points
616 IMPEMEXPORT DensityMap* get_segment(DensityMap *map_to_segment,
617  algebra::Vector3Ds vecs,float dist);
618 //! Get a segment of the map covered by another map
619 IMPEMEXPORT DensityMap* get_segment_by_masking(DensityMap *map_to_segment,
620  DensityMap *mask,
621  float mas_threshold);
622 
623 //! Return a map with 0 for all voxels below the
624 //! threshold and 1 for thoes above
625 /**
626 \param[in] orig_map the map to binarize
627 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
628 \param[in] reverse if true values above the threshold
629  are set to 0 and 1 otherwise
630  */
631 IMPEMEXPORT DensityMap* binarize(DensityMap *orig_map,
632  float threshold,bool reverse=false);
633 
634 
635 //! Return a map with 0 for all voxels below the
636 //! threshold
637 /**
638 \param[in] orig_map the map to binarize
639 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
640  */
641 IMPEMEXPORT DensityMap* get_threshold_map(DensityMap *orig_map,
642  float threshold);
643 
644 //! Return a density map for which voxel i contains the result of
645 //! m1[i]*m2[i]. The function assumes m1 and m2 are of the same dimensions.
646 IMPEMEXPORT DensityMap* multiply(const DensityMap *m1,
647  const DensityMap *m2);
648 //! Return a convolution between density maps m1 and m2.
649 //! The function assumes m1 and m2 are of the same dimensions.
650 IMPEMEXPORT double convolute(const DensityMap *m1,const DensityMap *m2);
651 //! Return the sum of all voxels
652 IMPEMEXPORT double get_sum(const DensityMap *m1);
653 //!Return a density map for which each voxel is the maximum value from
654 //!the input densities.
655 /**
656 \note all input maps should have the same extent
657  */
658 IMPEMEXPORT DensityMap* get_max_map(DensityMaps maps);
659 
660 //! Return a new map with an updated spacing and interpolate input
661 //! data accordinly
662 IMPEMEXPORT
663 DensityMap* interpolate_map(DensityMap *in_map,double new_spacing);
664 
665 /** Return a dense grid containing the voxels of the passed density map
666  as well as the same bounding box.
667 */
668 IMPEMEXPORT
671 get_grid(DensityMap *in_map);
672 
673 
674 /** Return a density map with the values taken from the grid.
675 */
676 IMPEMEXPORT
677 DensityMap* create_density_map(const algebra::GridD<3,
678  algebra::DenseGridStorageD<3, float>,
679  float >& grid);
680 //!Return a binaries density map with 1 for voxels that are internal
681 //in the input density map
682 IMPEMEXPORT
683 DensityMap* get_binarized_interior(DensityMap *dmap);
684 IMPEM_END_NAMESPACE
685 
686 #endif /* IMPEM_DENSITY_MAP_H */