corrected Example#2 and fixed a bug in the -discard-blobs argument

bin/filter_mrc/file_io.hpp

@@ -284,7 +284,6 @@ ReadBlobCoordsFile(string in_coords_file_name, //!< name of file we will read
                    vector<array<Coordinate, 3> > *pCrds=nullptr, //!< store the blob coordinates here (if !=nullptr)
                    vector<Scalar> *pDiameters=nullptr, //!< store the blob diameters here (if !=nullptr)
                    vector<Scalar> *pScores=nullptr, //!< store blob scores here (if !=nullptr)
-                   Scalar distance_scale=1.0, //!< divide all distances and coordinates by this value
                    Scalar diameter_override=-1.0, //!< use this diameter (useful if no 4th column is present)
                    Scalar score_default=0.0, //!< default "score" (if no 5th column is present)
                    Scalar diameter_factor=1.0, //!< multiply all diameters in the file by this number
@@ -325,22 +324,16 @@ ReadBlobCoordsFile(string in_coords_file_name, //!< name of file we will read
     ssLine >> x;
     ssLine >> y;
     ssLine >> z;
-    double ix, iy, iz;
-    ix = floor((x / distance_scale) + 0.5);
-    iy = floor((y / distance_scale) + 0.5);
-    iz = floor((z / distance_scale) + 0.5);
     array<Coordinate, 3> ixiyiz;
-    ixiyiz[0] = ix;
-    ixiyiz[1] = iy;
-    ixiyiz[2] = iz;
+    ixiyiz[0] = x;
+    ixiyiz[1] = y;
+    ixiyiz[2] = z;
 
     Scalar diameter = -1.0;
     Scalar score = score_default;
     if (ssLine) { // Does the file contain a 4th column? (the diameter)
       Scalar _diameter;
       ssLine >> _diameter;
-      // convert from physical distance to # of voxels:
-      _diameter /= distance_scale;
       if (ssLine)
         diameter = _diameter;
       custom_diameters = true;

bin/filter_mrc/filter_mrc.cpp

@@ -29,8 +29,8 @@ using namespace std;
 
 
 string g_program_name("filter_mrc");
-string g_version_string("0.29.15");
-string g_date_string("2021-8-18");
+string g_version_string("0.29.16");
+string g_date_string("2021-9-05");
 
 
 
@@ -115,7 +115,7 @@ int main(int argc, char **argv) {
       image_size[d] = tomo_in.header.nvoxels[d];
 
     // ---- Voxel width? ----
-    float voxel_width[3];
+    float voxel_width[3] = {-1.0 , -1.0, -1.0}; //initial impossible value
     DetermineVoxelWidth(settings, tomo_in, mask, voxel_width);
 
 
@@ -391,9 +391,13 @@ int main(int argc, char **argv) {
     cerr << "allocating space for new 3D image..." << endl;
     MrcSimple tomo_out = tomo_in; //this will take care of allocating the array
 
-    if ((voxel_width[0] <= 0.0) ||
-        (voxel_width[1] <= 0.0) ||
-        (voxel_width[2] <= 0.0))
+    if (((voxel_width[0] <= 0.0) ||
+         (voxel_width[1] <= 0.0) ||
+         (voxel_width[2] <= 0.0))
+        &&
+        (!
+         (settings.filter_type == settings.BLOB_NONMAX_SUPPRESSION) ||
+         (settings.filter_type == settings.BLOB_NONMAX_SUPERVISED_MULTI)))
       throw VisfdErr("Error in tomogram header: Invalid voxel width(s).\n"
                      "Use the -w argument to specify the voxel width.");
 
@@ -573,7 +577,6 @@ int main(int argc, char **argv) {
     }
 
 
-
     else if (settings.filter_type == settings.BLOB_NONMAX_SUPERVISED_MULTI) {
 
       // Use training data to choose the right threshold(s) for discarding blobs

bin/filter_mrc/handlers.cpp

@@ -445,11 +445,18 @@ HandleBlobsNonmaxSuppression(const Settings &settings,
                      &crds,
                      &diameters,
                      &scores,
-                     voxel_width[0],
                      settings.sphere_decals_diameter,
                      settings.sphere_decals_foreground,
                      settings.sphere_decals_scale);
 
+  if (voxel_width_ > 0.0) {
+    assert(crds.size() == diameters.size());
+    for (size_t i = 0; i < crds.size(); i++) {
+      diameters[i] /= voxel_width_;
+      for (int d=0; d < 3; d++)
+          crds[i][d] /= voxel_width_;
+    }
+  }
 
 
   cerr << " --- discarding blobs in file \n"
@@ -515,15 +522,29 @@ HandleBlobsNonmaxSuppression(const Settings &settings,
   }
 
   // Discard overlapping blobs?
-  cerr << "  discarding overlapping blobs" << endl;
-  DiscardOverlappingBlobs(crds,
-                          diameters, 
-                          scores,
-                          settings.nonmax_min_radial_separation_ratio,
-                          settings.nonmax_max_volume_overlap_large,
-                          settings.nonmax_max_volume_overlap_small,
-                          SORT_DECREASING_MAGNITUDE,
-                          &cerr);
+
+  if ((settings.nonmax_min_radial_separation_ratio > 0) ||
+      (settings.nonmax_max_volume_overlap_large !=
+       std::numeric_limits<float>::infinity()) ||
+      (settings.nonmax_max_volume_overlap_small !=
+       std::numeric_limits<float>::infinity())) {
+
+    if (voxel_width_ < 0)
+      throw VisfdErr("Error: Checking for overlapping blobs requires that you either specify the\n"
+                     "       voxel width (using the \"-w\" argument).  Alternatively, if the file\n"
+                     "       containing the original tomogram has accurate voxel width information,\n"
+                     "       you can specify that file instead (using the \"-in\" argument).\n");
+
+    cerr << "  discarding overlapping blobs" << endl;
+    DiscardOverlappingBlobs(crds,
+                            diameters, 
+                            scores,
+                            settings.nonmax_min_radial_separation_ratio,
+                            settings.nonmax_max_volume_overlap_large,
+                            settings.nonmax_max_volume_overlap_small,
+                            SORT_DECREASING_MAGNITUDE,
+                            &cerr);
+  }
 
   cerr << " " << crds.size() << " blobs remaining" << endl;
 
@@ -539,7 +560,6 @@ HandleBlobsNonmaxSuppression(const Settings &settings,
 
     cerr << "  discarding blobs based on score using training data" << endl;
 
-
     DiscardBlobsByScoreSupervised(crds,
                                   diameters,
                                   scores,
@@ -558,10 +578,10 @@ HandleBlobsNonmaxSuppression(const Settings &settings,
     fstream out_coords_file;
     out_coords_file.open(settings.out_coords_file_name, ios::out);
     for (size_t i=0; i < crds.size(); i++) {
-      out_coords_file << crds[i][0]*voxel_width[0] << " "
-                      << crds[i][1]*voxel_width[1] << " "
-                      << crds[i][2]*voxel_width[2] << " " 
-                      << diameters[i]*voxel_width[0] << " " 
+      out_coords_file << crds[i][0] << " "
+                      << crds[i][1] << " "
+                      << crds[i][2] << " " 
+                      << diameters[i] << " " 
                       << scores[i]
                       << endl;
     }
@@ -579,6 +599,14 @@ void
 HandleBlobScoreSupervisedMulti(const Settings &settings,
                                float voxel_width[3])
 {
+  // At some point I was trying to be as general as possible and allowed
+  // for the possibility that voxels need not be cubes (same width x,y,z)
+  // Now, I realize that allowing for this possibility would slow down some
+  // calculations considerably, so I just assume cube-shaped voxels:
+  float voxel_width_ = voxel_width[0];
+  assert((voxel_width[0] == voxel_width[1]) &&
+         (voxel_width[1] == voxel_width[2]));
+
   float threshold_lower_bound;
   float threshold_upper_bound;
 
@@ -588,12 +616,11 @@ HandleBlobScoreSupervisedMulti(const Settings &settings,
   vector<vector<float> > blob_diameters_multi(Nsets);
   vector<vector<float> > blob_scores_multi(Nsets);
 
-  for (int I = 0; I < settings.multi_in_coords_file_names.size(); ++I) {
+  for (int I = 0; I < Nsets; ++I) {
     ReadBlobCoordsFile(settings.multi_in_coords_file_names[I],
                        &blob_crds_multi[I],
                        &blob_diameters_multi[I],
                        &blob_scores_multi[I],
-                       voxel_width[0],
                        settings.sphere_decals_diameter,
                        settings.sphere_decals_foreground,
                        settings.sphere_decals_scale);
@@ -605,6 +632,16 @@ HandleBlobScoreSupervisedMulti(const Settings &settings,
   assert(Nsets == settings.multi_training_pos_crds.size());
   assert(Nsets == settings.multi_training_neg_crds.size());
 
+  if (voxel_width_ > 0.0) {
+    for (int I = 0; I < Nsets; ++I) {
+      for (size_t i = 0; i < blob_crds_multi[I].size(); i++) {
+        blob_diameters_multi[I][i] /= voxel_width_;
+        for (int d=0; d < 3; d++)
+          blob_crds_multi[I][i][d] /= voxel_width_;
+      }
+    }
+  }
+
   ChooseBlobScoreThresholdsMulti(blob_crds_multi,
                                  blob_diameters_multi,
                                  blob_scores_multi,
@@ -2333,7 +2370,17 @@ DetermineVoxelWidth(const Settings &settings,
     }
   }
   else {
-    // Otherwise, infer it from the header of the MRC file
+    if ((tomo_in.header.nvoxels[0] == 0) ||
+        (tomo_in.header.nvoxels[1] == 0) ||
+        (tomo_in.header.nvoxels[2] == 0)) {
+      // If the "tomo_in" image argument does not contain any voxels, then the
+      // user did not supply an input image.  Then there is no voxel width.
+      voxel_width[0] = -1.0; // Set the voxel width to an invalid number.
+      voxel_width[1] = -1.0; // This will let other parts of the code know
+      voxel_width[2] = -1.0; // that the voxel width is unknown.
+      return;
+    }
+    // Otherwise, infer the voxel width from the header metadata in the MRC file
     voxel_width[0] = tomo_in.header.cellA[0]/image_size[0];
     voxel_width[1] = tomo_in.header.cellA[1]/image_size[1];
     voxel_width[2] = tomo_in.header.cellA[2]/image_size[2];

doc/doc_filter_mrc.md

@@ -118,7 +118,8 @@ filter_mrc -in tomogram.rec \
 
 # Now discard the faint, noisy, or overlapping blobs.
 
-filter_mrc -discard-blobs tomogram_blobs.txt tomogram_ribosomes.txt \
+filter_mrc -in tomogram.rec \
+  -discard-blobs tomogram_blobs.txt tomogram_ribosomes.txt \
   -minima-threshold -70 \
   -radial-separation 0.8 \
   -mask cytoplasmic_volume.mrc   # <- optional

lib/visfd/feature.hpp

@@ -615,14 +615,27 @@ SortBlobs(vector<array<Scalar1,3> >& blob_crds,//!< x,y,z of each blob's center
 
 
 
-/// @brief  This variant of the function also takes care of discarding 
-///         training data which is not sufficiently close to any of the blobs.
-///         It also throws an exception if the remaining training data is empty.
+/// @brief Figure out the score for blobs located at each position in crds[].
+///        Terminology:  Every blob has a position, diameter, and a score.
+///        The j'th "blob" is a sphere, centerd at blob_crds[j], with diameter
+///        blob_diameters[j].  Associated with every blob is a "score"
+///        (a number) stored in blob_scores[j].  If the coordinates of the
+///        i'th data point (located in "in_training_crds[i]") lies within
+///        one of the spherical blobs, then store the corresponding coordinates,
+///        score, and the classification (accepted or rejected) for that blob
+///        in these vectors:
+///           out_training_crds[i],
+///           out_training_scores[i],
+///           out_training_accepted[i]
+///        (If crds[i] lies inside more than one blob, give priority to blobs
+///        occuring later in the list.)  If in_training_crds[i] lies outside
+///        any of the blobs, then discard this data point and do not append
+///        it to the output lists.)
 /// @return This function has no return value.
 ///         The results are stored in:
 ///           out_training_crds,
-///           out_training_accepted,
-///           out_training_scores
+///           out_training_scores,
+///           out_training_accepted
 
 template<typename Scalar>
 static void
@@ -647,7 +660,7 @@ FindBlobScores(const vector<array<Scalar,3> >& in_training_crds, //!< locations
   vector<size_t> in_training_which_blob; // which blob (sphere) contains this position?
   vector<Scalar> in_training_scores; // what is the score of that blob?
 
-  // The next function is defined in "feather_implementation.hpp"
+  // The next function is defined in "feature_implementation.hpp"
   _FindBlobScores(in_training_crds,
                   in_training_scores,
                   in_training_which_blob,

lib/visfd/feature_implementation.hpp

@@ -41,7 +41,7 @@ namespace visfd {
 ///        than one sphere, give priority spheres occuring later in the list.)
 ///        If crds[i] lies outside any of the spheres, store -infinity there.
 ///
-/// @returns void.  Results are stored in "scores".
+/// @returns void.  Results are stored in "scores" and "sphere_ids".
 ///
 /// @note  THIS FUNCTION WAS NOT INTENDED FOR PUBLIC USE