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