Split EntityToDensity into two functions, fixed bug

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

modules/mol/alg/pymod/export_entity_to_density.cc

@@ -21,6 +21,9 @@
  * Author Juergen Haas
  */
 #include <boost/python.hpp>
+#include <ost/config.hh>
+
+#if OST_IMG_ENABLED
 
 #include <ost/mol/alg/entity_to_density.hh>
 
@@ -28,48 +31,19 @@ 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;
-}
-
+BOOST_PYTHON_FUNCTION_OVERLOADS(etd_rosetta, EntityToDensityRosetta, 4, 6)
+BOOST_PYTHON_FUNCTION_OVERLOADS(etd_scattering, EntityToDensityScattering, 4, 6)
 
 void export_entity_to_density()
 {
-
-    def("EntityToDensity",EntityToDensity);
-    def("EntityToDensity",etd1);
-    def("EntityToDensity",etd2);
+    def("EntityToDensityRosetta",EntityToDensityRosetta,etd_rosetta());
+    def("EntityToDensityScattering",EntityToDensityScattering,etd_scattering());
 
     enum_<DensityType>("DensityType")
-    .value("SCATTERING_FACTORS", SCATTERING_FACTORS)
-    .value("ROSETTA_HIGH_RESOLUTION", ROSETTA_HIGH_RESOLUTION)
-    .value("ROSETTA_LOW_RESOLUTION", ROSETTA_LOW_RESOLUTION)
+    .value("HIGH_RESOLUTION", HIGH_RESOLUTION)
+    .value("LOW_RESOLUTION", LOW_RESOLUTION)
     .export_values()
     ;
- }
+}
 #endif

modules/mol/alg/src/entity_to_density.cc

@@ -42,8 +42,10 @@ void FillAtomScatteringPropsTable(AtomScatteringPropsTable& table)
   boost::filesystem::path loc(fullpath);
   boost::filesystem::ifstream infile(loc);
   if (!infile) {
-    LOGN_ERROR("Couldn't find " << fullpath);
-    return;
+    String msg;
+    std::stringstream msgstr(msg);
+    msgstr << "Couldn't find " << fullpath << std::endl;
+    throw ost::Error(msg);
   }
   String line;
   std::getline(infile, line);
@@ -67,7 +69,6 @@ void FillAtomScatteringPropsTable(AtomScatteringPropsTable& table)
                >> atom_props.b4;
     table.push_back(atom_props);
   }
-  std::cout << table.size() << std::endl;
 };
 
 
@@ -91,14 +92,17 @@ 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),
+                             Real source_wavelength,
+                             geom::Vec3 map_start,
+                             geom::Vec3 map_end
+                             ):
+      entity_view_(entity_view),
       falloff_start_frequency_(1.0/falloff_start),
       falloff_end_frequency_(1.0/falloff_end),
-      source_wavelength_(source_wavelength)
+      source_wavelength_(source_wavelength),
+      map_start_(map_start),map_end_(map_end)
   {
     if (scattering_props_table_loaded_ == false)
     {
@@ -109,8 +113,7 @@ public:
 
   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();
 
@@ -133,15 +136,15 @@ public:
       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])
+          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
@@ -149,7 +152,7 @@ public:
             // Eventually, when maps and images will be merged,
             // it will substituted by the map's xyz2uvw method
 
-            geom::Vec3 adjusted_coord  = coord-map_start;
+            geom::Vec3 adjusted_coord  = coord-map_start_;
 
             for (img::ExtentIterator mp_it(reduced_extent);
                !mp_it.AtEnd();++mp_it)
@@ -217,6 +220,9 @@ private:
   static AtomScatteringPropsTable scatt_props_table_;
   static bool scattering_props_table_loaded_;
   Real source_wavelength_;
+  geom::Vec3 map_start_;
+  geom::Vec3 map_end_;
+
 };
 
 AtomScatteringPropsTable EntityToDensityHelperBase::scatt_props_table_ =
@@ -253,7 +259,7 @@ public:
     mol::EntityView effective_entity_view=entity_view_;
     Real k,C;
 
