From 080f80b0e907654317bca3736f3b9d8c48499ddd Mon Sep 17 00:00:00 2001 From: Andrew Jewett Date: Mon, 6 Sep 2021 12:10:35 -0700 Subject: [PATCH] fixed more blob-detection bugs related to Example 2 --- bin/filter_mrc/filter_mrc.cpp | 4 ++-- bin/filter_mrc/handlers.cpp | 16 ++++++++++------ doc/doc_filter_mrc.md | 20 ++++++++++++++++---- tests/test_blob_detection.sh | 6 +++--- 4 files changed, 31 insertions(+), 15 deletions(-) diff --git a/bin/filter_mrc/filter_mrc.cpp b/bin/filter_mrc/filter_mrc.cpp index f52a98e..e33571f 100644 --- a/bin/filter_mrc/filter_mrc.cpp +++ b/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.16"); -string g_date_string("2021-9-05"); +string g_version_string("0.29.17"); +string g_date_string("2021-9-06"); diff --git a/bin/filter_mrc/handlers.cpp b/bin/filter_mrc/handlers.cpp index d9f963d..2143255 100644 --- a/bin/filter_mrc/handlers.cpp +++ b/bin/filter_mrc/handlers.cpp @@ -454,7 +454,7 @@ HandleBlobsNonmaxSuppression(const Settings &settings, 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_; + crds[i][d] = floor((crds[i][d] / voxel_width_) + 0.5); } } @@ -576,12 +576,15 @@ HandleBlobsNonmaxSuppression(const Settings &settings, // Write out the remaining blobs to a file: if (settings.out_coords_file_name != "") { fstream out_coords_file; + double voxel_width_positive = 1.0; + if (voxel_width_ > 0.0) + voxel_width_positive = voxel_width_; 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] << " " - << crds[i][1] << " " - << crds[i][2] << " " - << diameters[i] << " " + out_coords_file << crds[i][0] * voxel_width_positive << " " + << crds[i][1] * voxel_width_positive << " " + << crds[i][2] * voxel_width_positive << " " + << diameters[i] * voxel_width_positive << " " << scores[i] << endl; } @@ -637,7 +640,8 @@ HandleBlobScoreSupervisedMulti(const Settings &settings, 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_; + blob_crds_multi[I][i][d] = + floor((blob_crds_multi[I][i][d] / voxel_width_) + 0.5); } } } diff --git a/doc/doc_filter_mrc.md b/doc/doc_filter_mrc.md index 5edfdbc..b140c98 100644 --- a/doc/doc_filter_mrc.md +++ b/doc/doc_filter_mrc.md @@ -81,14 +81,22 @@ creates a new image file ("membranes.rec") with the membranes emphasized. View this file (eg. using IMOD) to see if the membranes are successfully detected. You may need to adjust this thickness parameter (and the other parameters) until the membrane is clearly visible. -**This is extremely slow, so try this on a small, cropped image first.** -(In my experience, 60-70 Angstroms is a reasonable default thickness parameter +In my experience, 60-70 Angstroms is a reasonable default thickess parameter to detect lipid bilayers. Most Cryo-EM tomogram files have voxel widths in units of Angstroms. You can override this by specifying -the physical width of each voxel using the [-w argument](#Voxel-Width).) +the physical width of each voxel using the [-w argument](#Voxel-Width). +**This operation is extremely slow, so try this on a small, +cropped image first.** +*(Note: Extra care must be taken if your image contains gold fiducial markers, +ice, or other extremely dark objects which can throw the detector off. +You can usually reduce the severity of this problem by including the +"-cl -0.4 0.4" argument to clip these extreme voxel brightnesses. +Alternatively you can remove gold fiducial markers from consideration +by finding their location beforehand using the "-blob" argument, +and runing "filter_mrc" again later using the "-mask-spheres" argument.)* Note: After the membrane features are detected, they must be analyzed. -This typically requires running the "filter_mrc" program multiple times. +This typically requires running the "filter_mrc" program again, multiple times. The optional ["-save-progress"](#-save-progress-FILENAME) argument used here allows us to skip the time consuming process of @@ -137,6 +145,10 @@ argument used in the second step is optional. I use it in this example because we are searching for ribosomes, so we want to restrict our search to blobs which lie *inside* the cell. +*(Note: Extra care must be taken if your image contains ice, gold fiducial +markers, or other extremely dark objects which can throw the detector off. +You can usually reduce the severity of this problem by including the +"-cl -0.4 0.4" argument to clip these extreme voxel brightnesses.)* Note: diff --git a/tests/test_blob_detection.sh b/tests/test_blob_detection.sh index 5cf381d..68a4f46 100755 --- a/tests/test_blob_detection.sh +++ b/tests/test_blob_detection.sh @@ -46,11 +46,11 @@ test_blob_detection() { # Test supervised learning of blob threshold parameter (single image) ../bin/filter_mrc/filter_mrc -w 19.6 -mask test_blob_detect_mask.rec -in test_blob_detect.rec -discard-blobs test_blobs.txt test_blobs_sep_${SEP}_SUPERVISED.txt -blob-separation ${SEP} -auto-thresh score -supervised test_supervised_pos.txt test_supervised_neg.txt >& test_log_e.txt - assertTrue "visualization failed. File test_blobs_sep_${SEP}_SUPERVISED.txt was not created" "[ -s test_blobs_sep_${SEP}_SUPERVISED.txt ]" + assertTrue "\"-supervised\" non-max suppression failed: File test_blobs_sep_${SEP}_SUPERVISED.txt was not created" "[ -s test_blobs_sep_${SEP}_SUPERVISED.txt ]" NBLOBS_SUPERVISED_SINGLE=`wc test_blobs_sep_${SEP}_SUPERVISED.txt | awk '{print $1}'` - assertTrue "supervised-learning failed: number of remaining blobs == 0" "[ $NBLOBS_SUPERVISED_SINGLE -gt 0 ]" + assertTrue "\"-supervised\" non-max suppression failed: Number of remaining blobs == 0" "[ $NBLOBS_SUPERVISED_SINGLE -gt 0 ]" THRESH_SUPERVISED_SINGLE=`grep 'threshold upper bound:' < test_log_e.txt | awk '{print $4}'` @@ -67,7 +67,7 @@ test_blob_detection() { echo "test_supervised_pos.txt test_supervised_neg.txt test_blobs_sep_${SEP}.txt" > test_supervised_multi.txt echo "test_supervised_pos.txt test_supervised_neg.txt test_blobs_sep_${SEP}.txt" >> test_supervised_multi.txt - + # run the learning algorithm to choose the threshold ../bin/filter_mrc/filter_mrc -w 19.6 -in test_blob_detect.rec -auto-thresh score -supervised-multi test_supervised_multi.txt >& test_log_e.txt