From 542e6f4160ba55aea5b6b20004e167a8ad2e17fb Mon Sep 17 00:00:00 2001
From: Markus Hollander <markus-hollander@web.de>
Date: Mon, 20 Jul 2020 18:59:14 +0200
Subject: [PATCH] added closeness to centre of mass as an additional axis end
 point

---
 src/grid.h    |  2 ++
 src/trace.cpp | 17 +++++++++++++++++
 2 files changed, 19 insertions(+)

diff --git a/src/grid.h b/src/grid.h
index df13cc8..2bcb2f4 100644
--- a/src/grid.h
+++ b/src/grid.h
@@ -440,6 +440,8 @@ struct PoreGrid : Grid {
     // HELP
     // small number to avoid divisions by 0
     const double perturb;
+    // true if the pore has a sufficiently large patch at the surface
+    bool has_surface_patch = false;
 
     // constructor
     PoreGrid(const PoreCluster &cluster, const double perturb) :
diff --git a/src/trace.cpp b/src/trace.cpp
index 089057d..b5cf202 100644
--- a/src/trace.cpp
+++ b/src/trace.cpp
@@ -156,6 +156,7 @@ void compute_starting_points(PoreGrid &grid, const size_t min_size) {
                 AxisPoint point = AxisPoint(grid.index(max_radius_vec), max_radius_vec,
                         "surface patch " + std::to_string(count));
                 grid.starting_points.push_back(point);
+                grid.has_surface_patch = true;
                 count++;
             }
         }
@@ -248,6 +249,10 @@ void trace(PoreGrid &grid, const size_t start_idx) {
     // initialise the minimum score as infinity
     double min_score = INFINITY;
     int min_index = -1;
+    // find the box on the outside of the pore/cavity with the best balance of closeness to the centre of mass and
+    // path length, that should hopefully go to and through the centre of the protein
+    double min_COM = INFINITY;
+    int COM_index = -1;
     // iterate over all boxes that are part of the pore/cavity and lay on its outer layer. make sure they are not
     // the start point and do not have a score of infinity
     for (size_t i = 0; i < grid.size(); i++) {
@@ -263,6 +268,13 @@ void trace(PoreGrid &grid, const size_t start_idx) {
             min_score = adjusted_score;
             min_index = i;
         }
+        // divide the distance to the centre of mass by the length of the path to identify the best balance
+        double COM_score = grid.distance_to_centre[i] / lengths[i];
+        // if the score is lower than the current minimum, set the bo as the candidate for the end point
+        if (COM_score < min_COM) {
+            min_COM = COM_score;
+            COM_index = i;
+        }
     }
 
     // add all not yet processed start points from surface patches to the list of possible axis end points
@@ -273,6 +285,11 @@ void trace(PoreGrid &grid, const size_t start_idx) {
         AxisPoint end = AxisPoint(min_index, grid.coordinates[min_index], "lowest score point");
         end_points.push_back(end);
     }
+    // add the end point with the best balance of length and closeness to centre of mass
+    if (COM_index > -1) {
+        AxisPoint end = AxisPoint(COM_index, grid.coordinates[COM_index], "centre of mass point");
+        end_points.push_back(end);
+    }
 
     // compute the axes between the current start point and all identified possible end points
     for (const AxisPoint &end: end_points) {