-    if (density_type_ == ROSETTA_HIGH_RESOLUTION)
+    if (density_type_ == HIGH_RESOLUTION)
     {
       k = (M_PI*M_PI)/(resolution_*resolution_);
       C = sqrt((M_PI*M_PI*M_PI)/(k*k*k));
@@ -292,9 +298,8 @@ public:
               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 adjusted_coord  = coord-map_start;
             geom::Vec3 pixel_coord=is.CoordToIndex(coord);
             img::Point rounded_pixel_coord(round(pixel_coord[0]),
                                             round(pixel_coord[1]),
@@ -366,53 +371,72 @@ typedef img::ImageStateConstModIPAlgorithm<EntityToDensityRosettaHelperBase>
 
 } // 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)
+void EntityToDensityScattering(const mol::EntityView& entity_view,
+                                     img::MapHandle& map,
+                                     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);
+ geom ::Vec3 rs_sampl = map.GetSpatialSampling();
+  geom ::Vec3 abs_orig = map.GetAbsoluteOrigin();
+  geom::Vec3 map_start = geom::Vec3(abs_orig[0]+map.GetExtent().GetStart()[0]*rs_sampl[0],
+                                    abs_orig[1]+map.GetExtent().GetStart()[1]*rs_sampl[1],
+                                    abs_orig[2]+map.GetExtent().GetStart()[2]*rs_sampl[2]);
+
+  geom::Vec3 map_end = geom::Vec3(abs_orig[0]+map.GetExtent().GetEnd()[0]*rs_sampl[0],
+                                  abs_orig[1]+map.GetExtent().GetEnd()[1]*rs_sampl[1],
+                                  abs_orig[2]+map.GetExtent().GetEnd()[2]*rs_sampl[2]);
+  detail::EntityToDensityHelper e_to_d_helper(entity_view,
+                                              falloff_start,
+                                              falloff_end,
+                                              source_wavelength, map_start,map_end);
+  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.SetSpatialSampling(map.GetSpatialSampling());
+    mm.SetAbsoluteOrigin(map.GetAbsoluteOrigin());
+    map.Swap(mm);
+  } else {
     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
+  }
+
+  map.ApplyIP(e_to_d_helper);
+  map.ApplyIP(img::alg::DFT());
+}
+
+
+void EntityToDensityRosetta(const mol::EntityView& entity_view,
+                            img::MapHandle& map,
+                            const DensityType& density_type,
+                            Real resolution,
+                            bool clear_map_flag,
+                            Real source_wavelength)
+
+{
+  if(resolution <=0.0) throw std::runtime_error("Invalid resolution");
+  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.SetSpatialSampling(map.GetSpatialSampling());
+    mm.SetAbsoluteOrigin(map.GetAbsoluteOrigin());
+    map.Swap(mm);
+  }
+
+  detail::EntityToDensityRosettaHelper e_to_d_r_helper
                                                 (entity_view,density_type,
                                                 resolution);
-    map.ApplyIP(e_to_d_r_helper);
-  }
+  map.ApplyIP(e_to_d_r_helper);
 }
 
+
 }}} // ns
 
 

modules/mol/alg/src/entity_to_density.hh

@@ -12,20 +12,43 @@ 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
+  HIGH_RESOLUTION,
+  LOW_RESOLUTION
 };
 
+/// \brief create a density representation of an entity in a density map (using electron scattering factors)
+///
+/// This functions creates a density representation of the entity provided by
+/// the user in a density map, in Fourier space using the correct scattering
+/// factors for the elements involved.
+///
+/// The user can also choose if the density map should be cleared of its
+/// previous content before creating the density representation.
+///
+/// The density is 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.
+///
+/// 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 EntityToDensityScattering(const mol::EntityView& entity_view,
+                                        img::MapHandle& map,
+                                        Real falloff_start,
+                                        Real falloff_end,
+                                        bool clear_map_flag = false,
+                                        Real source_wavelength = 1.5418);
+
 /// \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
@@ -38,27 +61,19 @@ enum DensityType
 /// 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.
+/// The user must also provide a 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,
+void DLLEXPORT_OST_MOL_ALG EntityToDensityRosetta(const mol::EntityView& entity_view,
                                         img::MapHandle& map,
                                         const DensityType& density_type,
-                                        Real falloff_start,
-                                        Real falloff_end,
+                                        Real resolution,
                                         bool clear_map_flag = false,
                                         Real source_wavelength = 1.5418);
-
 }}} // ns
 
 #endif // OST_ENTITY_TO_DENSITY