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