IMP logo
IMP Reference Guide  2.14.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-2020 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
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
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 
330  //! Set the map dimension and reset all voxels to 0
331  /** \param[in] nx x-dimension (voxels)
332  \param[in] ny y-dimension (voxels)
333  \param[in] nz z-dimension (voxels)
334  \note the origin and spacing remain unchanged
335  */
336  void set_void_map(int nx, int ny, int nz);
337 
338  //! Increase the dimension of the map
339  /** The function pads zeroes to the right-upper section on the map.
340  The original content of the map will be in the lower XYZ part of the map.
341  \param[in] nx the number of voxels on the X axis
342  \param[in] ny the number of voxels on the Y axis
343  \param[in] nz the number of voxels on the Z axis
344  \param[in] val all additional voxels will have this value
345  \exception UsageException if the input numbers of x/y/z voxels are
346  smaller than those currently in the map
347  */
348  void pad(int nx, int ny, int nz, float val = 0.0);
349 
350  //! Create a new padded map
351  /** Given this map of size [nx,ny,nz],
352  the new map is of size [2*mrg_x+nx,2*mrg_y+ny,2*mrg_z+nz].
353  The new map will consist of the values of the old map,
354  with a padded margin on all sides.
355 
356  \param[in] mrg_x number of margin voxels to add on both right and
357  left on the X axis
358  \param[in] mrg_y number of margin voxels to add on both right and
359  left on the Y axis
360  \param[in] mrg_z number of margin voxels to add on both right and
361  left on the Z axis
362  \param[in] val ignored
363  \return the new padded map.
364  \exception UsageException if the input numbers of x/y/z voxels are
365  smaller than those currently in the map
366  */
367  DensityMap *pad_margin(int mrg_x, int mrg_y, int mrg_z, float val = 0.0);
368 
369  //! Create and return a new cropped map
370  /** The margins are determined to be the bounding box
371  with density values below the input
372  \param[in] threshold used for cropping
373  \return the new cropped map.
374  */
375  DensityMap *get_cropped(float threshold);
376 
377  //! Create and return a new cropped map with the given bounding box extent
378  /** \param[in] bb the bounding box
379  \return the new cropped map.
380  \note If the input bounding box is larger than the density box,
381  it is snapped to the right size.
382  */
383  DensityMap *get_cropped(const algebra::BoundingBox3D &bb);
384 
385  //! Get the maximum value in a XY plane indicated by a Z index
386  float get_maximum_value_in_xy_plane(int z_ind);
387  //! Get the maximum value in a XZ plane indicated by a Y index
388  float get_maximum_value_in_xz_plane(int y_ind);
389  //! Get the maximum value in a YZ plane indicated by a X index
390  float get_maximum_value_in_yz_plane(int x_ind);
391 
392  //! Multiply in place each voxel in the map by the input factor
393  /** The result is kept in the map.
394  \param[in] factor the multiplication factor
395  */
396  void multiply(float factor);
397 
398  //! Print the locations of all voxels with value above a given threshold
399  std::string get_locations_string(float t);
400 
401  //! Update the voxel size of the map
402  void update_voxel_size(float new_apix);
403 
404  //! Get the voxel size of the map
405  /** \note Use update_voxel_size() to set the spacing value.
406  */
407  Float get_spacing() const { return header_.get_spacing(); }
408  //! Calculates the coordinates that correspond to all voxels.
409  void calc_all_voxel2loc();
410  //! Copy map into this map
411  void copy_map(const DensityMap *other);
413 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
414  //! Convolute a kernel with a map and write results into current map
415  /**
416  \param[in] other the map to convolute
417  \param[in] kernel an array of kernel values. The data is in ZYX
418  order, Z is the slowest.
419  \param[in] dim_len the kernel array
420  */
421  void convolute_kernel(DensityMap *other, double *kernel, int dim_len);
422  //! Convolute a kernel with this map and write results to this map
423  /**
424  \param[in] kernel an array of kernel values. The data is in ZYX
425  order, Z is the slowest.
426  \param[in] dim_len the kernel array
427  */
428  void convolute_kernel(double *kernel, int dim_len) {
430  cmap->set_was_used(true);
431  convolute_kernel(cmap, kernel, dim_len);
432  cmap = static_cast<DensityMap *>(nullptr);
433  }
434 #endif
435  int lower_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
436  int upper_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
437  inline bool get_rms_calculated() const { return rms_calculated_; }
438  int get_dim_index_by_location(float loc_val, int ind) const;
439 
440  protected:
441  //!update the header values -- still in work
442  void update_header();
443  void reset_all_voxel2loc();
444 
445  void allocated_data();
446  void float2real(float *f_data, boost::scoped_array<emreal> &r_data);
447  void real2float(emreal *r_data, boost::scoped_array<float> &f_data);
448 
449  DensityHeader header_; // holds all the info about the map
450  boost::scoped_array<emreal> data_; // the order is ZYX (Z-slowest)
451  bool data_allocated_;
452 
453  //! Locations (centers) for each of the voxels of the map (they are
454  //! precomputed and each one is of size nvox, where nvox is the
455  //! size of the map)
456  boost::scoped_array<float> x_loc_, y_loc_, z_loc_;
457  //! true if the locations have already been computed
459 
460  bool normalized_;
461  bool rms_calculated_;
462 };
463 
465  const DensityHeader *h = m->get_header();
466  float hspace = m->get_spacing() / 2.0;
467  algebra::Vector3D lb = m->get_origin()
468  - algebra::Vector3D(hspace, hspace, hspace);
470  lb,
471  lb + algebra::Vector3D(m->get_spacing() * h->get_nx(),
472  m->get_spacing() * h->get_ny(),
473  m->get_spacing() * h->get_nz()));
474 }
475 
476 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
477 //! Calculate a bounding box around a 3D point within the EM grid
478 /**
479 \param[in] d_map the density map
480 \param[in] x coordinates of the point to sample around
481 \param[in] y coordinates of the point to sample around
482 \param[in] z coordinates of the point to sample around
483 \param[in] kdist the length of the box
484 \param[out] iminx the minimum index on the X axis of the output bounding box
485 \param[out] iminy the minimum index on the Y axis of the output bounding box
486 \param[out] iminz the minimum index on the Z axis of the output bounding box
487 \param[out] imaxx the maximum index on the X axis of the output bounding box
488 \param[out] imaxy the maximum index on the Y axis of the output bounding box
489 \param[out] imaxz the maximum index on the Z axis of the output bounding box
490  */
491 inline void calc_local_bounding_box(const em::DensityMap *d_map, double x,
492  double y, double z, float kdist, int &iminx,
493  int &iminy, int &iminz, int &imaxx,
494  int &imaxy, int &imaxz) {
495  const DensityHeader *h = d_map->get_header();
496  iminx = d_map->lower_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
497  iminy = d_map->lower_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
498  iminz = d_map->lower_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
499  imaxx = d_map->upper_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
500  imaxy = d_map->upper_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
501  imaxz = d_map->upper_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
502 }
503 #endif
504 
506 
507 //! Get density value at point v, interpolating linearly from the sample values.
508 /** The resulting function is C0 over R3.
509  \relates DensityMap
510 */
511 IMPEMEXPORT double get_density(const DensityMap *m, const algebra::Vector3D &v);
512 
513 //! Return a new density map containing a rotated version of the old one.
514 /** Only voxels whose value is above threshold are considered when
515  computing the bounding box of the new map (set IMP::em::get_bounding_box()).
516  \relates DensityMap
517 */
518 IMPEMEXPORT DensityMap *get_transformed(const DensityMap *input,
519  const algebra::Transformation3D &tr,
520  double threshold);
521 
522 //! Return a new density map containing a rotated version of the old one.
523 /** The dimension of the new map is the same as the old one.
524  \relates DensityMap
525 */
526 IMPEMEXPORT DensityMap *get_transformed(DensityMap *input,
527  const algebra::Transformation3D &tr);
528 
529 //! Get a resampled version of the map.
530 /** The spacing is multiplied by scaling, i.e. scaling values greater
531  than 1 increase the voxel size.
532  \relates DensityMap
533 */
534 IMPEMEXPORT DensityMap *get_resampled(DensityMap *input, double scaling);
535 
536 //! Rotate a density map into another map
537 /** \param[in] source the map to transform
538  \param[in] tr transform the from density map by this transformation
539  \param[out] into the map to transform into
540  \param[in] calc_rms if true RMS is calculated on the transformed map
541  \relates DensityMap
542 */
543 IMPEMEXPORT void get_transformed_into(const DensityMap *source,
544  const algebra::Transformation3D &tr,
545  DensityMap *into, bool calc_rms = true);
546 IMPEMEXPORT void get_transformed_into2(const DensityMap *source,
547  const algebra::Transformation3D &tr,
548  DensityMap *into);
549 
550 inline bool get_interiors_intersect(const DensityMap *d1,
551  const DensityMap *d2) {
552  return get_interiors_intersect(get_bounding_box(d1), get_bounding_box(d2));
553 }
554 #if 0
555 //! Get a histogram of density values
556 /**
557 \param[in] dmap the density map to analyze
558 \param[in] threshold only add voxels with value above this threshold
559  to the histogram
560 \param[in] num_bins the number of bins to have in the histogram
561  \see DensityMap
562 */
563 // IMPEMEXPORT statistics::Histogram
564 // get_density_histogram(const DensityMap *dmap, float threshold,int num_bins);
565 #endif
566 
567 //! Get a segment of the map according to xyz indexes
568 /** \note the output map will cover the region
569  [[nx_start,nx_end],[]ny_start,ny_end,[nz_start,nz_end]]
570  */
571 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment, int nx_start,
572  int nx_end, int ny_start, int ny_end,
573  int nz_start, int nz_end);
574 
575 //! Get a segment of the map covered by the input points
576 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment,
577  algebra::Vector3Ds vecs, float dist);
578 //! Get a segment of the map covered by another map
579 IMPEMEXPORT DensityMap *get_segment_by_masking(DensityMap *map_to_segment,
580  DensityMap *mask,
581  float mas_threshold);
582 
583 //! Return a map with 0 for all voxels below the threshold and 1 for those above
584 /** \param[in] orig_map the map to binarize
585  \param[in] threshold values below the threshold are set to 0 and 1 otherwise
586  \param[in] reverse if true values above the threshold
587  are set to 0 and 1 otherwise
588  */
589 IMPEMEXPORT DensityMap *binarize(DensityMap *orig_map, float threshold,
590  bool reverse = false);
591 
592 //! Return a map with 0 for all voxels below the threshold
593 /** \param[in] orig_map the map to binarize
594  \param[in] threshold values below the threshold are set to 0 and
595  1 otherwise
596  */
597 IMPEMEXPORT DensityMap *get_threshold_map(const DensityMap *orig_map,
598  float threshold);
599 
600 //! Return a density map for which voxel i contains the result of m1[i]*m2[i].
601 /** The function assumes m1 and m2 are of the same dimensions.
602  */
603 IMPEMEXPORT DensityMap *multiply(const DensityMap *m1, const DensityMap *m2);
604 
605 //! Return a convolution between density maps m1 and m2.
606 /** The function assumes m1 and m2 are of the same dimensions.
607  */
608 IMPEMEXPORT double convolute(const DensityMap *m1, const DensityMap *m2);
609 
610 //! Return the sum of all voxels
611 IMPEMEXPORT double get_sum(const DensityMap *m1);
612 
613 //! Return a map where each voxel is the maximum value from the input maps.
614 /** \note all input maps should have the same extent
615  */
616 IMPEMEXPORT DensityMap *get_max_map(DensityMaps maps);
617 
618 //! Return a new map with an updated spacing
619 /** Input data is interpolated accordingly.
620  */
621 IMPEMEXPORT
622 DensityMap *interpolate_map(DensityMap *in_map, double new_spacing);
623 
624 //! Return a dense grid containing the voxels of the passed density map
625 /** The returned grid has the same bounding box as the map.
626  */
627 IMPEMEXPORT
629  DensityMap *in_map);
630 
631 //! Create a density map from an arbitrary IMP::algebra::GridD
632 template <class S, class V, class E>
634  const IMP::algebra::GridD<3, S, V, E> &arg) {
636  typedef IMP::algebra::GridD<3, S, V, E> Grid;
638  std::abs(arg.get_unit_cell()[0] - arg.get_unit_cell()[1]) < .01,
639  "The passed grid does not seem to have cubic voxels");
641  algebra::get_bounding_box(arg), arg.get_unit_cell()[0]);
642  IMP_USAGE_CHECK(arg.get_number_of_voxels(0) ==
643  static_cast<unsigned int>(ret->get_header()->get_nx()),
644  "X voxels don't match");
645  IMP_USAGE_CHECK(arg.get_number_of_voxels(1) ==
646  static_cast<unsigned int>(ret->get_header()->get_ny()),
647  "Y voxels don't match");
648  IMP_USAGE_CHECK(arg.get_number_of_voxels(2) ==
649  static_cast<unsigned int>(ret->get_header()->get_nz()),
650  "Z voxels don't match");
651  IMP_FOREACH(typename Grid::Index i, arg.get_all_indexes()) {
652  long vi = ret->xyz_ind2voxel(i[0], i[1], i[2]);
653  ret->set_value(vi, arg[vi]);
654  }
655  return ret.release();
656 }
657 
658 //! Return a density map with the values taken from the grid.
659 /** \relates DensityMap
660  */
661 IMPEMEXPORT
663 
664 //! Return a density map with the values taken from the grid.
665 /** \relates DensityMap
666  */
667 IMPEMEXPORT
669 
670 
671 //! Return a binarized map with 1 for voxels that are internal in the input map
672 /** \relates DensityMap
673  */
674 IMPEMEXPORT
675 DensityMap *get_binarized_interior(DensityMap *dmap);
676 
677 //! Rasterize the particles into an existing density map.
678 /** \relates DensityMap
679  */
680 IMPEMEXPORT
681 void add_to_map(DensityMap *dm, const Particles &pis); // defined in
682  // SampledDensityMap.cpp
683 
684 IMPEM_END_NAMESPACE
685 
686 #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:505
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:458
A more IMP-like version of the std::vector.
Definition: Vector.h:39
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.
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:456
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
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:407
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.