diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt
index 6c68c3fc4263a7f8cd634c633bb425337cfa68f8..30fae2eb199a4ff4f4c461eefcd5bfc4f00d12c8 100644
--- a/modules/mol/alg/pymod/CMakeLists.txt
+++ b/modules/mol/alg/pymod/CMakeLists.txt
@@ -7,6 +7,16 @@ set(OST_MOL_ALG_PYMOD_MODULES
"__init__.py"
views.py
)
+
+if (ENABLE_IMG)
+
+ set(OST_MOL_ALG_PYMOD_SOURCES
+ ${OST_MOL_ALG_PYMOD_SOURCES}
+ export_entity_to_density.cc
+ )
+
+endif()
+
pymod(NAME mol_alg OUTPUT_DIR ost/mol/alg CPP ${OST_MOL_ALG_PYMOD_SOURCES}
PY ${OST_MOL_ALG_PYMOD_MODULES})
diff --git a/modules/mol/alg/pymod/export_entity_to_density.cc b/modules/mol/alg/pymod/export_entity_to_density.cc
new file mode 100644
index 0000000000000000000000000000000000000000..858740eeb9500d82ddcaa57d45bcbd0296ce0c80
--- /dev/null
+++ b/modules/mol/alg/pymod/export_entity_to_density.cc
@@ -0,0 +1,74 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+//------------------------------------------------------------------------------
+
+/*
+ * Author Juergen Haas
+ */
+#include <boost/python.hpp>
+
+#include <ost/mol/alg/entity_to_density.hh>
+
+using namespace boost::python;
+using namespace ost;
+using namespace ost::mol::alg;
+
+
+#if OST_IMG_ENABLED
+#include <ost/mol/alg/entity_to_density.hh>
+
+//"thin wrappers" for default parameters
+int etd1(const mol::EntityView& entity_view,
+ img::MapHandle& map,
+ const DensityType& density_type,
+ Real falloff_start,
+ Real falloff_end,
+ bool clear_map_flag)
+{
+ EntityToDensity(entity_view, map, density_type, falloff_start,
+ falloff_end, clear_map_flag);
+ return 0;
+}
+
+int etd2(const mol::EntityView& entity_view,
+ img::MapHandle& map,
+ const DensityType& density_type,
+ Real falloff_start,
+ Real falloff_end)
+{
+ EntityToDensity(entity_view, map, density_type, falloff_start,
+ falloff_end);
+ return 0;
+}
+
+
+void export_entity_to_density()
+{
+
+ def("EntityToDensity",EntityToDensity);
+ def("EntityToDensity",etd1);
+ def("EntityToDensity",etd2);
+
+ enum_<DensityType>("DensityType")
+ .value("SCATTERING_FACTORS", SCATTERING_FACTORS)
+ .value("ROSETTA_HIGH_RESOLUTION", ROSETTA_HIGH_RESOLUTION)
+ .value("ROSETTA_LOW_RESOLUTION", ROSETTA_LOW_RESOLUTION)
+ .export_values()
+ ;
+ }
+#endif
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 1670b122e20cf6053bef0866d0dbb794ecb65315..14785ddf753a30194aa995ea0b195d56f08b9dcf 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -24,7 +24,12 @@
using namespace boost::python;
void export_svdSuperPose();
+void export_entity_to_density();
+
BOOST_PYTHON_MODULE(_mol_alg)
{
export_svdSuperPose();
+ #if OST_IMG_ENABLED
+ export_entity_to_density();
+ #endif
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 0eff2e75f86b2c1b0f911ae85e938d21d70f1fbd..adb666c1bf47005b5654c94aa4f01ed7b0ec6def 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -9,8 +9,24 @@ set(OST_MOL_ALG_SOURCES
sec_structure_segments.cc
)
+set(MOL_ALG_DEPS mol)
+
+if (ENABLE_IMG)
+ set(OST_MOL_ALG_HEADERS
+ ${OST_MOL_ALG_HEADERS}
+ entity_to_density.hh
+ )
+
+ set(OST_MOL_ALG_SOURCES
+ ${OST_MOL_ALG_SOURCES}
+ entity_to_density.cc
+ )
+
+ set(MOL_ALG_DEPS ${MOL_ALG_DEPS} img img_alg)
+endif()
+
module(NAME mol_alg SOURCES ${OST_MOL_ALG_SOURCES}
HEADERS ${OST_MOL_ALG_HEADERS}
HEADER_OUTPUT_DIR ost/mol/alg
- DEPENDS_ON mol)
+ DEPENDS_ON ${MOL_ALG_DEPS})
diff --git a/modules/mol/alg/src/entity_to_density.cc b/modules/mol/alg/src/entity_to_density.cc
new file mode 100644
index 0000000000000000000000000000000000000000..af52b64b17b11ac351e6e9ee3f1d57ae43eba0de
--- /dev/null
+++ b/modules/mol/alg/src/entity_to_density.cc
@@ -0,0 +1,418 @@
+#include <sstream>
+#include <cmath>
+
+#include <boost/filesystem/fstream.hpp>
+#include <boost/filesystem/convenience.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <ost/log.hh>
+#include <ost/platform.hh>
+#include <ost/img/alg/dft.hh>
+#include <ost/mol/mol.hh>
+
+#include "entity_to_density.hh"
+
+namespace ost { namespace mol { namespace alg {
+
+namespace detail {
+
+struct AtomScatteringProps {
+ String element;
+ Real atomic_weight;
+ Real c;
+ Real a1;
+ Real a2;
+ Real a3;
+ Real a4;
+ Real b1;
+ Real b2;
+ Real b3;
+ Real b4;
+};
+
+typedef std::vector<AtomScatteringProps> AtomScatteringPropsTable;
+
+void FillAtomScatteringPropsTable(AtomScatteringPropsTable& table)
+{
+ AtomScatteringProps atom_props;
+
+ String ost_root=GetSharedDataPath();
+ String filename = "/atom_scattering_properties.txt";
+ String fullpath = ost_root+filename;
+ boost::filesystem::path loc(fullpath);
+ boost::filesystem::ifstream infile(loc);
+ if (!infile) {
+ LOGN_ERROR("Couldn't find " << fullpath);
+ return;
+ }
+ String line;
+ std::getline(infile, line);
+ std::getline(infile, line);
+ std::getline(infile, line);
+ std::getline(infile, line);
+ std::getline(infile, line);
+ while (std::getline(infile, line))
+ {
+ std::stringstream line_stream(line);
+ line_stream >> atom_props.element
+ >> atom_props.atomic_weight
+ >> atom_props.c
+ >> atom_props.a1
+ >> atom_props.a2
+ >> atom_props.a3
+ >> atom_props.a4
+ >> atom_props.b1
+ >> atom_props.b2
+ >> atom_props.b3
+ >> atom_props.b4;
+ table.push_back(atom_props);
+ }
+ std::cout << table.size() << std::endl;
+};
+
+
+Real ScatteringFactor (Real frequency, const AtomScatteringProps& props,
+ Real source_wavelength)
+{
+ Real scattering_factor =
+ props.a1*exp(-props.b1*frequency*frequency) +
+ props.a2*exp(-props.b2*frequency*frequency) +
+ props.a3*exp(-props.b3*frequency*frequency) +
+ props.a4*exp(-props.b4*frequency*frequency) +
+ props.c;
+
+ return scattering_factor;
+}
+
+
+class EntityToDensityHelperBase
+{
+
+public:
+
+ EntityToDensityHelperBase (const ost::mol::EntityView entity_view,
+ const ost::mol::alg::DensityType& density_type,
+ Real falloff_start,
+ Real falloff_end,
+ Real source_wavelength):
+ density_type_(density_type), entity_view_(entity_view),
+ falloff_start_frequency_(1.0/falloff_start),
+ falloff_end_frequency_(1.0/falloff_end),
+ source_wavelength_(source_wavelength)
+ {
+ if (scattering_props_table_loaded_ == false)
+ {
+ FillAtomScatteringPropsTable(scatt_props_table_);
+ scattering_props_table_loaded_ = true;
+ }
+ }
+
+ void VisitState(img::ComplexHalfFrequencyImageState& is) const
+ {
+ geom::Vec3 map_start=is.IndexToCoord(is.GetExtent().GetStart());
+ geom::Vec3 map_end=is.IndexToCoord(is.GetExtent().GetEnd());
+ geom::Vec3 frequency_sampling =
+ is.GetSampling().GetFrequencySampling();
+
+ Real minus_2_pi = -2.0*M_PI;
+
+ uint x_limit = ceil(falloff_end_frequency_ / frequency_sampling[0]);
+ uint y_limit = ceil(falloff_end_frequency_ / frequency_sampling[1]);
+ uint z_limit = ceil(falloff_end_frequency_ / frequency_sampling[2]);
+ img::Extent reduced_extent = img::Extent
+ (img::Point(-x_limit,-y_limit,0),
+ img::Point(x_limit,y_limit,z_limit));
+
+ mol::AtomViewIter iterator_end = entity_view_.AtomsEnd();
+ for (mol::AtomViewIter iterator = entity_view_.AtomsBegin();
+ iterator!=iterator_end; ++iterator)
+ {
+ AtomScatteringPropsTable::iterator table_iter =
+ scatt_props_table_.begin();
+ bool found = false;
+ while (found != true && table_iter!=scatt_props_table_.end())
+ {
+ if ( (*table_iter).element == (*iterator).GetAtomProps().element)
+ {
+ geom::Vec3 coord = (*iterator).GetPos();
+
+ if (coord[0] >= map_start[0] &&
+ coord[0] <= map_end[0] &&
+ coord[1] >= map_start[1] &&
+ coord[1] <= map_end[1] &&
+ coord[2] >= map_start[2] &&
+ coord[2] <= map_end[2])
+ {
+
+ // This part of the code assumes that the three axes
+ // of the map are at right angles and the origin at 0,0,0.
+ // Eventually, when maps and images will be merged,
+ // it will substituted by the map's xyz2uvw method
+
+ geom::Vec3 adjusted_coord = coord-map_start;
+
+ for (img::ExtentIterator mp_it(reduced_extent);
+ !mp_it.AtEnd();++mp_it)
+ {
+ img::Point mp_it_point = img::Point(mp_it);
+
+ geom::Vec3 frequency_sampling =
+ is.GetSampling().GetFrequencySampling();
+
+ geom::Vec3 mp_it_vec = geom::Vec3(
+ (mp_it_point[0]*frequency_sampling[0]),
+ (mp_it_point[1]*frequency_sampling[1]),
+ (mp_it_point[2]*frequency_sampling[2])
+ );
+
+ Real frequency = Length(mp_it_vec);
+
+ Real sigma = (falloff_end_frequency_ -
+ falloff_start_frequency_)/3.0;
+
+ if (frequency <= falloff_end_frequency_)
+ {
+ Real falloff_term = 1.0;
+ if (sigma!=0 && frequency >= falloff_start_frequency_)
+ {
+ falloff_term=exp(-(frequency-falloff_start_frequency_)*
+ (frequency-falloff_start_frequency_)/
+ (2.0*sigma*sigma));
+ }
+ AtomScatteringProps scatt_props = (*table_iter);
+ Real scatt_fact = ScatteringFactor(frequency,
+ scatt_props,
+ source_wavelength_);
+
+ Real amp_term = scatt_fact*falloff_term;
+ Real exp = minus_2_pi * geom::Dot(mp_it_vec,adjusted_coord);
+ Real real = amp_term *std::cos(exp);
+ Real imag = amp_term *std::sin(exp);
+ is.Value(mp_it) = is.Value(mp_it)+Complex(real,imag);
+
+ }
+ }
+ }
+ }
+ ++table_iter;
+ }
+ }
+ }
+
+ template <typename T, class D>
+ void VisitState(img::ImageStateImpl<T,D>& is) const
+ {
+ assert(false);
+ }
+
+ static String GetAlgorithmName() {return "EntityToDensityHelper"; }
+
+private:
+
+ ost::mol::alg::DensityType density_type_;
+ Real limit_;
+ ost::mol::EntityView entity_view_;
+ Real falloff_start_frequency_;
+ Real falloff_end_frequency_;
+ static AtomScatteringPropsTable scatt_props_table_;
+ static bool scattering_props_table_loaded_;
+ Real source_wavelength_;
+};
+
+AtomScatteringPropsTable EntityToDensityHelperBase::scatt_props_table_ =
+ AtomScatteringPropsTable();
+bool EntityToDensityHelperBase::scattering_props_table_loaded_ = false;
+
+typedef img::ImageStateConstModIPAlgorithm<EntityToDensityHelperBase>
+ EntityToDensityHelper;
+
+
+
+class EntityToDensityRosettaHelperBase
+{
+
+public:
+
+ EntityToDensityRosettaHelperBase (const ost::mol::EntityView entity_view,
+ const ost::mol::alg::DensityType& density_type,
+ Real resolution):
+ density_type_(density_type), entity_view_(entity_view),
+ resolution_(resolution)
+ {
+ if (scattering_props_table_loaded_ == false)
+ {
+ FillAtomScatteringPropsTable(scatt_props_table_);
+ scattering_props_table_loaded_ = true;
+ }
+ }
+
+ void VisitState(img::RealSpatialImageState& is) const
+ {
+ geom::Vec3 map_start=is.IndexToCoord(is.GetExtent().GetStart());
+ geom::Vec3 map_end=is.IndexToCoord(is.GetExtent().GetEnd());
+ mol::EntityView effective_entity_view=entity_view_;
+ Real k,C;
+
+ if (density_type_ == ROSETTA_HIGH_RESOLUTION)
+ {
+ k = (M_PI*M_PI)/(resolution_*resolution_);
+ C = sqrt((M_PI*M_PI*M_PI)/(k*k*k));
+ effective_entity_view = entity_view_;
+
+ } else { // ROSETTA LOW RESOLUTION
+
+ k = (M_PI*M_PI)/((2.4+0.8*resolution_)*(2.4+0.8*resolution_));
+ C = sqrt((M_PI*M_PI*M_PI)/(k*k*k));
+
+ effective_entity_view = entity_view_.Select("aname=CA and peptide=true");
+ }
+
+ Real sigma = sqrt(1.0/(2.0*k));
+ Real four_sigma_squared = 16.0*sigma*sigma;
+ Real four_sigma = 4.0*sigma;
+
+ geom::Vec3 sampling = is.GetSampling().GetPixelSampling();
+ img::Extent is_extent = is.GetExtent();
+
+ mol::AtomViewIter iterator_end = effective_entity_view.AtomsEnd();
+ for (mol::AtomViewIter iterator = effective_entity_view.AtomsBegin();
+ iterator!=iterator_end; ++iterator) {
+ AtomScatteringPropsTable::iterator table_iter =
+ scatt_props_table_.begin();
+ bool found = false;
+ while (found != true && table_iter!=scatt_props_table_.end()) {
+ if ((*table_iter).element == (*iterator).GetAtomProps().element) {
+ found = true;
+ Real a = (*table_iter).atomic_weight;
+ geom::Vec3 coord = (*iterator).GetPos();
+ if (coord[0] >= map_start[0] &&
+ coord[0] <= map_end[0] &&
+ coord[1] >= map_start[1] &&
+ coord[1] <= map_end[1] &&
+ coord[2] >= map_start[2] &&
+ coord[2] <= map_end[2])
+ {
+ geom::Vec3 adjusted_coord = geom::Vec3(coord[0]-map_start[0],
+ coord[1]-map_start[1],
+ coord[2]-map_start[2]);
+ geom::Vec3 pixel_coord=is.CoordToIndex(coord);
+ img::Point rounded_pixel_coord(round(pixel_coord[0]),
+ round(pixel_coord[1]),
+ round(pixel_coord[2]));
+
+ uint x_limit = ceil(2.0*four_sigma/sampling[0])+1;
+ uint y_limit = ceil(2.0*four_sigma/sampling[1])+1;
+ uint z_limit = ceil(2.0*four_sigma/sampling[2])+1;
+
+ img::Extent reduced_extent = img::Extent
+ (img::Size(x_limit,y_limit,z_limit),
+ rounded_pixel_coord);
+
+ img::Extent iteration_extent=
+ img::Overlap(is_extent,reduced_extent);
+
+ for (img::ExtentIterator mp_it(iteration_extent);
+ !mp_it.AtEnd();++mp_it)
+ {
+
+ img::Point mp_it_point = img::Point(mp_it);
+
+ geom::Vec3 mp_it_vec = geom::Vec3(
+ (mp_it_point[0]*sampling[0]),
+ (mp_it_point[1]*sampling[1]),
+ (mp_it_point[2]*sampling[2])
+ );
+
+ Real distance_squared = Length2(mp_it_vec-adjusted_coord);
+
+ if (distance_squared <= four_sigma_squared)
+ {
+ Real exp_term = exp(-k*distance_squared);
+ Real value = C * a * exp_term;
+ is.Value(mp_it_point) = is.Value(mp_it_point) + value;
+ }
+ }
+ }
+ }
+ ++table_iter;
+ }
+ }
+ }
+
+ template <typename T, class D>
+ void VisitState(img::ImageStateImpl<T,D>& is) const
+ {
+ assert(false);
+ }
+
+ static String GetAlgorithmName() { return "EntityToDensityRosettaHelper"; }
+
+private:
+
+ ost::mol::alg::DensityType density_type_;
+ ost::mol::EntityView entity_view_;
+ Real resolution_;
+ static AtomScatteringPropsTable scatt_props_table_;
+ static bool scattering_props_table_loaded_;
+};
+
+AtomScatteringPropsTable EntityToDensityRosettaHelperBase::scatt_props_table_ =
+ AtomScatteringPropsTable();
+bool EntityToDensityRosettaHelperBase::scattering_props_table_loaded_ = false;
+
+typedef img::ImageStateConstModIPAlgorithm<EntityToDensityRosettaHelperBase>
+ EntityToDensityRosettaHelper;
+
+
+} // namespace
+
+void EntityToDensity (const mol::EntityView& entity_view,
+ img::MapHandle& map,
+ const DensityType& density_type,
+ Real falloff_start,
+ Real falloff_end,
+ bool clear_map_flag,
+ Real source_wavelength)
+{
+ if(falloff_start<=0.0) throw std::runtime_error("Invalid falloff start");
+ if(falloff_end<=0.0 || falloff_end>falloff_start)
+ throw std::runtime_error("Invalid falloff end");
+
+ if (density_type == SCATTERING_FACTORS)
+ {
+ if (clear_map_flag==true) {
+ img::MapHandle mm=img::CreateImage(img::Extent(map.GetSize(),
+ img::Point(0,0)),
+ img::COMPLEX,img::HALF_FREQUENCY);
+ // swap newly created map into place
+ mm.SetAbsoluteOrigin(map.GetAbsoluteOrigin());
+ map.Swap(mm);
+ } else {
+ map.ApplyIP(img::alg::DFT());
+ }
+
+ detail::EntityToDensityHelper e_to_d_helper(entity_view,density_type,
+ falloff_start,
+ falloff_end,
+ source_wavelength);
+ map.ApplyIP(e_to_d_helper);
+ map.ApplyIP(img::alg::DFT());
+ } else {
+ if (clear_map_flag==true) {
+ img::MapHandle mm=img::CreateImage(img::Extent(img::Point(0,0),
+ map.GetSize()));
+ // swap newly created map into place
+ mm.SetAbsoluteOrigin(map.GetAbsoluteOrigin());
+ map.Swap(mm);
+ }
+ Real resolution = (falloff_start+falloff_end)/2.0;
+ detail::EntityToDensityRosettaHelper e_to_d_r_helper
+ (entity_view,density_type,
+ resolution);
+ map.ApplyIP(e_to_d_r_helper);
+ }
+}
+
+}}} // ns
+
+
diff --git a/modules/mol/alg/src/entity_to_density.hh b/modules/mol/alg/src/entity_to_density.hh
new file mode 100644
index 0000000000000000000000000000000000000000..04cf2719619d43592a4eef5a64c72ab9761b8197
--- /dev/null
+++ b/modules/mol/alg/src/entity_to_density.hh
@@ -0,0 +1,64 @@
+#ifndef OST_ENTITY_TO_DENSITY_HH
+#define OST_ENTITY_TO_DENSITY_HH
+
+#import <ost/mol/entity_view.hh>
+#import <ost/img/map.hh>
+
+#include <ost/mol/alg/module_config.hh>
+
+
+namespace ost { namespace mol { namespace alg {
+
+/// \brief type of density being created by the EntityToDensity function
+enum DensityType
+{
+ SCATTERING_FACTORS,
+ ROSETTA_HIGH_RESOLUTION,
+ ROSETTA_LOW_RESOLUTION
+};
+
+/// \brief create a density representation of an entity in a density map
+///
+/// This function creates a density representation of the entity provided by
+/// the user in a density map, also provided by the user. The user can choose
+/// the type of density of the output map:
+///
+/// SCATTERING_FACTORS density, generated in Fourier space using the
+/// correct scattering factors for the elements
+/// involved.
+/// ROSETTA_HIGH_RESOLUTION gaussian spheres in real space to represent
+/// density, one per atom, see Dimaio et al.,
+/// Refinement of Protein Structures into Low-Resolution
+/// Density Maps Using Rosetta. Journal of Molecular
+/// Biology (2009) pp. 1-10
+/// ROSETTA_LOW_RESOLUTION guassian spheres in real space to represent
+/// density, one per residue. See reference above. Only
+/// useful at low resolution
+///
+/// The user can also choose if the density map should be cleared of its
+/// previous content before creating the density representation.
+///
+/// The first two types of density are generated in Fourier space. In order to
+/// avoid artifacts in the final density representation, the function avoids
+/// sharp frequency cutoffs by applying a Gaussian falloff. The user must
+/// provide the resolutions at which the cutoff should begin and end, as opposed
+/// to a single resolution cutoff value. For the other two density types, which
+/// are generated in real space, the average of the two falloff values is used
+/// as resolution parameter.
+///
+/// This function will only create a density represenation of the entities
+/// (or portion of entities ) that fall within the borders of the map.
+/// The user must take care that this condition is verified for all
+/// entities for which he wants a representation.
+///
+void DLLEXPORT_OST_MOL_ALG EntityToDensity(const mol::EntityView& entity_view,
+ img::MapHandle& map,
+ const DensityType& density_type,
+ Real falloff_start,
+ Real falloff_end,
+ bool clear_map_flag = false,
+ Real source_wavelength = 1.5418);
+
+}}} // ns
+
+#endif // OST_ENTITY_TO_DENSITY