IMP  2.3.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-2014 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPEM_DENSITY_MAP_H
10 #define IMPEM_DENSITY_MAP_H
11 #include <IMP/base/Pointer.h>
12 #include <IMP/em/em_config.h>
13 #include "DensityHeader.h"
14 #include "MapReaderWriter.h"
15 #include <IMP/base/Object.h>
16 #include <IMP/algebra/Vector3D.h>
19 #include <IMP/base/Object.h>
20 #include <boost/scoped_array.hpp>
21 #include <iostream>
22 #include <iomanip>
24 #include <IMP/kernel/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::base::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  */
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  /**
210  \param[in] index voxel number in physical sense, NOT logical
211  */
212  emreal get_value(long index) const;
213 
214  //! Set the value of the voxel at a given index
215  /**
216  \param[in] index voxel number in physical sense, NOT logical
217  \param[in] value value
218  */
219  void set_value(long index, emreal value);
220 
221  //! Set the value of the voxel at a given index
222  /**
223  index voxel number in physical sense, NOT logical
224  */
225  void set_value(float x, float y, float z, emreal value);
226 
227  //! Sets the origin of the header
228  /**
229  \param x the new x (angstroms)
230  \param y the new y (angstroms)
231  \param z the new z (angstroms)
232  */
233  void set_origin(float x, float y, float z);
234  void set_origin(const IMP::algebra::Vector3D &v) {
235  set_origin(v[0], v[1], v[2]);
236  }
237 
238  algebra::Vector3D get_origin() const {
239  return algebra::Vector3D(header_.get_origin(0), header_.get_origin(1),
240  header_.get_origin(2));
241  }
242 
243  algebra::Vector3D get_top() const {
244  return algebra::Vector3D(header_.get_top(0), header_.get_top(1),
245  header_.get_top(2));
246  }
247 
248  // inspection functions
249  const DensityHeader *get_header() const { return &header_; }
250  //! Returns a pointer to the header of the map in a writable version
251  DensityHeader *get_header_writable() { return &header_; }
252 
253 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
254  //! Returns the x-location of the map
255  /**
256  \exception InvalidStateException The locations have not been calculated.
257  */
258  float *get_x_loc() const {
259  IMP_USAGE_CHECK(loc_calculated_,
260  "x location requested before being calculated");
261  return x_loc_.get();
262  }
263  //! Returns the y-location of the map
264  /**
265  \exception InvalidStateException The locations have not been calculated.
266  */
267  float *get_y_loc() const {
268  IMP_USAGE_CHECK(loc_calculated_,
269  "y location requested before being calculated");
270  return y_loc_.get();
271  }
272  //! Returns the z-location of the map
273  /**
274  \exception InvalidStateException The locations have not been calculated.
275  */
276  float *get_z_loc() const {
277  IMP_USAGE_CHECK(loc_calculated_,
278  "z location requested before being calculated");
279  return z_loc_.get();
280  }
281 
282  emreal *get_data() const { return data_.get(); }
283 #endif
284 
285  //! Checks if two maps have the same origin
286  /** \param[in] other the map to compare with
287  \return true if the two maps have the same origin
288  */
289  bool same_origin(const DensityMap *other) const;
290 
291  //! Checks if two maps have the same dimensions
292  /** \param[in] other the map to compare with
293  \return true if the two maps have the same dimensions
294  */
295  bool same_dimensions(const DensityMap *other) const;
296 
297  //! Checks if two maps have the same voxel size
298  /** \param[in] other the map to compare with
299  \return true if the two maps have the same voxel size
300  */
301  bool same_voxel_size(const DensityMap *other) const;
302  //! Calculates the centroid of all the voxels with
303  //! density above a given threshold
304  /** \param[in] threshold the input threshold
305  */
306  algebra::Vector3D get_centroid(emreal threshold = 0.0) const;
307  //! Returns the value of the voxel with the highest density.
308  emreal get_max_value() const;
309  //! Returns the value of the voxel with the lowest density.
310  emreal get_min_value() const;
311  //! Sums two grids.
312  //! The result is kept in the map.
313  /** \param[in] other the other map
314  \note The shared extend is summed
315  \note The two maps should have the same voxelsize and other
316  should be contained within this map
317  */
318  void add(const DensityMap *other);
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  //! number of map voxels
328  long get_number_of_voxels() const;
329 
330  //! Set the map dimension and reset all voxels to 0
331  /**
332  \param[in] nx x-dimension (voxels)
333  \param[in] ny y-dimension (voxels)
334  \param[in] nz z-dimension (voxels)
335  \note the origin and spacing remain unchanged
336  */
337  void set_void_map(int nx, int ny, int nz);
338 
339  //! Increase the dimension of the map
340  //! The function pads zeros to the right-upper section on the map.
341  //! The original content of the map will be in the lower XYZ part of the map
342  /** \param[in] nx the number of voxels on the X axis
343  \param[in] ny the number of voxels on the Y axis
344  \param[in] nz the number of voxels on the Z axis
345  \param[in] val all additional voxels will have this value
346  \exception if the input x/y/z voxels is smaller than the one
347  currently in the map
348  */
349  void pad(int nx, int ny, int nz, float val = 0.0);
350 
351  //! Create a new padded map
352  /** \brief Given this map of size [nx,ny,nz],
353  the new map is of size [2*mrg_x+nx,2*mrg_y+ny,2*mrg_z+nz].
354  The new map will consist of the values of the old map,
355  padded margin on all sides.
356  \param[in] mrg_x
357  number of margin voxels to add on both right and left on the X axis
358  \param[in] mrg_y
359  number of margin voxels to add on both right and left on the Y axis
360  \param[in] mrg_z
361  number of margin voxels to add on both right and left on the Z axis
362  \param[in] val ignored
363  \exception if the input x/y/z voxels is smaller than the one
364  currently in the map
365  */
366  DensityMap *pad_margin(int mrg_x, int mrg_y, int mrg_z, float val = 0.0);
367 
368  //! Create a new cropped map
369  /** \brief The margins are determined to be the bounding box
370  with density values below the input
371  \param[in] threshold used for cropping
372  */
373  DensityMap *get_cropped(float threshold);
374  //! Create a new cropped map with the bounding box extent
375  /**
376  \param[in] bb the bounding box
377  \note If the input bounding box is larger than the density box,
378  it is snapped to the right size.
379  */
380  DensityMap *get_cropped(const algebra::BoundingBox3D &bb);
381 
382  //! Get the maximum value in a XY plane indicated by a Z index
383  float get_maximum_value_in_xy_plane(int z_ind);
384  //! Get the maximum value in a XZ plane indicated by a Y index
385  float get_maximum_value_in_xz_plane(int y_ind);
386  //! Get the maximum value in a YZ plane indicated by a X index
387  float get_maximum_value_in_yz_plane(int x_ind);
388 
389  //! Multiply each voxel in the map by the input factor
390  //! The result is kept in the map.
391  /** \param[in] factor the multiplication factor
392  */
393  void multiply(float factor);
394 
395  //! Prints the locations of all of the voxels with value above a given
396  //! threshold into the input stream.
397  std::string get_locations_string(float t);
398 
399  //! Updated the voxel size of the map
400  void update_voxel_size(float new_apix);
401 
402  //! Updated the voxel size of the map
403  /** \note Use update_voxel_size() to set the spacing value.
404  */
405  Float get_spacing() const { return header_.get_spacing(); }
406  //! Calculates the coordinates that correspond to all voxels.
407  void calc_all_voxel2loc();
408  //! copy map into this map
409  void copy_map(const DensityMap *other);
411 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
412  //! Convolution a kernel with a map and write results into current map
413  /**
414  \param[in] other the map to convolute
415  \param[in] kernel an array of kernel values. The data is in ZYX
416  order, Z is the slowest.
417  \param[in] dim_len the kernel array
418  */
419  void convolute_kernel(DensityMap *other, double *kernel, int dim_len);
420  //! Convolution a kernel with this map and write results to this map
421  /**
422  \param[in] kernel an array of kernel values. The data is in ZYX
423  order, Z is the slowest.
424  \param[in] dim_len the kernel array
425  */
426  void convolute_kernel(double *kernel, int dim_len) {
428  cmap->set_was_used(true);
429  convolute_kernel(cmap, kernel, dim_len);
430  cmap = static_cast<DensityMap *>(nullptr);
431  }
432 #endif
433  int lower_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
434  int upper_voxel_shift(emreal loc, emreal kdist, emreal orig, int ndim) const;
435  inline bool get_rms_calculated() const { return rms_calculated_; }
436  int get_dim_index_by_location(float loc_val, int ind) const;
437 
438  protected:
439  //!update the header values -- still in work
440  void update_header();
441  void reset_all_voxel2loc();
442 
443  void allocated_data();
444  void float2real(float *f_data, boost::scoped_array<emreal> &r_data);
445  void real2float(emreal *r_data, boost::scoped_array<float> &f_data);
446 
447  DensityHeader header_; // holds all the info about the map
448  boost::scoped_array<emreal> data_; // the order is ZYX (Z-slowest)
449  bool data_allocated_;
450 
451  //! Locations for each of the voxels of the map (they are precomputed and
452  //! each one is of size nvox, being nvox the size of the map)
453  boost::scoped_array<float> x_loc_, y_loc_, z_loc_;
454  //! true if the locations have already been computed
456 
457  bool normalized_;
458  bool rms_calculated_;
459 };
460 
462  const DensityHeader *h = m->get_header();
464  m->get_origin(),
465  m->get_origin() + algebra::Vector3D(m->get_spacing() * h->get_nx(),
466  m->get_spacing() * h->get_ny(),
467  m->get_spacing() * h->get_nz()));
468 }
469 
470 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
471 //! Calculate a bounding box around a 3D point within the EM grid
472 /**
473 \param[in] d_map the density map
474 \param[in] x coordinates of the point to sample around
475 \param[in] y coordinates of the point to sample around
476 \param[in] z coordinates of the point to sample around
477 \param[in] kdist the length of the box
478 \param[out] iminx the minimum index on the X axis of the output bounding box
479 \param[out] iminy the minimum index on the Y axis of the output bounding box
480 \param[out] iminz the minimum index on the Z axis of the output bounding box
481 \param[out] imaxx the maximum index on the X axis of the output bounding box
482 \param[out] imaxy the maximum index on the Y axis of the output bounding box
483 \param[out] imaxz the maximum index on the Z axis of the output bounding box
484  */
485 inline void calc_local_bounding_box(const em::DensityMap *d_map, double x,
486  double y, double z, float kdist, int &iminx,
487  int &iminy, int &iminz, int &imaxx,
488  int &imaxy, int &imaxz) {
489  const DensityHeader *h = d_map->get_header();
490  iminx = d_map->lower_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
491  iminy = d_map->lower_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
492  iminz = d_map->lower_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
493  imaxx = d_map->upper_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
494  imaxy = d_map->upper_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
495  imaxz = d_map->upper_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
496 }
497 #endif
498 
500 
501 /** Return the value for the density map, m, at point v, interpolating linearly
502  from the sample values. The resulting function is C0 over R3.
503  \see DensityMap
504 */
505 IMPEMEXPORT double get_density(const DensityMap *m, const algebra::Vector3D &v);
506 
507 //! Return a new density map containing a rotated version of the old one.
508 /** Only voxels whose value is above threshold are considered when
509  computing the bounding box of the new map (set IMP::em::get_bounding_box()).
510  \see DensityMap
511 */
512 IMPEMEXPORT DensityMap *get_transformed(const DensityMap *input,
513  const algebra::Transformation3D &tr,
514  double threshold);
515 
516 //! Return a new density map containing a rotated version of the old one.
517 /** The dimension of the new map is the same as the old one.
518  \see DensityMap
519 */
520 IMPEMEXPORT DensityMap *get_transformed(DensityMap *input,
521  const algebra::Transformation3D &tr);
522 
523 //! Get a resampled version of the map.
524 /** The spacing is multiplied by scaling.
525  That means, scaling values greater than 1 increase the voxel size.
526  \see DensityMap
527 */
528 IMPEMEXPORT DensityMap *get_resampled(DensityMap *input, double scaling);
529 
530 //! Rotate a density map into another map
531 /**
532 \param[in] source the map to transform
533 \param[in] tr transform the from density map by this transformation
534 \param[out] into the map to transform into
535 \param[in] calc_rms if true RMS is calculated on the transformed map
536  \see DensityMap
537 */
538 IMPEMEXPORT void get_transformed_into(const DensityMap *source,
539  const algebra::Transformation3D &tr,
540  DensityMap *into, bool calc_rms = true);
541 IMPEMEXPORT void get_transformed_into2(const DensityMap *source,
542  const algebra::Transformation3D &tr,
543  DensityMap *into);
544 
545 inline bool get_interiors_intersect(const DensityMap *d1,
546  const DensityMap *d2) {
547  return get_interiors_intersect(get_bounding_box(d1), get_bounding_box(d2));
548 }
549 #if 0
550 //! Get a histogram of density values
551 /**
552 \param[in] dmap the density map to analyze
553 \param[in] threshold only add voxels with value above this threshold
554  to the histogram
555 \param[in] num_bins the number of bins to have in the histogram
556  \see DensityMap
557 */
558 // IMPEMEXPORT statistics::Histogram
559 // get_density_histogram(const DensityMap *dmap, float threshold,int num_bins);
560 #endif
561 
562 //! Get a segment of the map according to xyz indexes
563 /**
564 \note the output map will be cover
565 the region [[nx_start,nx_end],[]ny_start,ny_end,[nz_start,nz_end]]
566  */
567 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment, int nx_start,
568  int nx_end, int ny_start, int ny_end,
569  int nz_start, int nz_end);
570 
571 //! Get a segment of the map covered by the input points
572 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment,
573  algebra::Vector3Ds vecs, float dist);
574 //! Get a segment of the map covered by another map
575 IMPEMEXPORT DensityMap *get_segment_by_masking(DensityMap *map_to_segment,
576  DensityMap *mask,
577  float mas_threshold);
578 
579 //! Return a map with 0 for all voxels below the
580 //! threshold and 1 for those above
581 /**
582 \param[in] orig_map the map to binarize
583 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
584 \param[in] reverse if true values above the threshold
585  are set to 0 and 1 otherwise
586  */
587 IMPEMEXPORT DensityMap *binarize(DensityMap *orig_map, float threshold,
588  bool reverse = false);
589 
590 //! Return a map with 0 for all voxels below the
591 //! threshold
592 /**
593 \param[in] orig_map the map to binarize
594 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
595  */
596 IMPEMEXPORT DensityMap *get_threshold_map(const DensityMap *orig_map,
597  float threshold);
598 
599 //! Return a density map for which voxel i contains the result of
600 //! m1[i]*m2[i]. The function assumes m1 and m2 are of the same dimensions.
601 IMPEMEXPORT DensityMap *multiply(const DensityMap *m1, const DensityMap *m2);
602 //! Return a convolution between density maps m1 and m2.
603 //! The function assumes m1 and m2 are of the same dimensions.
604 IMPEMEXPORT double convolute(const DensityMap *m1, const DensityMap *m2);
605 //! Return the sum of all voxels
606 IMPEMEXPORT double get_sum(const DensityMap *m1);
607 //!Return a density map for which each voxel is the maximum value from
608 //!the input densities.
609 /**
610 \note all input maps should have the same extent
611  */
612 IMPEMEXPORT DensityMap *get_max_map(DensityMaps maps);
613 
614 //! Return a new map with an updated spacing and interpolate input
615 //! data accordingly
616 IMPEMEXPORT
617 DensityMap *interpolate_map(DensityMap *in_map, double new_spacing);
618 
619 /** Return a dense grid containing the voxels of the passed density map
620  as well as the same bounding box.
621 */
622 IMPEMEXPORT
624  DensityMap *in_map);
625 
626 /** Create a density map from an arbitrary IMP::algebra::GridD */
627 template <class S, class V, class E>
629  const IMP::algebra::GridD<3, S, V, E> &arg) {
631  typedef IMP::algebra::GridD<3, S, V, E> Grid;
633  std::abs(arg.get_unit_cell()[0] - arg.get_unit_cell()[1]) < .01,
634  "The passed grid does not seem to have cubic voxels");
636  algebra::get_bounding_box(arg), arg.get_unit_cell()[0]);
637  IMP_USAGE_CHECK(arg.get_number_of_voxels(0) ==
638  static_cast<unsigned int>(ret->get_header()->get_nx()),
639  "X voxels don't match");
640  IMP_USAGE_CHECK(arg.get_number_of_voxels(1) ==
641  static_cast<unsigned int>(ret->get_header()->get_ny()),
642  "Y voxels don't match");
643  IMP_USAGE_CHECK(arg.get_number_of_voxels(2) ==
644  static_cast<unsigned int>(ret->get_header()->get_nz()),
645  "Z voxels don't match");
646  IMP_FOREACH(typename Grid::Index i, arg.get_all_indexes()) {
647  long vi = ret->xyz_ind2voxel(i[0], i[1], i[2]);
648  ret->set_value(vi, arg[vi]);
649  }
650  return ret.release();
651 }
652 
653 /** Return a density map with the values taken from the grid.
654 */
655 IMPEMEXPORT
656 DensityMap *create_density_map(
658 
659 //!Return a binaries density map with 1 for voxels that are internal
660 // in the input density map
661 IMPEMEXPORT
662 DensityMap *get_binarized_interior(DensityMap *dmap);
663 
664 /** Rasterize the particles into an existing density map. */
665 IMPEMEXPORT
666 void add_to_map(DensityMap *dm, const kernel::Particles &pis); // defined in
667 // SampledDensityMap.cpp
668 IMPEM_END_NAMESPACE
669 
670 #endif /* IMPEM_DENSITY_MAP_H */
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)
IMP::base::Vector< IMP::base::Pointer< DensityMap > > DensityMaps
Definition: DensityMap.h:499
DensityMap * create_density_map(const algebra::GridD< 3, algebra::DenseGridStorageD< 3, float >, float > &grid)
BoundingBoxD< 3 > BoundingBox3D
Typedef for Python.
Definition: BoundingBoxD.h:160
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
void add_to_map(DensityMap *dm, const kernel::Particles &pis)
Basic types used by IMP.
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.
A smart pointer to a reference counted object.
Definition: Pointer.h:87
bool loc_calculated_
true if the locations have already been computed
Definition: DensityMap.h:455
O * release()
Relinquish control of the raw pointer stored in the Pointer.
void write_map(DensityMap *m, std::string filename)
Write a density map to a file.
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)
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:453
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:251
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
Common base class for heavy weight IMP objects.
Definition: Object.h:106
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:52
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:405
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:170
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(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.