9 #ifndef IMPEM_DENSITY_MAP_H
10 #define IMPEM_DENSITY_MAP_H
12 #include <IMP/em/em_config.h>
20 #include <boost/scoped_array.hpp>
21 #include <cereal/access.hpp>
22 #include <cereal/types/base_class.hpp>
47 IMPEMEXPORT DensityMap *
read_map(std::string filename, MapReaderWriter *reader);
57 IMPEMEXPORT DensityMap *
read_map(std::string filename);
62 IMPEMEXPORT
void write_map(DensityMap *m, std::string filename,
63 MapReaderWriter *writer);
74 IMPEMEXPORT
void write_map(DensityMap *m, std::string filename);
100 std::string filename,
104 DensityMap(std::string name =
"DensityMap%1%");
111 void reset_data(
float value = 0.0);
121 void std_normalize();
123 inline bool is_normalized()
const {
return normalized_; }
131 float get_location_in_dim_by_voxel(
long index,
int dim)
const;
140 return z * header_.get_nx() * header_.get_ny() + y * header_.get_nx() + x;
150 long get_voxel_by_location(
float x,
float y,
float z)
const;
158 return get_voxel_by_location(v[0], v[1], v[2]);
172 "invalid map index");
175 "locations should be calculated prior to calling this function");
179 bool is_xyz_ind_part_of_volume(
int ix,
int iy,
int iz)
const;
187 bool is_part_of_volume(
float x,
float y,
float z)
const;
193 return is_part_of_volume(v[0], v[1], v[2]);
205 double get_value(
float x,
float y,
float z)
const;
207 return get_value(point[0], point[1], point[2]);
213 double get_value(
long index)
const;
219 void set_value(
long index,
double value);
222 void set_value(
float x,
float y,
float z,
double value);
229 void set_origin(
float x,
float y,
float z);
231 set_origin(v[0], v[1], v[2]);
236 header_.get_origin(2));
250 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
255 float *get_x_loc()
const {
257 "x location requested before being calculated");
264 float *get_y_loc()
const {
266 "y location requested before being calculated");
273 float *get_z_loc()
const {
275 "z location requested before being calculated");
279 double *get_data()
const {
return data_.get(); }
286 bool same_origin(
const DensityMap *other)
const;
292 bool same_dimensions(
const DensityMap *other)
const;
298 bool same_voxel_size(
const DensityMap *other)
const;
304 double get_max_value()
const;
306 double get_min_value()
const;
314 void add(
const DensityMap *other);
327 void pick_max(
const DensityMap *other);
331 return header_.get_number_of_voxels();
340 void set_void_map(
int nx,
int ny,
int nz);
352 void pad(
int nx,
int ny,
int nz,
float val = 0.0);
371 DensityMap *pad_margin(
int mrg_x,
int mrg_y,
int mrg_z,
float val = 0.0);
394 DensityMap *get_cropped(
Particles ps,
double distance,
bool inverse =
false,
bool keep_map_dimensions =
false);
405 float get_maximum_value_in_xy_plane(
int z_ind);
407 float get_maximum_value_in_xz_plane(
int y_ind);
409 float get_maximum_value_in_yz_plane(
int x_ind);
418 std::string get_locations_string(
float t);
421 void update_voxel_size(
float new_apix);
428 void calc_all_voxel2loc();
432 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
440 void convolute_kernel(
DensityMap *other,
double *kernel,
int dim_len);
447 void convolute_kernel(
double *kernel,
int dim_len) {
449 cmap->set_was_used(
true);
450 convolute_kernel(cmap, kernel, dim_len);
454 int lower_voxel_shift(
double loc,
double kdist,
double orig,
int ndim)
const;
455 int upper_voxel_shift(
double loc,
double kdist,
double orig,
int ndim)
const;
456 inline bool get_rms_calculated()
const {
return rms_calculated_; }
457 int get_dim_index_by_location(
float loc_val,
int ind)
const;
461 void update_header();
462 void reset_all_voxel2loc();
464 void allocated_data();
465 void float2real(
float *f_data, boost::scoped_array<double> &r_data);
466 void real2float(
double *r_data, boost::scoped_array<float> &f_data);
468 DensityHeader header_;
469 boost::scoped_array<double> data_;
470 bool data_allocated_;
475 boost::scoped_array<float>
x_loc_, y_loc_, z_loc_;
480 bool rms_calculated_;
482 friend class cereal::access;
484 template<
class Archive>
void serialize(Archive &ar) {
485 ar(cereal::base_class<Object>(
this),
486 header_, data_allocated_, loc_calculated_, normalized_,
488 long size = get_number_of_voxels();
490 if (std::is_base_of<cereal::detail::InputArchiveBase, Archive>::value) {
491 data_.reset(
new double[size]);
492 if (loc_calculated_) {
494 loc_calculated_ =
false;
495 calc_all_voxel2loc();
499 for (
long i = 0; i < size; ++i) {
506 const DensityHeader *h = m->get_header();
507 float hspace = m->get_spacing() / 2.0;
510 return algebra::BoundingBoxD<3>(
513 m->get_spacing() * h->get_ny(),
514 m->get_spacing() * h->get_nz()));
517 #if !defined(IMP_DOXYGEN) && !defined(SWIG)
532 inline void calc_local_bounding_box(
const em::DensityMap *d_map,
double x,
533 double y,
double z,
float kdist,
int &iminx,
534 int &iminy,
int &iminz,
int &imaxx,
535 int &imaxy,
int &imaxz) {
536 const DensityHeader *h = d_map->get_header();
537 iminx = d_map->lower_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
538 iminy = d_map->lower_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
539 iminz = d_map->lower_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
540 imaxx = d_map->upper_voxel_shift(x, kdist, h->get_xorigin(), h->get_nx());
541 imaxy = d_map->upper_voxel_shift(y, kdist, h->get_yorigin(), h->get_ny());
542 imaxz = d_map->upper_voxel_shift(z, kdist, h->get_zorigin(), h->get_nz());
587 IMPEMEXPORT
void get_transformed_into2(
const DensityMap *source,
591 inline bool get_interiors_intersect(
const DensityMap *d1,
612 IMPEMEXPORT DensityMap *
get_segment(DensityMap *map_to_segment,
int nx_start,
613 int nx_end,
int ny_start,
int ny_end,
614 int nz_start,
int nz_end);
617 IMPEMEXPORT DensityMap *
get_segment(DensityMap *map_to_segment,
622 float mas_threshold);
630 IMPEMEXPORT DensityMap *
binarize(DensityMap *orig_map,
float threshold,
631 bool reverse =
false);
644 IMPEMEXPORT DensityMap *
multiply(
const DensityMap *m1,
const DensityMap *m2);
649 IMPEMEXPORT
double convolute(
const DensityMap *m1,
const DensityMap *m2);
652 IMPEMEXPORT
double get_sum(
const DensityMap *m1);
673 template <
class S,
class V,
class E>
679 std::abs(arg.get_unit_cell()[0] - arg.get_unit_cell()[1]) < .01,
680 "The passed grid does not seem to have cubic voxels");
682 algebra::get_bounding_box(arg), arg.get_unit_cell()[0]);
684 static_cast<unsigned int>(ret->get_header()->get_nx()),
685 "X voxels don't match");
687 static_cast<unsigned int>(ret->get_header()->get_ny()),
688 "Y voxels don't match");
690 static_cast<unsigned int>(ret->get_header()->get_nz()),
691 "Z voxels don't match");
692 for(
typename Grid::Index i : arg.get_all_indexes()) {
693 long vi = ret->xyz_ind2voxel(i[0], i[1], i[2]);
694 ret->set_value(vi, arg[vi]);
O * release()
Relinquish control of the raw pointer stored in the Pointer.
IMP::Vector< IMP::Pointer< DensityMap > > DensityMaps
The base class to handle reading and writing of density maps.
const DensityHeader * get_header() const
Returns a read-only pointer to the header of the map.
DensityMap * interpolate_map(DensityMap *in_map, double new_spacing)
Return a new map with an updated spacing.
BoundingBoxD< 3 > BoundingBox3D
Typedef for Python.
DensityMap * get_segment(DensityMap *map_to_segment, algebra::Vector3Ds vecs, float dist)
Get a segment of the map covered by the input points.
DensityMap * binarize(DensityMap *orig_map, float threshold, bool reverse=false)
Return a map with 0 for all voxels below the threshold and 1 for those above.
#define IMP_FUNCTION_LOG
Beginning logging for a non-member function.
double get_sum(const DensityMap *m1)
Return the sum of all voxels.
DensityMap * read_map(std::string filename, MapReaderWriter *reader)
Read a density map from a file (using the given reader) and return it.
#define IMP_OBJECT_METHODS(Name)
Define the basic things needed by any Object.
void add_to_map(DensityMap *dm, const Particles &pis)
Rasterize the particles into an existing density map.
A voxel grid in d-dimensional space.
DensityMap * get_threshold_map(const DensityMap *orig_map, float threshold)
Return a map with 0 for all voxels below the threshold.
long get_voxel_by_location(const algebra::Vector3D &v) const
Calculate the voxel of a given location.
DensityMap * get_segment_by_masking(DensityMap *map_to_segment, DensityMap *mask, float mas_threshold)
Get a segment of the map covered by another map.
DensityMap * create_density_map(algebra::DenseGrid3D< float > &grid)
Return a density map with the values taken from the grid.
bool loc_calculated_
true if the locations have already been computed
A smart pointer to a reference counted object.
algebra::GridD< 3, algebra::DenseGridStorageD< 3, float >, float > get_grid(DensityMap *in_map)
Return a dense grid containing the voxels of the passed density map.
Class for handling density maps.
double convolute(const DensityMap *m1, const DensityMap *m2)
Return a convolution between density maps m1 and m2.
DensityMap * get_resampled(DensityMap *input, double scaling)
Get a resampled version of the map.
void get_transformed_into(const DensityMap *source, const algebra::Transformation3D &tr, DensityMap *into, bool calc_rms=true)
Rotate a density map into another map.
void write_map(DensityMap *m, std::string filename, MapReaderWriter *writer)
Write a density map to a file using the given writer.
Common base class for heavy weight IMP objects.
DensityMap * multiply(const DensityMap *m1, const DensityMap *m2)
Return a density map for which voxel i contains the result of m1[i]*m2[i].
An abstract class for reading a map.
BoundingBoxD< 3 > get_bounding_box(const Cone3D &g)
boost::scoped_array< float > x_loc_
A bounding box in D dimensions.
algebra::Vector3D get_location_by_voxel(long index) const
Calculate the location of a given voxel.
double get_density(const DensityMap *m, const algebra::Vector3D &v)
Get density value at point v, interpolating linearly from the sample values.
DensityMap * get_max_map(DensityMaps maps)
Return a map where each voxel is the maximum value from the input maps.
DensityMap * create_density_map(const DensityMap *other)
Create a copy of another map.
DensityHeader * get_header_writable()
Returns a pointer to the header of the map in a writable version.
#define IMP_NO_SWIG(x)
Hide the line when SWIG is compiled or parses it.
long get_number_of_voxels() const
Get the number of map voxels.
DensityMap * get_binarized_interior(DensityMap *dmap)
Return a binarized map with 1 for voxels that are internal in the input map.
bool is_part_of_volume(const algebra::Vector3D &v) const
Checks whether a given point is in the grid.
Float approximate_molecular_mass(DensityMap *m, Float threshold)
Estimate the molecular mass from a map.
Vector3D get_centroid(const Vector3Ds &ps)
Return the centroid of a set of vectors.
#define IMP_OBJECTS(Name, PluralName)
Define the types for storing lists of object pointers.
A nullptr-initialized pointer to an IMP Object.
A shared base class to help in debugging and things.
Float get_spacing() const
Get the voxel size of the map.
double Float
Basic floating-point value (could be float, double...)
long xyz_ind2voxel(int x, int y, int z) const
Calculate the voxel of a given xyz index.
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
All grids that are in the Python API should be defined here.
DensityMap * get_transformed(const DensityMap *input, const algebra::Transformation3D &tr, double threshold)
Return a new density map containing a rotated version of the old one.