IMP logo
IMP Reference Guide  2.5.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-2015 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 and return it.
42 /** \see 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  \see DensityMap
53 */
54 IMPEMEXPORT DensityMap *read_map(std::string filename);
55 
56 //! Write a density map to a file.
57 /** \see 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  \see DensityMap
70 */
71 IMPEMEXPORT void write_map(DensityMap *m, std::string filename);
72 
73 //!
74 /**
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 //!
82 /**
83 \param[in] m a density map
84 \param[in] threshold consider volume of only voxels above this threshold
85 \note The method assumes 1.21 cubic A per Dalton (Harpaz 1994).
86  */
87 IMPEMEXPORT Float approximate_molecular_mass(DensityMap *m, Float threshold);
88 
89 //! Class for handling density maps.
90 /** \note The location of a voxel is its center. That is important
91  for sampling function as well as for functions
92  like get_location_in_dim_by_voxel.
93  */
94 class IMPEMEXPORT DensityMap : public IMP::Object {
96  IMP_NO_SWIG(friend IMPEMEXPORT DensityMap *read_map(std::string filename,
97  MapReaderWriter *reader));
98  IMP_NO_SWIG(friend IMPEMEXPORT void write_map(DensityMap *m,
99  std::string filename,
100  MapReaderWriter *writer));
101 
102  public:
103  DensityMap(std::string name = "DensityMap%1%");
104  //! Construct a density map as instructed in the input header
105  DensityMap(const DensityHeader &header, std::string name = "DensityMap%1%");
106 
107  //! Set the density voxels to some value and reset the management flags.
108  /**
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 map is normalized.
115  */
116  emreal 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 indexes
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 the voxel of a given location
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 the voxel of a given location
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  emreal get_value(float x, float y, float z) const;
206  emreal 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  /**
212  \param[in] index voxel number in physical sense, NOT logical
213  */
214  emreal get_value(long index) const;
215 
216  //! Set the value of the voxel at a given index
217  /**
218  \param[in] index voxel number in physical sense, NOT logical
219  \param[in] value value
220  */
221  void set_value(long index, emreal value);
222 
223  //! Set the value of the voxel at a given index
224  /**
225  index voxel number in physical sense, NOT logical
226  */
227  void set_value(float x, float y, float z, emreal value);
228 
229  //! Sets the origin of the header
230  /**
231  \param x the new x (angstroms)
232  \param y the new y (angstroms)
233  \param z the new z (angstroms)
234  */
235  void set_origin(float x, float y, float z);
236  void set_origin(const IMP::algebra::Vector3D &v) {
237  set_origin(v[0], v[1], v[2]);
238  }
239 
240  algebra::Vector3D get_origin() const {
241  return algebra::Vector3D(header_.get_origin(0), header_.get_origin(1),
242  header_.get_origin(2));
243  }
244 
245  algebra::Vector3D get_top() const {
246  return algebra::Vector3D(header_.get_top(0), header_.get_top(1),
247  header_.get_top(2));
248  }
249 
250  // inspection functions
251  const DensityHeader *get_header() const { return &header_; }
252  //! Returns a pointer to the header of the map in a writable version
253  DensityHeader *get_header_writable() { return &header_; }
254 
255 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
256  //! Returns the x-location of the map
257  /**
258  \exception InvalidStateException The locations have not been calculated.
259  */
260  float *get_x_loc() const {
261  IMP_USAGE_CHECK(loc_calculated_,
262  "x location requested before being calculated");
263  return x_loc_.get();
264  }
265  //! Returns the y-location of the map
266  /**
267  \exception InvalidStateException The locations have not been calculated.
268  */
269  float *get_y_loc() const {
270  IMP_USAGE_CHECK(loc_calculated_,
271  "y location requested before being calculated");
272  return y_loc_.get();
273  }
274  //! Returns the z-location of the map
275  /**
276  \exception InvalidStateException The locations have not been calculated.
277  */
278  float *get_z_loc() const {
279  IMP_USAGE_CHECK(loc_calculated_,
280  "z location requested before being calculated");
281  return z_loc_.get();
282  }
283 
284  emreal *get_data() const { return data_.get(); }
285 #endif
286 
287  //! Checks if two maps have the same origin
288  /** \param[in] other the map to compare with
289  \return true if the two maps have the same origin
290  */
291  bool same_origin(const DensityMap *other) const;
292 
293  //! Checks if two maps have the same dimensions
294  /** \param[in] other the map to compare with
295  \return true if the two maps have the same dimensions
296  */
297  bool same_dimensions(const DensityMap *other) const;
298 
299  //! Checks if two maps have the same voxel size
300  /** \param[in] other the map to compare with
301  \return true if the two maps have the same voxel size
302  */
303  bool same_voxel_size(const DensityMap *other) const;
304  //! Calculates the centroid of all the voxels with
305  //! density above a given threshold
306  /** \param[in] threshold the input threshold
307  */
308  algebra::Vector3D get_centroid(emreal threshold = 0.0) const;
309  //! Returns the value of the voxel with the highest density.
310  emreal get_max_value() const;
311  //! Returns the value of the voxel with the lowest density.
312  emreal get_min_value() const;
313  //! Sums two grids.
314  //! The result is kept in the map.
315  /** \param[in] other the other map
316  \note The shared extend is summed
317  \note The two maps should have the same voxelsize and other
318  should be contained within this map
319  */
320  void add(const DensityMap *other);
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  //! number of map voxels
330  long get_number_of_voxels() const;
331 
332  //! Set the map dimension and reset all voxels to 0
333  /**
334  \param[in] nx x-dimension (voxels)
335  \param[in] ny y-dimension (voxels)
336  \param[in] nz z-dimension (voxels)
337  \note the origin and spacing remain unchanged
338  */
339  void set_void_map(int nx, int ny, int nz);
340 
341  //! Increase the dimension of the map
342  //! The function pads zeros to the right-upper section on the map.
343  //! The original content of the map will be in the lower XYZ part of the map
344  /** \param[in] nx the number of voxels on the X axis
345  \param[in] ny the number of voxels on the Y axis
346  \param[in] nz the number of voxels on the Z axis
347  \param[in] val all additional voxels will have this value
348  \exception if the input x/y/z voxels is smaller than the one
349  currently in the map
350  */
351  void pad(int nx, int ny, int nz, float val = 0.0);
352 
353  //! Create a new padded map
354  /** \brief Given this map of size [nx,ny,nz],
355  the new map is of size [2*mrg_x+nx,2*mrg_y+ny,2*mrg_z+nz].
356  The new map will consist of the values of the old map,
357  padded margin on all sides.
358  \param[in] mrg_x
359  number of margin voxels to add on both right and left on the X axis
360  \param[in] mrg_y
361  number of margin voxels to add on both right and left on the Y axis
362  \param[in] mrg_z
363  number of margin voxels to add on both right and left on the Z axis
364  \param[in] val ignored
365  \exception if the input x/y/z voxels is smaller than the one
366  currently in the map
367  */
368  DensityMap *pad_margin(int mrg_x, int mrg_y, int mrg_z, float val = 0.0);
369 
370  //! Create a new cropped map
371  /** \brief The margins are determined to be the bounding box
372  with density values below the input
373  \param[in] threshold used for cropping
374  */
375  DensityMap *get_cropped(float threshold);
376  //! Create a new cropped map with the bounding box extent
377  /**
378  \param[in] bb the bounding box
379  \note If the input bounding box is larger than the density box,
380  it is snapped to the right size.
381  */
382  DensityMap *get_cropped(const algebra::BoundingBox3D &bb);
383 
384  //! Get the maximum value in a XY plane indicated by a Z index
385  float get_maximum_value_in_xy_plane(int z_ind);
386  //! Get the maximum value in a XZ plane indicated by a Y index
387  float get_maximum_value_in_xz_plane(int y_ind);
388  //! Get the maximum value in a YZ plane indicated by a X index
389  float get_maximum_value_in_yz_plane(int x_ind);
390 
391  //! Multiply each voxel in the map by the input factor
392  //! The result is kept in the map.
393  /** \param[in] factor the multiplication factor
394  */
395  void multiply(float factor);
396 
397  //! Prints the locations of all of the voxels with value above a given
398  //! threshold into the input stream.
399  std::string get_locations_string(float t);
400 
401  //! Updated the voxel size of the map
402  void update_voxel_size(float new_apix);
403 
404  //! Updated 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  //! Convolution 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  //! Convolution 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 /** Return the value for the density map, m, at point v, interpolating linearly
508  from the sample values. The resulting function is C0 over R3.
509  \see 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  \see 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  \see 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.
531  That means, scaling values greater than 1 increase the voxel size.
532  \see DensityMap
533 */
534 IMPEMEXPORT DensityMap *get_resampled(DensityMap *input, double scaling);
535 
536 //! Rotate a density map into another map
537 /**
538 \param[in] source the map to transform
539 \param[in] tr transform the from density map by this transformation
540 \param[out] into the map to transform into
541 \param[in] calc_rms if true RMS is calculated on the transformed map
542  \see DensityMap
543 */
544 IMPEMEXPORT void get_transformed_into(const DensityMap *source,
545  const algebra::Transformation3D &tr,
546  DensityMap *into, bool calc_rms = true);
547 IMPEMEXPORT void get_transformed_into2(const DensityMap *source,
548  const algebra::Transformation3D &tr,
549  DensityMap *into);
550 
551 inline bool get_interiors_intersect(const DensityMap *d1,
552  const DensityMap *d2) {
553  return get_interiors_intersect(get_bounding_box(d1), get_bounding_box(d2));
554 }
555 #if 0
556 //! Get a histogram of density values
557 /**
558 \param[in] dmap the density map to analyze
559 \param[in] threshold only add voxels with value above this threshold
560  to the histogram
561 \param[in] num_bins the number of bins to have in the histogram
562  \see DensityMap
563 */
564 // IMPEMEXPORT statistics::Histogram
565 // get_density_histogram(const DensityMap *dmap, float threshold,int num_bins);
566 #endif
567 
568 //! Get a segment of the map according to xyz indexes
569 /**
570 \note the output map will be cover
571 the region [[nx_start,nx_end],[]ny_start,ny_end,[nz_start,nz_end]]
572  */
573 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment, int nx_start,
574  int nx_end, int ny_start, int ny_end,
575  int nz_start, int nz_end);
576 
577 //! Get a segment of the map covered by the input points
578 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment,
579  algebra::Vector3Ds vecs, float dist);
580 //! Get a segment of the map covered by another map
581 IMPEMEXPORT DensityMap *get_segment_by_masking(DensityMap *map_to_segment,
582  DensityMap *mask,
583  float mas_threshold);
584 
585 //! Return a map with 0 for all voxels below the
586 //! threshold and 1 for those above
587 /**
588 \param[in] orig_map the map to binarize
589 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
590 \param[in] reverse if true values above the threshold
591  are set to 0 and 1 otherwise
592  */
593 IMPEMEXPORT DensityMap *binarize(DensityMap *orig_map, float threshold,
594  bool reverse = false);
595 
596 //! Return a map with 0 for all voxels below the
597 //! threshold
598 /**
599 \param[in] orig_map the map to binarize
600 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
601  */
602 IMPEMEXPORT DensityMap *get_threshold_map(const DensityMap *orig_map,
603  float threshold);
604 
605 //! Return a density map for which voxel i contains the result of
606 //! m1[i]*m2[i]. The function assumes m1 and m2 are of the same dimensions.
607 IMPEMEXPORT DensityMap *multiply(const DensityMap *m1, const DensityMap *m2);
608 //! Return a convolution between density maps m1 and m2.
609 //! The function assumes m1 and m2 are of the same dimensions.
610 IMPEMEXPORT double convolute(const DensityMap *m1, const DensityMap *m2);
611 //! Return the sum of all voxels
612 IMPEMEXPORT double get_sum(const DensityMap *m1);
613 //!Return a density map for which each voxel is the maximum value from
614 //!the input densities.
615 /**
616 \note all input maps should have the same extent
617  */
618 IMPEMEXPORT DensityMap *get_max_map(DensityMaps maps);
619 
620 //! Return a new map with an updated spacing and interpolate input
621 //! data accordingly
622 IMPEMEXPORT
623 DensityMap *interpolate_map(DensityMap *in_map, double new_spacing);
624 
625 /** Return a dense grid containing the voxels of the passed density map
626  as well as the same bounding box.
627 */
628 IMPEMEXPORT
630  DensityMap *in_map);
631 
632 /** Create a density map from an arbitrary IMP::algebra::GridD */
633 template <class S, class V, class E>
635  const IMP::algebra::GridD<3, S, V, E> &arg) {
637  typedef IMP::algebra::GridD<3, S, V, E> Grid;
639  std::abs(arg.get_unit_cell()[0] - arg.get_unit_cell()[1]) < .01,
640  "The passed grid does not seem to have cubic voxels");
642  algebra::get_bounding_box(arg), arg.get_unit_cell()[0]);
643  IMP_USAGE_CHECK(arg.get_number_of_voxels(0) ==
644  static_cast<unsigned int>(ret->get_header()->get_nx()),
645  "X voxels don't match");
646  IMP_USAGE_CHECK(arg.get_number_of_voxels(1) ==
647  static_cast<unsigned int>(ret->get_header()->get_ny()),
648  "Y voxels don't match");
649  IMP_USAGE_CHECK(arg.get_number_of_voxels(2) ==
650  static_cast<unsigned int>(ret->get_header()->get_nz()),
651  "Z voxels don't match");
652  IMP_FOREACH(typename Grid::Index i, arg.get_all_indexes()) {
653  long vi = ret->xyz_ind2voxel(i[0], i[1], i[2]);
654  ret->set_value(vi, arg[vi]);
655  }
656  return ret.release();
657 }
658 
659 /** Return a density map with the values taken from the grid.
660 */
661 IMPEMEXPORT
662 DensityMap *create_density_map(
664 
665 //!Return a binaries density map with 1 for voxels that are internal
666 // in the input density map
667 IMPEMEXPORT
668 DensityMap *get_binarized_interior(DensityMap *dmap);
669 
670 /** Rasterize the particles into an existing density map. */
671 IMPEMEXPORT
672 void add_to_map(DensityMap *dm, const Particles &pis); // defined in
673 // SampledDensityMap.cpp
674 IMPEM_END_NAMESPACE
675 
676 #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
void get_transformed_into(const DensityMap *source, const algebra::Transformation3D &tr, DensityMap *into, bool calc_rms=true)
Rotate a density map into another map.
The base class to handle reading and writing of density maps.
Simple 3D transformation class.
DensityMap * interpolate_map(DensityMap *in_map, double new_spacing)
Basic types used by IMP.
DensityMap * create_density_map(const algebra::GridD< 3, algebra::DenseGridStorageD< 3, float >, float > &grid)
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)
#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.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
Definition: object_macros.h:25
Metadata for a density file.
A voxel grid in d-dimensional space space.
Definition: GridD.h:79
DensityMap * get_threshold_map(const DensityMap *orig_map, float 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.
bool loc_calculated_
true if the locations have already been computed
Definition: DensityMap.h:458
void write_map(DensityMap *m, std::string filename)
Write a density map to a file.
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)
Class for handling density maps.
Definition: DensityMap.h:94
double convolute(const DensityMap *m1, const DensityMap *m2)
Common base class for heavy weight IMP objects.
Definition: Object.h:106
DensityMap * multiply(const DensityMap *m1, const DensityMap *m2)
An abstract class for reading a map.
BoundingBoxD< 3 > get_bounding_box(const Cone3D &g)
Definition: Cone3D.h:64
DensityMap * get_binarized_interior(DensityMap *dmap)
Return a binaries density map with 1 for voxels that are internal.
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:170
double get_density(const DensityMap *m, const algebra::Vector3D &v)
DensityMap * read_map(std::string filename)
Read a density map from a file and return it.
DensityMap * get_max_map(DensityMaps 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:253
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
bool is_part_of_volume(const algebra::Vector3D &v) const
Checks whether a given point is in the grid the voxel of a given location.
Definition: DensityMap.h:192
Float approximate_molecular_mass(DensityMap *m, Float threshold)
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
Updated the voxel size of the map.
Definition: DensityMap.h:407
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 indexes.
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.
void add_to_map(DensityMap *dm, const Particles &pis)
A dense grid of values.
int get_nz() const
Get the number of voxels in the z dimension.
DensityMap * get_transformed(DensityMap *input, const algebra::Transformation3D &tr)
Return a new density map containing a rotated version of the old one.
DensityMap * get_resampled(DensityMap *input, double scaling)
Get a resampled version of the map.