Moved Entity to Density to the main repository

git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2282 5a81b35b-ba03-0410-adc8-b2c5c5119f08

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})
 

modules/mol/alg/pymod/export_entity_to_density.cc

+//------------------------------------------------------------------------------
+// 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

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
 }

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})
 

modules/mol/alg/src/entity_to_density.cc

+#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
+
+

modules/mol/alg/src/entity_to_density.hh

+#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