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