IMP  2.2.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 boudning 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. Guess the file type from the
47  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 boudning 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 intructed 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 voxles 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 sumed
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 
461 inline algebra::BoundingBoxD<3> get_bounding_box(const DensityMap *m) {
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 
499 //! rotate a grid
500 /**
501 /param[in] orig_dens the density map to rotate
502 /param[in] trans the transformation
503 \note this is a low resolution operation.
504 IMPEMEXPORT DensityMap* rotate_grid(const DensityMap *orig_dens,
505  const algebra::Transformation3D &trans);
506 */
507 
509 
510 /** Return the value for the density map, m, at point v, interpolating linearly
511  from the sample values. The resulting function is C0 over R3.
512  See DensityMap
513 */
514 IMPEMEXPORT double get_density(const DensityMap *m, const algebra::Vector3D &v);
515 
516 /** Return a new density map containing a rotated version of the old
517  one. Only voxels whose value is above threshold are considered when
518  computing the bounding box of the new map (set IMP::em::get_bounding_box()).
519  See DensityMap
520 */
521 IMPEMEXPORT DensityMap *get_transformed(const DensityMap *input,
522  const algebra::Transformation3D &tr,
523  double threshold);
524 
525 /** Return a new density map containing a rotated version of the old
526  one. The dimension of the new map is the same as the old one.
527  See DensityMap
528 */
529 IMPEMEXPORT DensityMap *get_transformed(DensityMap *input,
530  const algebra::Transformation3D &tr);
531 
532 //! Get a resampled version of the map.
533 /** The spacing is multiplied by scaling.
534  That means, scaling values greater than 1 increase the voxel size.
535  See DensityMap
536 */
537 IMPEMEXPORT DensityMap *get_resampled(DensityMap *input, double scaling);
538 
539 //! Rotate a density map into another map
540 /**
541 \param[in] source the map to transform
542 \param[in] tr transform the from density map by this transformation
543 \param[out] into the map to tranform into
544 \param[in] calc_rms if true RMS is calculated on the transformed map
545  See DensityMap
546 */
547 IMPEMEXPORT void get_transformed_into(const DensityMap *source,
548  const algebra::Transformation3D &tr,
549  DensityMap *into, bool calc_rms = true);
550 IMPEMEXPORT void get_transformed_into2(const DensityMap *source,
551  const algebra::Transformation3D &tr,
552  DensityMap *into);
553 
554 inline bool get_interiors_intersect(const DensityMap *d1,
555  const DensityMap *d2) {
556  return get_interiors_intersect(get_bounding_box(d1), get_bounding_box(d2));
557 }
558 #if 0
559 //! Get a histrogram of density values
560 /**
561 \param[in] dmap the density map to analyse
562 \param[in] threshold only add voxels with value above this threshold
563  to the histogram
564 \param[in] num_bins the number of bins to have in the histogram
565  See DensityMap
566 */
567 // IMPEMEXPORT statistics::Histogram
568 // get_density_histogram(const DensityMap *dmap, float threshold,int num_bins);
569 #endif
570 
571 //! Get a segment of the map according to xyz indexes
572 /**
573 \note the output map will be cover
574 the region [[nx_start,nx_end],[]ny_start,ny_end,[nz_start,nz_end]]
575  */
576 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment, int nx_start,
577  int nx_end, int ny_start, int ny_end,
578  int nz_start, int nz_end);
579 
580 //! Get a segment of the map covered by the input points
581 IMPEMEXPORT DensityMap *get_segment(DensityMap *map_to_segment,
582  algebra::Vector3Ds vecs, float dist);
583 //! Get a segment of the map covered by another map
584 IMPEMEXPORT DensityMap *get_segment_by_masking(DensityMap *map_to_segment,
585  DensityMap *mask,
586  float mas_threshold);
587 
588 //! Return a map with 0 for all voxels below the
589 //! threshold and 1 for thoes above
590 /**
591 \param[in] orig_map the map to binarize
592 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
593 \param[in] reverse if true values above the threshold
594  are set to 0 and 1 otherwise
595  */
596 IMPEMEXPORT DensityMap *binarize(DensityMap *orig_map, float threshold,
597  bool reverse = false);
598 
599 //! Return a map with 0 for all voxels below the
600 //! threshold
601 /**
602 \param[in] orig_map the map to binarize
603 \param[in] threshold values below the threshold are set to 0 and 1 otherwise
604  */
605 IMPEMEXPORT DensityMap *get_threshold_map(const DensityMap *orig_map,
606  float threshold);
607 
608 //! Return a density map for which voxel i contains the result of
609 //! m1[i]*m2[i]. The function assumes m1 and m2 are of the same dimensions.
610 IMPEMEXPORT DensityMap *multiply(const DensityMap *m1, const DensityMap *m2);
611 //! Return a convolution between density maps m1 and m2.
612 //! The function assumes m1 and m2 are of the same dimensions.
613 IMPEMEXPORT double convolute(const DensityMap *m1, const DensityMap *m2);
614 //! Return the sum of all voxels
615 IMPEMEXPORT double get_sum(const DensityMap *m1);
616 //!Return a density map for which each voxel is the maximum value from
617 //!the input densities.
618 /**
619 \note all input maps should have the same extent
620  */
621 IMPEMEXPORT DensityMap *get_max_map(DensityMaps maps);
622 
623 //! Return a new map with an updated spacing and interpolate input
624 //! data accordinly
625 IMPEMEXPORT
626 DensityMap *interpolate_map(DensityMap *in_map, double new_spacing);
627 
628 /** Return a dense grid containing the voxels of the passed density map
629  as well as the same bounding box.
630 */
631 IMPEMEXPORT
633  DensityMap *in_map);
634 
635 /** Create a density map from an arbitrary IMP::algebra::GridD */
636 template <class S, class V, class E>
638  const IMP::algebra::GridD<3, S, V, E> &arg) {
640  typedef IMP::algebra::GridD<3, S, V, E> Grid;
642  std::abs(arg.get_unit_cell()[0] - arg.get_unit_cell()[1]) < .01,
643  "The passed grid does not seem to have cubic voxels");
645  algebra::get_bounding_box(arg), arg.get_unit_cell()[0]);
646  IMP_USAGE_CHECK(arg.get_number_of_voxels(0) ==
647  static_cast<unsigned int>(ret->get_header()->get_nx()),
648  "X voxels don't match");
649  IMP_USAGE_CHECK(arg.get_number_of_voxels(1) ==
650  static_cast<unsigned int>(ret->get_header()->get_ny()),
651  "Y voxels don't match");
652  IMP_USAGE_CHECK(arg.get_number_of_voxels(2) ==
653  static_cast<unsigned int>(ret->get_header()->get_nz()),
654  "Z voxels don't match");
655  IMP_FOREACH(typename Grid::Index i, arg.get_all_indexes()) {
656  long vi = ret->xyz_ind2voxel(i[0], i[1], i[2]);
657  ret->set_value(vi, arg[vi]);
658  }
659  return ret.release();
660 }
661 
662 /** Return a density map with the values taken from the grid.
663 */
664 IMPEMEXPORT
665 DensityMap *create_density_map(
667 
668 //!Return a binaries density map with 1 for voxels that are internal
669 // in the input density map
670 IMPEMEXPORT
671 DensityMap *get_binarized_interior(DensityMap *dmap);
672 
673 /** Rasterize the particles into an existing density map. */
674 IMPEMEXPORT
675 void add_to_map(DensityMap *dm, const kernel::Particles &pis); // defined in
676 // SampledDensityMap.cpp
677 IMPEM_END_NAMESPACE
678 
679 #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
rotate a grid
Definition: DensityMap.h:508
DensityMap * create_density_map(const algebra::GridD< 3, algebra::DenseGridStorageD< 3, float >, float > &grid)
BoundingBoxD< 3 > BoundingBox3D
Definition: BoundingBoxD.h:160
#define IMP_FUNCTION_LOG
Beginning logging for a non-member function.
A nullptr-initialized pointer to an IMP Object.
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)
double get_sum(const DensityMap *m1)
Return the sum of all voxels.
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: base/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)
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
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.
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.
#define IMP_NO_SWIG(x)
Hide the line when SWIG is compiled or parses it.
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)
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
DensityMap * read_map(std::string filename)
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.
Common base class for heavy weight IMP objects.
Definition: base/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.
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing sets of objects.
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.
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: base/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
All grids that are in the python API should be defined here.
int get_nz() const
Get the number of voxels in the z dimension.
A shared base class to help in debugging and things.
DensityMap * get_transformed(DensityMap *input, const algebra::Transformation3D &tr)
DensityMap * get_resampled(DensityMap *input, double scaling)
Get a resampled version of the map.
BoundingBoxD< D > get_bounding_box(const BoundingBoxD< D > &g)
Definition: BoundingBoxD.h:160