diff --git a/bin/filter_mrc/filter_mrc.cpp b/bin/filter_mrc/filter_mrc.cpp index f8833d0..49ef274 100644 --- a/bin/filter_mrc/filter_mrc.cpp +++ b/bin/filter_mrc/filter_mrc.cpp @@ -29,7 +29,7 @@ using namespace std; string g_program_name("filter_mrc"); -string g_version_string("0.28.1"); +string g_version_string("0.28.2"); string g_date_string("2021-7-11"); @@ -367,76 +367,118 @@ int main(int argc, char **argv) { "Use the -w argument to specify the voxel width."); - // ---- filtering ---- - // perform an operation which generates a new image based on the old image + + // ---- Select the primary operation we will perform on the image ---- + if (settings.filter_type == Settings::NONE) { + cerr << "filter_type = Intensity Map \n"; - // Not needed: - //for (int iz = 0; iz < size[2]; iz++) - // for (int iy = 0; iy < size[1]; iy++) - // for (int ix = 0; ix < size[0]; ix++) - // tomo_out.aaafI[iz][iy][ix]=tomo_in.aaafI[iz][iy][ix]; - // (We have copied the contents from tomo_in into tomo_out already.) + } + else if (settings.filter_type == Settings::DILATION) { + + // Apply a greyscale dilation filter to the image. + HandleDilation(settings, tomo_in, tomo_out, mask, voxel_width); + + } + + + + else if (settings.filter_type == Settings::EROSION) { + + // Apply a greyscale erosion filter to the image. + HandleErosion(settings, tomo_in, tomo_out, mask, voxel_width); + + } + + + + else if (settings.filter_type == Settings::OPENING) { + + // Apply a greyscale opening filter to the image. + HandleOpening(settings, tomo_in, tomo_out, mask, voxel_width); + + } + + + + else if (settings.filter_type == Settings::CLOSING) { + + // Apply a greyscale closing filter to the image. + HandleClosing(settings, tomo_in, tomo_out, mask, voxel_width); + + } + + + else if (settings.filter_type == Settings::GAUSS) { + // Apply a Gaussian filter to the image. HandleGauss(settings, tomo_in, tomo_out, mask, voxel_width); - - } // if (settings.filter_type == Settings::GAUSS) + + } else if (settings.filter_type == Settings::GGAUSS) { + // Apply a generalized Gaussian filter HandleGGauss(settings, tomo_in, tomo_out, mask, voxel_width); - } //else if (settings.filter_type == Settings::GGAUSS) + } else if (settings.filter_type == Settings::DOG) { + // Apply a Difference-of-Gaussians filter HandleDog(settings, tomo_in, tomo_out, mask, voxel_width); - } //if (settings.filter_type == Settings::DOG) + } else if (settings.filter_type == Settings::DOGG) { + // Apply a generalized DoG filter. + // (Note: This kind of operation is probably not useful. + // I may remove this feature in the future. -Andrew 2021-7-11) HandleDogg(settings, tomo_in, tomo_out, mask, voxel_width); - } //if (settings.filter_type == Settings::DOGG) + } #ifndef DISABLE_DOGGXY else if (settings.filter_type == Settings::DOGGXY) { + // Apply a generalized DoG filter in the XY direction + // and a Gaussian filter in the Z direction. HandleDoggXY(settings, tomo_in, tomo_out, mask, voxel_width); - } //else if (settings.filter_type = Settings::DOGGXY) + } #endif + else if (settings.filter_type == Settings::LOG_DOG) { + // Apply a LoG filter (implemented using the DoG approximation). HandleLoGDoG(settings, tomo_in, tomo_out, mask, voxel_width); - } //if (settings.filter_type == Settings::DOG_SCALE_FREE) - + } else if (settings.filter_type == Settings::BLOB) { HandleBlobDetector(settings, tomo_in, tomo_out, mask, voxel_width); - } //if (settings.filter_type == Settings::DOG) + } @@ -444,116 +486,140 @@ int main(int argc, char **argv) { (settings.filter_type == Settings::SURFACE_RIDGE) || (settings.filter_type == Settings::CURVE)) { - // find surface ridges (ie membranes or wide tubes) + // Detect 2-D surfaces or 1-D curves HandleTV(settings, tomo_in, tomo_out, mask, voxel_width); } + else if (settings.filter_type == Settings::WATERSHED) { - // perform watershed segmentation + // Perform watershed segmentation HandleWatershed(settings, tomo_in, tomo_out, mask, voxel_width); } + else if (settings.filter_type == Settings::CLUSTER_CONNECTED) { - // cluster adjacent nearby voxels into disconnected "islands" - // (this is similar to watershed segmentation) + // Cluster adjacent nearby voxels with high "saliency" into "islands" + // neighboring voxels (this is similar to watershed segmentation). HandleClusterConnected(settings, tomo_in, tomo_out, mask, voxel_width); } + else if (settings.filter_type == Settings::LOCAL_FLUCTUATIONS) { + // Apply a filter that measures the amount of brightness fluctuations + // in the local neighborhood. HandleLocalFluctuations(settings, tomo_in, tomo_out, mask, voxel_width); - } //else if (settings.filter_type == Settings::TEMPLATE_GGAUSS) + } + // ----- distance_filter ----- + + else if (settings.filter_type == settings.DISTANCE_TO_POINTS) { + // Apply a filter which replaces the voxel brightness with the distance + // to the nearest poing in a (user-supplied) point cloud. + HandleDistanceToPoints(settings, tomo_in, tomo_out, mask, voxel_width); - // ----- template matching with error reporting (probably not useful) ----- + } - else if (settings.filter_type == Settings::TEMPLATE_GGAUSS) { + // ----- sphere_decals_filter ----- - #ifndef DISABLE_TEMPLATE_MATCHING - HandleTemplateGGauss(settings, tomo_in, tomo_out, mask, voxel_width); - #endif //#ifndef DISABLE_TEMPLATE_MATCHING + else if (settings.filter_type == settings.SPHERE_DECALS) { - } //else if (settings.filter_type == Settings::TEMPLATE_GGAUSS) + // Draw a new image or annotate (draw on top of) an existing image. + HandleDrawSpheres(settings, tomo_in, tomo_out, mask, voxel_width); + } - else if (settings.filter_type == Settings::TEMPLATE_GAUSS) { + else if (settings.filter_type == settings.SPHERE_NONMAX_SUPPRESSION) { - #ifndef DISABLE_TEMPLATE_MATCHING - HandleTemplateGauss(settings, tomo_in, tomo_out, mask, voxel_width); - #endif //#ifndef DISABLE_TEMPLATE_MATCHING + vector > crds; + vector diameters; + vector scores; - } //else if (settings.filter_type == Settings::TEMPLATE_GAUSS) + // Discard overlapping or poor scoring blobs from a (user-supplied) list. + // (Note: In this case, we are not analyzing the image. + // The list of blobs is supplied by the user, most likely by + // running this program at an earlier time using the "-blob" + // argument.) + HandleBlobsNonmaxSuppression(settings, + mask, + voxel_width, + crds, + diameters, + scores); + } - else if (settings.filter_type == Settings::BLOB_RADIAL_INTENSITY) { - #ifndef DISABLE_INTENSITY_PROFILES - HandleBlobRadialIntensity(settings, tomo_in, tomo_out, mask, voxel_width); - #endif //#ifndef DISABLE_TEMPLATE_MATCHING + else if (settings.filter_type == settings.SPHERE_NONMAX_SUPERVISED_MULTI) { - } //else if (settings.filter_type == Settings::TEMPLATE_GAUSS) + // Use training data to choose the right threshold(s) for discarding blobs + HandleBlobScoreSupervisedMulti(settings, + voxel_width); + } - else if (settings.filter_type == Settings::BOOTSTRAP_DOGG) { + // ------- DEPRECIATED FEATURES ------- - #ifdef DISABLE_BOOSTRAPPING - HandleBootstrappDogg(settings, tomo_in, tomo_out, mask, voxel_width); - #endif //#ifndef DISABLE_BOOTSTRAPPING - } //else if (settings.filter_type == Settings::BOOTSTRAP_DOGG) { + else if (settings.filter_type == Settings::TEMPLATE_GGAUSS) { + #ifndef DISABLE_TEMPLATE_MATCHING + // Spherical template matching with RMS error reporting. + // (Not very useful. I may remove this feature in the future -A 2021-7-11) + HandleTemplateGGauss(settings, tomo_in, tomo_out, mask, voxel_width); + #endif //#ifndef DISABLE_TEMPLATE_MATCHING - // ----- distance_filter ----- + } - else if (settings.filter_type == settings.DISTANCE_TO_POINTS) { - HandleDistanceToPoints(settings, tomo_in, tomo_out, mask, voxel_width); + else if (settings.filter_type == Settings::TEMPLATE_GAUSS) { + + #ifndef DISABLE_TEMPLATE_MATCHING + // Spherical template matching with RMS error reporting. + // (Not very useful. I may remove this feature in the future -A 2021-7-11) + HandleTemplateGauss(settings, tomo_in, tomo_out, mask, voxel_width); + #endif //#ifndef DISABLE_TEMPLATE_MATCHING } - // ----- sphere_decals_filter ----- - else if (settings.filter_type == settings.SPHERE_DECALS) { + else if (settings.filter_type == Settings::BOOTSTRAP_DOGG) { - HandleDrawSpheres(settings, tomo_in, tomo_out, mask, voxel_width); + #ifdef DISABLE_BOOSTRAPPING + // (I may remove this feature in the future -Andrew 2021-7-11) + HandleBootstrappDogg(settings, tomo_in, tomo_out, mask, voxel_width); + #endif //#ifndef DISABLE_BOOTSTRAPPING } - else if (settings.filter_type == settings.SPHERE_NONMAX_SUPPRESSION) { - vector > crds; - vector diameters; - vector scores; + else if (settings.filter_type == Settings::BLOB_RADIAL_INTENSITY) { - HandleBlobsNonmaxSuppression(settings, - mask, - voxel_width, - crds, - diameters, - scores); + #ifndef DISABLE_INTENSITY_PROFILES + // This is a feature that a user requested for a specific project. + // It is probably not relevant to most users. + // (I may remove this feature in the future -Andrew 2021-7-11) + HandleBlobRadialIntensity(settings, tomo_in, tomo_out, mask, voxel_width); + #endif //#ifndef DISABLE_TEMPLATE_MATCHING } - else if (settings.filter_type == settings.SPHERE_NONMAX_SUPERVISED_MULTI) { - HandleBlobScoreSupervisedMulti(settings, - voxel_width); - } else { assert(false); //should be one of the choices above diff --git a/bin/filter_mrc/handlers.cpp b/bin/filter_mrc/handlers.cpp index ea529e0..003c83c 100644 --- a/bin/filter_mrc/handlers.cpp +++ b/bin/filter_mrc/handlers.cpp @@ -37,6 +37,75 @@ using namespace visfd; + +void +HandleDilation(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]) +{ + DilateSphere(settings.thickness_morphology, + tomo_in.header.nvoxels, + tomo_in.aaafI, + tomo_out.aaafI, + mask.aaafI, + false, + &cerr); +} + + +void +HandleErosion(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]) +{ + ErodeSphere(settings.thickness_morphology, + tomo_in.header.nvoxels, + tomo_in.aaafI, + tomo_out.aaafI, + mask.aaafI, + false, + &cerr); +} + + +void +HandleOpening(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]) +{ + OpenSphere(settings.thickness_morphology, + tomo_in.header.nvoxels, + tomo_in.aaafI, + tomo_out.aaafI, + mask.aaafI, + false, + &cerr); +} + + +void +HandleClosing(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]) +{ + CloseSphere(settings.thickness_morphology, + tomo_in.header.nvoxels, + tomo_in.aaafI, + tomo_out.aaafI, + mask.aaafI, + false, + &cerr); +} + + void HandleGGauss(const Settings &settings, MrcSimple &tomo_in, @@ -1057,8 +1126,6 @@ HandleLocalFluctuations(const Settings &settings, settings.filter_truncate_threshold, settings.normalize_near_boundaries, &cerr); - - tomo_out.FindMinMaxMean(); } @@ -1208,8 +1275,6 @@ HandleClusterConnected(const Settings &settings, - - void HandleTV(const Settings &settings, MrcSimple &tomo_in, diff --git a/bin/filter_mrc/handlers.hpp b/bin/filter_mrc/handlers.hpp index e1b0c27..f9a698d 100644 --- a/bin/filter_mrc/handlers.hpp +++ b/bin/filter_mrc/handlers.hpp @@ -3,6 +3,37 @@ #include "settings.hpp" +void +HandleDilation(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]); + + +void +HandleErosion(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]); + + +void +HandleOpening(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]); + + +void +HandleClosing(const Settings &settings, + MrcSimple &tomo_in, + MrcSimple &tomo_out, + MrcSimple &mask, + float voxel_width[3]); + void HandleGGauss(const Settings &settings, @@ -12,8 +43,6 @@ HandleGGauss(const Settings &settings, float voxel_width[3]); - - void HandleGauss(const Settings &settings, MrcSimple &tomo_in, @@ -22,8 +51,6 @@ HandleGauss(const Settings &settings, float voxel_width[3]); - - void HandleDogg(const Settings &settings, MrcSimple &tomo_in, @@ -32,8 +59,6 @@ HandleDogg(const Settings &settings, float voxel_width[3]); - - void HandleDog(const Settings &settings, MrcSimple &tomo_in, @@ -42,8 +67,6 @@ HandleDog(const Settings &settings, float voxel_width[3]); - - void HandleLoGDoG(const Settings &settings, MrcSimple &tomo_in, @@ -52,7 +75,6 @@ HandleLoGDoG(const Settings &settings, float voxel_width[3]); - void HandleBlobsNonmaxSuppression(const Settings &settings, MrcSimple &mask, @@ -62,8 +84,6 @@ HandleBlobsNonmaxSuppression(const Settings &settings, vector& scores); - - void HandleBlobScoreSupervisedMulti(const Settings &settings, float voxel_width[3]); @@ -78,7 +98,6 @@ HandleDrawSpheres(const Settings &settings, float voxel_width[3]); - void HandleBlobDetector(const Settings &settings, MrcSimple &tomo_in, @@ -87,7 +106,6 @@ HandleBlobDetector(const Settings &settings, float voxel_width[3]); - void HandleThresholds(const Settings &settings, MrcSimple &tomo_in, @@ -105,7 +123,6 @@ HandleExtrema(const Settings &settings, float voxel_width[3]); - void HandleLocalFluctuations(const Settings &settings, MrcSimple &tomo_in, diff --git a/bin/filter_mrc/settings.cpp b/bin/filter_mrc/settings.cpp index 353b4e4..1efec11 100644 --- a/bin/filter_mrc/settings.cpp +++ b/bin/filter_mrc/settings.cpp @@ -260,6 +260,7 @@ Settings::ParseArgs(vector& vArgs) } // if ((vArgs[i] == "-in") || (vArgs[i] == "-i")) + else if (vArgs[i] == "-image-size") { try { @@ -280,6 +281,7 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-image-size") + else if ((vArgs[i] == "-out") || (vArgs[i] == "-o")) { if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) @@ -290,6 +292,7 @@ Settings::ParseArgs(vector& vArgs) } // if ((vArgs[i] == "-out") || (vArgs[i] == "-o")) + else if ((vArgs[i] == "-outf") || (vArgs[i] == "-out-force")) { if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) @@ -301,6 +304,7 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-outf") + else if (vArgs[i] == "-mask") { if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) @@ -345,6 +349,7 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-mask-select") + else if ((vArgs[i] == "-mask-rect") || (vArgs[i] == "-mask-rectangle")) { @@ -376,6 +381,7 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-mask-rect") + else if ((vArgs[i] == "-mask-rect-subtract") || (vArgs[i] == "-mask-rectangle-subtract")) { @@ -486,6 +492,7 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-mask-rect-units") + else if (vArgs[i] == "-mask-out") { try { @@ -502,6 +509,7 @@ Settings::ParseArgs(vector& vArgs) } // if (vArgs[i] == "-mask-out") + else if (vArgs[i] == "-w") { try { @@ -544,7 +552,47 @@ Settings::ParseArgs(vector& vArgs) - else if ((vArgs[i] == "-dilate-gauss") || + else if ((vArgs[i] == "-dilation") || + (vArgs[i] == "-dilate")) + { + try { + if ((i+1 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) + throw invalid_argument(""); + thickness_morphology = stof(vArgs[i+1]); + filter_type = DILATION; + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by a positive number (\"s\"),\n"); + } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-dilation") + + + + else if ((vArgs[i] == "-erosion") || + (vArgs[i] == "-erode")) + { + try { + if ((i+1 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) + throw invalid_argument(""); + thickness_morphology = stof(vArgs[i+1]); + filter_type = EROSION; + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by a positive number (\"s\"),\n"); + } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-erosion") + + + + else if ((vArgs[i] == "-dilation-gauss") || + (vArgs[i] == "-dilate-gauss") || + (vArgs[i] == "-erosion-gauss") || (vArgs[i] == "-erode-gauss")) { float blur_distance; @@ -574,761 +622,568 @@ Settings::ParseArgs(vector& vArgs) - - - else if ((vArgs[i] == "-gauss-aniso") || (vArgs[i] == "-ggauss-aniso")) + else if ((vArgs[i] == "-opening") || + (vArgs[i] == "-open")) { try { - if ((i+3 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || - (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+3] == "") || (vArgs[i+3][0] == '-')) + if ((i+1 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) throw invalid_argument(""); - width_a[0] = stof(vArgs[i+1]); - width_a[1] = stof(vArgs[i+2]); - width_a[2] = stof(vArgs[i+3]); - width_b[0] = -1.0; - width_b[1] = -1.0; - width_b[2] = -1.0; - filter_type = GAUSS; - if (vArgs[i] == "-ggauss-aniso") - filter_type = GGAUSS; + thickness_morphology = stof(vArgs[i+1]); + filter_type = OPENING; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 3 positive numbers:\n" - " s_x s_y s_z\n" - " the Gaussian widths in the X, Y, and Z direction.)\n"); + " argument must be followed by a positive number (\"s\"),\n"); } - num_arguments_deleted = 4; - } //if (vArgs[i] == "-gauss-aniso") + num_arguments_deleted = 2; + } //if (vArgs[i] == "-opening") - else if ((vArgs[i] == "-gauss") || (vArgs[i] == "-ggauss")) + else if ((vArgs[i] == "-closing") || + (vArgs[i] == "-close")) { try { if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) throw invalid_argument(""); - width_a[0] = stof(vArgs[i+1]); - width_a[1] = width_a[0]; - width_a[2] = width_a[0]; - width_b[0] = -1.0; - width_b[1] = -1.0; - width_b[2] = -1.0; - filter_type = GAUSS; - if (vArgs[i] == "-ggauss") - filter_type = GGAUSS; + thickness_morphology = stof(vArgs[i+1]); + filter_type = CLOSING; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a positive number (\"s\"),\n" - " the Gaussian width\n"); + " argument must be followed by a positive number (\"s\"),\n"); } num_arguments_deleted = 2; - } //if (vArgs[i] == "-gauss") + } //if (vArgs[i] == "-closing") - else if ((vArgs[i] == "-dog-aniso") || - (vArgs[i] == "-dogg-aniso")) + else if (vArgs[i] == "-truncate") { try { - if ((i+6 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || - (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+3] == "") || (vArgs[i+3][0] == '-') || - (vArgs[i+4] == "") || (vArgs[i+4][0] == '-') || - (vArgs[i+5] == "") || (vArgs[i+5][0] == '-') || - (vArgs[i+6] == "") || (vArgs[i+6][0] == '-')) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) throw invalid_argument(""); - width_a[0] = stof(vArgs[i+1]); - width_a[1] = stof(vArgs[i+2]); - width_a[2] = stof(vArgs[i+3]); - width_b[0] = stof(vArgs[i+4]); - width_b[1] = stof(vArgs[i+5]); - width_b[2] = stof(vArgs[i+6]); - //if (width_b[0] <= width_a[0]) - // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" - // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); - filter_type = DOG; - if (vArgs[i] == "-dogg-aniso") - filter_type = DOGG; + filter_truncate_ratio = stof(vArgs[i+1]); + filter_truncate_threshold=-1.0; //(disables)override any filter_truncate_threshold settings + //filter_truncate_ratio_exp = exp(-pow(filter_truncate_ratio, + // n_exp)); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 6 positive numbers.\n"); + " argument must be followed by a number.\n"); } - num_arguments_deleted = 7; - } //if (vArgs[i] == "-dog-aniso") + num_arguments_deleted = 2; + } // if (vArgs[i] == "-truncate") - else if ((vArgs[i] == "-dog") || - (vArgs[i] == "-dogg")) + else if (vArgs[i] == "-truncate-threshold") { try { - if ((i+2 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || - (vArgs[i+2] == "") || (vArgs[i+2][0] == '-')) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) throw invalid_argument(""); - width_a[0] = stof(vArgs[i+1]); - width_a[1] = width_a[0]; - width_a[2] = width_a[0]; - width_b[0] = stof(vArgs[i+2]); - width_b[1] = width_b[0]; - width_b[2] = width_b[0]; - //if (width_b[0] <= width_a[0]) - // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" - // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); - filter_type = DOG; - if (vArgs[i] == "-dogg") - filter_type = DOGG; + filter_truncate_threshold = stof(vArgs[i+1]); + filter_truncate_ratio = -1.0; //(disables) override any filter_truncate_ratio settings } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 2 positive numbers.\n"); + " argument must be followed by a number.\n"); } - num_arguments_deleted = 3; - } //if (vArgs[i] == "-dog") + num_arguments_deleted = 2; + } // if (vArgs[i] == "-truncate-thresold") - #ifndef DISABLE_DOGGXY - else if (vArgs[i] == "-doggxy-aniso") + + else if (vArgs[i] == "-rescale") { try { - if ((i+5 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || - (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+3] == "") || (vArgs[i+3][0] == '-') || - (vArgs[i+4] == "") || (vArgs[i+4][0] == '-') || - (vArgs[i+5] == "") || (vArgs[i+5][0] == '-')) + if ((i+2 >= vArgs.size()) || + (vArgs[i+1] == "") || + (vArgs[i+2] == "")) throw invalid_argument(""); - width_a[0] = stof(vArgs[i+1]); - width_a[1] = stof(vArgs[i+2]); - width_b[0] = stof(vArgs[i+3]); - width_b[1] = stof(vArgs[i+4]); - //The "-doggxy" filter is a Difference-of-Generalized-Gaussians - //in the X,Y directions, multiplied by an ordinary Gaussian in - //the Z direction. - width_a[2] = stof(vArgs[i+5]); - width_b[2] = -1.0; //(This disables the second (negative) Gaussian) - //if (width_b[0] <= width_a[0]) - // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" - // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); - filter_type = DOGGXY; + use_intensity_map = true; + use_rescale_multiply = true; + out_rescale_multiply = stof(vArgs[i+1]); + out_rescale_offset = stof(vArgs[i+2]); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 5 positive numbers:\n" - " a_x a_y b_x b_y a_z\n" - " (I.E., the \"A\" and \"B\" Gaussian widths in X and Y directions,\n" - " followed by the Gaussian width in the Z direction.)\n"); + " argument must be followed by 2 numbers:\n" + " outA outB\n" + " (the desired minimum and maximum voxel intensity values for the final image)\n"); } - num_arguments_deleted = 6; - } //if (vArgs[i] == "-doggxy-aniso") + num_arguments_deleted = 3; + } // if (vArgs[i] == "-rescale") - else if (vArgs[i] == "-doggxy") + + else if (vArgs[i] == "-fill") { try { - if ((i+5 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || - (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+3] == "") || (vArgs[i+3][0] == '-')) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); - width_a[0] = stof(vArgs[i+1]); - width_a[1] = width_a[0]; - width_b[0] = stof(vArgs[i+2]); - width_b[1] = width_b[0]; - //The "-doggxy" filter is a Difference-of-Generalized-Gaussians - //in the X,Y directions, multiplied by an ordinary Gaussian in - //the Z direction. - width_a[2] = stof(vArgs[i+3]); - width_b[2] = -1.0; //(This disables the second (negative) Gaussian) - //if (width_b[0] <= width_a[0]) - // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" - // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); - filter_type = DOGGXY; + use_intensity_map = true; + use_rescale_multiply = true; + out_rescale_multiply = 0.0; + out_rescale_offset = stof(vArgs[i+1]); } catch (invalid_argument& exc) { - throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 3 positive numbers:\n" - " a_xy b_xy a_z\n" - " (I.E., the \"A\" and \"-B\" Gaussian widths in the XY plane,\n" - " followed by the Gaussian width in the Z direction.)\n"); + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by a number."); } - num_arguments_deleted = 4; - } //if (vArgs[i] == "-doggxy") - #endif //#ifndef DISABLE_DOGGXY + num_arguments_deleted = 2; + } // if (vArgs[i] == "-rescale") - else if ((vArgs[i] == "-log") || - (vArgs[i] == "-log-r") || - (vArgs[i] == "-log-d")) + else if ((vArgs[i] == "-thresh-range") || + (vArgs[i] == "-thresh-range-out")) { try { - if ((i+1 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) + if ((i+2 >= vArgs.size()) || + (vArgs[i+1] == "") || + (vArgs[i+2] == "")) throw invalid_argument(""); - float blob_sigma_multiplier = 1.0; - blob_sigma_multiplier = 1.0; - if (vArgs[i] == "-log") - blob_sigma_multiplier = 1.0; - if (vArgs[i] == "-log-r") - // (for a solid uniform 3-D sphere) - blob_sigma_multiplier = 1.0/sqrt(3.0); - if (vArgs[i] == "-log-d") - // (for a solid uniform 3-D sphere) - blob_sigma_multiplier = 1.0/(2.0*sqrt(3.0)); - log_width[0] = stof(vArgs[i+1]) * blob_sigma_multiplier; - log_width[1] = log_width[0]; - log_width[2] = log_width[0]; - m_exp = 2.0; - n_exp = 2.0; - filter_type = LOG_DOG; + out_thresh_a_value = stof(vArgs[i+1]); + out_thresh_b_value = stof(vArgs[i+2]); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a positive number:\n" - " the (approximate) width of the objects of interest.\n"); + " argument must be followed by 2 numbers:\n" + " outA outB\n" + " (the desired minimum and maximum voxel intensity values for the final image)\n"); } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-log") + num_arguments_deleted = 3; + } // if (vArgs[i] == "-outab") - else if (vArgs[i] == "-log-aniso") + + else if (vArgs[i] == "-rescale-min-max") { try { - if ((i+3 >= vArgs.size()) || - (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || - (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+3] == "") || (vArgs[i+3][0] == '-')) + if ((i+2 >= vArgs.size()) || + (vArgs[i+1] == "") || + (vArgs[i+2] == "")) throw invalid_argument(""); - log_width[0] = stof(vArgs[i+1]); - log_width[1] = stof(vArgs[i+2]); - log_width[2] = stof(vArgs[i+3]); - m_exp = 2.0; - n_exp = 2.0; - filter_type = LOG_DOG; + rescale_min_max_out = true; + out_rescale_max = stof(vArgs[i+1]); + out_rescale_min = stof(vArgs[i+2]); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 3 positive numbers:\n" - " the approximate widths of the objects of interest in the XYZ directions\n" - " (Due to the missing-wedge artifact they might differ.\n"); + " argument must be followed by 2 numbers:\n" + " outA outB\n" + " (the desired minimum and maximum voxel intensity values for the final image)\n"); } - num_arguments_deleted = 4; - } //if (vArgs[i] == "-log-aniso") + num_arguments_deleted = 3; + } // if (vArgs[i] == "-rescale-min-max") - else if ((vArgs[i] == "-spheres-nonmax-overlap") || - (vArgs[i] == "-max-volume-overlap") || - (vArgs[i] == "-max-overlap")) + else if ((vArgs[i] == "-no-rescale") || (vArgs[i] == "-norescale")) + { + rescale_min_max_out = false; + in_threshold_01_a = 1.0; + in_threshold_01_b = 1.0; + num_arguments_deleted = 1; + } // if (vArgs[i] == "-norescale") + + + + else if ((vArgs[i] == "-invert") || + (vArgs[i] == "-inv")) { + invert_output = true; + num_arguments_deleted = 1; + } // if (vArgs[i] == "-invert") + + + + else if ((vArgs[i] == "-thresh") || + (vArgs[i] == "-thresh-out")) { try { if (i+1 >= vArgs.size()) throw invalid_argument(""); - nonmax_max_volume_overlap_large = stof(vArgs[i+1]); - nonmax_min_radial_separation_ratio = 0.0; + use_intensity_map = true; + use_dual_thresholds = false; + in_threshold_01_a = stof(vArgs[i+1]); + in_threshold_01_b = in_threshold_01_a; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number:\n" - " the maximum overlap allowed between a pair of detected blobs\n" - " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + " argument must be followed by 1 number.\n"); } num_arguments_deleted = 2; - } //if (vArgs[i] == "-spheres-nonmax-overlap") + } - else if ((vArgs[i] == "-spheres-nonmax-overlap-small") || - (vArgs[i] == "-max-volume-overlap-small") || - (vArgs[i] == "-max-overlap-small")) - { + else if ((vArgs[i] == "-thresh2") || + (vArgs[i] == "-thresh2-out")) { try { - if (i+1 >= vArgs.size()) + if (i+2 >= vArgs.size()) throw invalid_argument(""); - nonmax_max_volume_overlap_small = stof(vArgs[i+1]); - nonmax_min_radial_separation_ratio = 0.0; + use_intensity_map = true; + use_dual_thresholds = false; + in_threshold_01_a = stof(vArgs[i+1]); + in_threshold_01_b = stof(vArgs[i+2]); + out_thresh2_use_clipping = false; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number:\n" - " the maximum overlap allowed between a pair of detected blobs\n" - " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + " argument must be followed by 2 numbers.\n"); } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-spheres-nonmax-overlap-small") + num_arguments_deleted = 3; + } - else if ((vArgs[i] == "-spheres-nonmax-overlap-radial") || - (vArgs[i] == "-max-overlap-radial")) - { + else if ((vArgs[i] == "-clip") || + (vArgs[i] == "-cl")) { try { - if (i+1 >= vArgs.size()) + if (i+2 >= vArgs.size()) throw invalid_argument(""); - float max_overlap = stof(vArgs[i+1]); - nonmax_min_radial_separation_ratio = (1.0 - max_overlap); - //REMOVE THIS CRUFT: - //nonmax_min_radial_separation_ratio = (1.0 - max_overlap) * sqrt(3.0); + use_intensity_map = true; + use_dual_thresholds = false; + in_threshold_01_a = stof(vArgs[i+1]); + in_threshold_01_b = stof(vArgs[i+2]); + out_thresh2_use_clipping = true; + if (vArgs[i] == "-cl") + out_thresh2_use_clipping_sigma = true; + else + out_thresh2_use_clipping_sigma = false; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number:\n" - " the maximum overlap allowed between a pair of detected blobs\n" - " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + " argument must be followed by 2 numbers.\n"); } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-spheres-nonmax-overlap-radial") + num_arguments_deleted = 3; + } - - // REMOVE THIS CRUFT - //else if ((vArgs[i] == "-spheres-nonmax-separation-sigma") || - // (vArgs[i] == "-blob-s-separation")) - //{ - // try { - // if (i+1 >= vArgs.size()) - // throw invalid_argument(""); - // nonmax_min_radial_separation_ratio = stof(vArgs[i+1]); - // //REMOVE THIS CRUFT - // //find_extrema_occlusion_ratio = stof(vArgs[i+1]); - // } - // catch (invalid_argument& exc) { - // throw InputErr("Error: The " + vArgs[i] + - // " argument must be followed by a number:\n" - // " the minimum allowed separation between a pair of detected blobs\n" - // " (as a fraction of the sum of their Gaussian widths (σ1+σ2))."); - // } - // num_arguments_deleted = 2; - //} //if (vArgs[i] == "-blob-s-separation") + else if ((vArgs[i] == "-thresh4") || + (vArgs[i] == "-thresh4-out")) { + try { + if (i+4 >= vArgs.size()) + throw invalid_argument(""); + use_intensity_map = true; + use_dual_thresholds = true; + in_threshold_01_a = stof(vArgs[i+1]); + in_threshold_01_b = stof(vArgs[i+2]); + in_threshold_10_a = stof(vArgs[i+3]); + in_threshold_10_b = stof(vArgs[i+4]); + bool all_increasing = ((in_threshold_01_a <= in_threshold_01_b) && + (in_threshold_01_b <= in_threshold_10_a) && + (in_threshold_10_a <= in_threshold_10_b)); + bool all_decreasing = ((in_threshold_01_a >= in_threshold_01_b) && + (in_threshold_01_b >= in_threshold_10_a) && + (in_threshold_10_a >= in_threshold_10_b)); + if (! (all_increasing || all_decreasing)) + throw invalid_argument(""); + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by 4 numbers\n" + " (These numbers must be either in increasing or decreasing order.)\n"); + } + num_arguments_deleted = 5; + } - else if ((vArgs[i] == "-radial-separation") || - (vArgs[i] == "-blob-separation") || - (vArgs[i] == "-blob-r-separation") || - (vArgs[i] == "-blobr-separation") || - (vArgs[i] == "-spheres-nonmax-separation-radius")) - { + else if ((vArgs[i] == "-thresh-interval") || + (vArgs[i] == "-thresh-interval-out")) { try { - if (i+1 >= vArgs.size()) + if (i+2 >= vArgs.size()) throw invalid_argument(""); - nonmax_min_radial_separation_ratio = stof(vArgs[i+1]); - //REMOVE THIS CRUFT - //nonmax_min_radial_separation_ratio = stof(vArgs[i+1]) * sqrt(3.0); - //find_extrema_occlusion_ratio = stof(vArgs[i+1]); + use_intensity_map = true; + use_dual_thresholds = true; + in_threshold_01_a = stof(vArgs[i+1]); + in_threshold_01_b = stof(vArgs[i+1]); + in_threshold_10_a = stof(vArgs[i+2]); + in_threshold_10_b = stof(vArgs[i+2]); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number:\n" - " the minimum allowed separation between a pair of detected blobs\n" - " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + " argument must be followed by 4 numbers.\n"); } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-radial-separation") - - + num_arguments_deleted = 3; + } - // REMOVE THIS CRUFT - //else if ((vArgs[i] == "-spheres-nonmax-separation-diameter") || - // (vArgs[i] == "-blob-d-separation") || - // (vArgs[i] == "-blob-separation")) - //{ - // try { - // if (i+1 >= vArgs.size()) - // throw invalid_argument(""); - // nonmax_min_radial_separation_ratio = stof(vArgs[i+1]) * 2.0; - // //REMOVE THIS CRUFT - // //nonmax_min_radial_separation_ratio = stof(vArgs[i+1]) * 2.0 * sqrt(3.0); - // //find_extrema_occlusion_ratio = stof(vArgs[i+1]); - // } - // catch (invalid_argument& exc) { - // throw InputErr("Error: The " + vArgs[i] + - // " argument must be followed by a number:\n" - // " the minimum allowed separation between a pair of detected blobs\n" - // " (as a fraction of the sum of their diameters, estimated using d=2r≈2σ√3)."); - // } - // num_arguments_deleted = 2; - //} //if (vArgs[i] == "-blobd-separation") + else if ((vArgs[i] == "-thresh-gauss") || + (vArgs[i] == "-thresh-gauss-out")) { + try { + if (i+2 >= vArgs.size()) + throw invalid_argument(""); + use_intensity_map = true; + use_gauss_thresholds = true; + out_thresh_gauss_x0 = stof(vArgs[i+1]); + out_thresh_gauss_sigma = stof(vArgs[i+2]); + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by 4 numbers.\n"); + } + num_arguments_deleted = 3; + } - else if ((vArgs[i] == "-blob") || - (vArgs[i] == "-blob-sigma") || - (vArgs[i] == "-blob-s") || - (vArgs[i] == "-blobs") || - (vArgs[i] == "-blob-radii") || - (vArgs[i] == "-blob-r") || - (vArgs[i] == "-blobr") || - (vArgs[i] == "-blob-diameters") || - (vArgs[i] == "-blob-d") || - (vArgs[i] == "-blob")) + else if ((vArgs[i] == "-gauss-aniso") || (vArgs[i] == "-ggauss-aniso")) { try { - if ((i+5 >= vArgs.size()) || + if ((i+3 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+3] == "") || (vArgs[i+3][0] == '-') || - (vArgs[i+4] == "") || (vArgs[i+4][0] == '-') || - (vArgs[i+5] == "") || (vArgs[i+5][0] == '-')) + (vArgs[i+3] == "") || (vArgs[i+3][0] == '-')) throw invalid_argument(""); + width_a[0] = stof(vArgs[i+1]); + width_a[1] = stof(vArgs[i+2]); + width_a[2] = stof(vArgs[i+3]); + width_b[0] = -1.0; + width_b[1] = -1.0; + width_b[2] = -1.0; + filter_type = GAUSS; + if (vArgs[i] == "-ggauss-aniso") + filter_type = GGAUSS; + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by 3 positive numbers:\n" + " s_x s_y s_z\n" + " the Gaussian widths in the X, Y, and Z direction.)\n"); + } + num_arguments_deleted = 4; + } //if (vArgs[i] == "-gauss-aniso") - string blob_file_name_base = vArgs[i+2]; - //if (EndsWith(blob_file_name_base, ".txt")) - // blob_file_name_base = - // blob_file_name_base.substr(0, - // blob_file_name_base.length()-4); - if ((vArgs[i+1] == "minima") || (vArgs[i+1] == "min")) { - blob_minima_file_name = blob_file_name_base; - blob_maxima_file_name = ""; // disable searching for maxima - score_upper_bound = 0.0; // disable searching for maxima - } - else if ((vArgs[i+1] == "maxima") || (vArgs[i+1] == "max")) { - blob_maxima_file_name = blob_file_name_base; - blob_minima_file_name = ""; // disable searching for minima - score_lower_bound = 0.0; // disable searching for maxima - } - else if (vArgs[i+1] == "all") { - blob_minima_file_name = blob_file_name_base + string(".minima.txt"); - blob_maxima_file_name = blob_file_name_base + string(".maxima.txt"); - if (score_lower_bound == 0.0) - // if searching for minima was disabled, then enable it - score_lower_bound = -std::numeric_limits::infinity(); - if (score_upper_bound == 0.0) - // if searching for maxima was disabled, then enable it - score_upper_bound = std::numeric_limits::infinity(); - } - else { - throw InputErr("Error: The 1st parameter to the \"" + vArgs[i] + "\" argument must be one of:\n" - " \"minima\", \"maxima\", (or) \"all\"\n" - "\n" - " (It indicates whether you are try to detect\n" - " dark blobs, bright blobs, or both.)\n"); - } - - float blob_width_min = stof(vArgs[i+3]); - float blob_width_max = stof(vArgs[i+4]); - if ((blob_width_min <= 0.0) || - (blob_width_max <= 0.0) || - (blob_width_min >= blob_width_max)) - throw invalid_argument(""); - // The old version of the code allows user to set "N" directly - //int N = stof(vArgs[i+5]); - //if (N < 3) - // throw InputErr("The 4th parameter to the \"" + vArgs[i] + "\" argument must be >= 3."); - //float growth_ratio = pow(blob_width_max / blob_width_min, 1.0/(N-1)); - - // In the new version of the code, the user specifies "growth_ratio" - float growth_ratio = stof(vArgs[i+5]); - if (growth_ratio <= 1.0) + else if ((vArgs[i] == "-gauss") || (vArgs[i] == "-ggauss")) + { + try { + if ((i+1 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) throw invalid_argument(""); - int N = 1+ceil(log(blob_width_max/blob_width_min) / log(growth_ratio)); - growth_ratio = pow(blob_width_max/blob_width_min, 1.0/N); - blob_diameters.resize(N); - - // REMOVE THIS CRUFT: - // Optional: - // Make sure the difference between the width of the Gaussians used - // to approximate the Laplacian-of-Gaussian (delta_sigma_over_sigma), - // is no larger than the difference between the successive widths - // in the list of Gaussian blurs we will apply to our source image. - // (Actually, I make sure it is no larger than 1/3 this difference.) - //if (((growth_ratio-1.0)/3 < delta_sigma_over_sigma) - // && - // (! user_set_delta_manually)) - // delta_sigma_over_sigma = (growth_ratio-1.0)/3; - - blob_width_multiplier = 1.0; - if ((vArgs[i] == "-blob-sigma") || (vArgs[i] == "-blob-s")) - // (for a solid uniform 3-D sphere) - blob_width_multiplier = 2.0*sqrt(3.0); - //REMOVE THIS CRUFT - //blob_width_multiplier = 1.0/sqrt(3.0); - if ((vArgs[i] == "-blob-radii") || (vArgs[i] == "-blob-r") || (vArgs[i] == "-blobr")) - // (for a solid uniform 3-D sphere) - blob_width_multiplier = 2.0; - //REMOVE THIS CRUFT - //blob_width_multiplier = 1.0/sqrt(3.0); - else if ((vArgs[i] == "-blob-diameters") || (vArgs[i] == "-blob-d")) - // (for a solid uniform 3-D sphere) - blob_width_multiplier = 1.0; - //REMOVE THIS CRUFT - //blob_width_multiplier = 1.0/(2.0*sqrt(3.0)); - blob_diameters[0] = blob_width_min * blob_width_multiplier; - for (int n = 1; n < N; n++) { - blob_diameters[n] = blob_diameters[n-1] * growth_ratio; - } - m_exp = 2.0; // (<--not necessary, we ignore these parameters anyway) - n_exp = 2.0; // (<--not necessary, we ignore these parameters anyway) - filter_type = BLOB; + width_a[0] = stof(vArgs[i+1]); + width_a[1] = width_a[0]; + width_a[2] = width_a[0]; + width_b[0] = -1.0; + width_b[1] = -1.0; + width_b[2] = -1.0; + filter_type = GAUSS; + if (vArgs[i] == "-ggauss") + filter_type = GGAUSS; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a type (either\n" - " \"minima\", \"maxima\", or \"all\")\n" - " as well as a file name, followed by 3 positive numbers.\n" - "For example:\n" - " minima file_name.txt 200.0 280.0 1.02\n"); + " argument must be followed by a positive number (\"s\"),\n" + " the Gaussian width\n"); } - num_arguments_deleted = 6; - } //if (vArgs[i] == "-blob") - + num_arguments_deleted = 2; + } //if (vArgs[i] == "-gauss") - else if ((vArgs[i] == "-discard-blobs") || - (vArgs[i] == "-blob-nonmax") || - (vArgs[i] == "-blobs-nonmax")) { + else if ((vArgs[i] == "-dog-aniso") || + (vArgs[i] == "-dogg-aniso")) + { try { - if ((i+2 >= vArgs.size()) || + if ((i+6 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || - (vArgs[i+1] == vArgs[i+2])) + (vArgs[i+3] == "") || (vArgs[i+3][0] == '-') || + (vArgs[i+4] == "") || (vArgs[i+4][0] == '-') || + (vArgs[i+5] == "") || (vArgs[i+5][0] == '-') || + (vArgs[i+6] == "") || (vArgs[i+6][0] == '-')) throw invalid_argument(""); - in_coords_file_name = vArgs[i+1]; - out_coords_file_name = vArgs[i+2]; + width_a[0] = stof(vArgs[i+1]); + width_a[1] = stof(vArgs[i+2]); + width_a[2] = stof(vArgs[i+3]); + width_b[0] = stof(vArgs[i+4]); + width_b[1] = stof(vArgs[i+5]); + width_b[2] = stof(vArgs[i+6]); + //if (width_b[0] <= width_a[0]) + // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" + // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); + filter_type = DOG; + if (vArgs[i] == "-dogg-aniso") + filter_type = DOGG; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by two different file names\n"); + " argument must be followed by 6 positive numbers.\n"); } - filter_type = SPHERE_NONMAX_SUPPRESSION; - num_arguments_deleted = 3; - } //if (vArgs[i] == "-discard-blobs") - - + num_arguments_deleted = 7; + } //if (vArgs[i] == "-dog-aniso") - else if (vArgs[i] == "-auto-thresh") - { - auto_thresh_score = true; - try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) - throw invalid_argument(""); - string list_of_attributes = vArgs[i+1]; - if (list_of_attributes != "score") - throw invalid_argument(""); - } - catch (invalid_argument& exc) { - throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by \"score\"\n" - " (In the future, it may be followed by a list of attributes surrounded\n" - " in quotes. But currently only the \"score\" attribute is supported.\n" - " Thresholds for these attributes will be determined automatically.)\n"); - } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-auto-thresh") - else if (vArgs[i] == "-supervised") + else if ((vArgs[i] == "-dog") || + (vArgs[i] == "-dogg")) { try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || - (i+2 >= vArgs.size()) || (vArgs[i+2] == "")) + if ((i+2 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || + (vArgs[i+2] == "") || (vArgs[i+2][0] == '-')) throw invalid_argument(""); - training_pos_fname = vArgs[i+1]; - training_neg_fname = vArgs[i+2]; + width_a[0] = stof(vArgs[i+1]); + width_a[1] = width_a[0]; + width_a[2] = width_a[0]; + width_b[0] = stof(vArgs[i+2]); + width_b[1] = width_b[0]; + width_b[2] = width_b[0]; + //if (width_b[0] <= width_a[0]) + // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" + // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); + filter_type = DOG; + if (vArgs[i] == "-dogg") + filter_type = DOGG; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by two file names:\n" - " FILE_POS.txt FILE_NEG.txt\n" - " containing positive and negative training data, respectively.\n"); + " argument must be followed by 2 positive numbers.\n"); } num_arguments_deleted = 3; - } //if (vArgs[i] == "-supervised") - - - else if (vArgs[i] == "-supervised-multi") - { - filter_type = SPHERE_NONMAX_SUPERVISED_MULTI; - try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) - throw invalid_argument(""); - string fname = vArgs[i+1]; - - { //Read multi_training_meta_fname - fstream f; - f.open(fname, ios::in); - const char comment_char = '#'; - string strLine; - size_t i_line = 1; - while (getline(f, strLine)) - { - // ignore text after comments - size_t ic = strLine.find(comment_char); - if (ic != string::npos) - strLine = strLine.substr(0, ic); - stringstream ssLine(strLine); - vector tokens; - string x; - try { - while (ssLine >> x) { - tokens.push_back(x); - } - } // try { - catch ( ... ) { - stringstream err_msg; - err_msg << "Error: Read error (invalid entry?) in file:\n" - << " \"" << fname << "\", line: " << i_line << "\n"; - throw VisfdErr(err_msg.str()); - } - if (tokens.size() == 0) - continue; - else if (tokens.size() != 3) { - stringstream err_msg; - err_msg << "Error: Format error in file:\n" - << " \"" << fname << "\", line: " << i_line << "\n" - << "\n" - << " Expected 3 file names in every line in this file:\n" - << "\n" - << " training_pos_file training_neg_file blob_info_list.txt\n"; - throw VisfdErr(err_msg.str()); - } - multi_training_pos_fnames.push_back(tokens[0]); - multi_training_neg_fnames.push_back(tokens[1]); - multi_in_coords_file_names.push_back(tokens[2]); - } //while (getline(f, strLine)) - f.close(); - } //Read multi_training_meta_fname - } - catch (invalid_argument& exc) { - throw InputErr("Error: The " + vArgs[i] + - " argument must be followed a file name. Example: FILE_TRAINING_SETS.txt\n" - "\n" - "File format details:\n" - " This should be a 3-column text file containing file names. Example:\n" - " training_pos_1.txt training_crds_neg_1.txt blob_info_1.txt\n" - " training_pos_2.txt training_crds_neg_2.txt blob_info_2.txt\n" - " training_pos_3.txt training_crds_neg_3.txt blob_info_3.txt\n" - " : : :\n" - " training_pos_N.txt training_crds_neg_N.txt blob_info_N.txt\n" - "\n" - " containing positive and negative training data, and a list of detected blobs\n" - " for each of the N images you want to simultaneously analyze.\n"); - } - - num_arguments_deleted = 2; + } //if (vArgs[i] == "-dog") - } //if (vArgs[i] == "-supervised-multi") - else if ((vArgs[i] == "-minima-threshold") || - (vArgs[i] == "-score-upper-bound")) + #ifndef DISABLE_DOGGXY + else if (vArgs[i] == "-doggxy-aniso") { try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) + if ((i+5 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || + (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || + (vArgs[i+3] == "") || (vArgs[i+3][0] == '-') || + (vArgs[i+4] == "") || (vArgs[i+4][0] == '-') || + (vArgs[i+5] == "") || (vArgs[i+5][0] == '-')) throw invalid_argument(""); - score_upper_bound = stof(vArgs[i+1]); - //score_lower_bound = -std::numeric_limits::infinity(); - score_bounds_are_ratios = false; + width_a[0] = stof(vArgs[i+1]); + width_a[1] = stof(vArgs[i+2]); + width_b[0] = stof(vArgs[i+3]); + width_b[1] = stof(vArgs[i+4]); + //The "-doggxy" filter is a Difference-of-Generalized-Gaussians + //in the X,Y directions, multiplied by an ordinary Gaussian in + //the Z direction. + width_a[2] = stof(vArgs[i+5]); + width_b[2] = -1.0; //(This disables the second (negative) Gaussian) + //if (width_b[0] <= width_a[0]) + // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" + // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); + filter_type = DOGGXY; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number number.\n"); + " argument must be followed by 5 positive numbers:\n" + " a_x a_y b_x b_y a_z\n" + " (I.E., the \"A\" and \"B\" Gaussian widths in X and Y directions,\n" + " followed by the Gaussian width in the Z direction.)\n"); } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-minima-threshold") - + num_arguments_deleted = 6; + } //if (vArgs[i] == "-doggxy-aniso") - else if ((vArgs[i] == "-maxima-threshold") || - (vArgs[i] == "-score-lower-bound")) - { - try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) - throw invalid_argument(""); - score_lower_bound = stof(vArgs[i+1]); - //score_upper_bound = std::numeric_limits::infinity(); - score_bounds_are_ratios = false; - } - catch (invalid_argument& exc) { - throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number number.\n"); - } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-maxima-threshold") - - else if ((vArgs[i] == "-minima-ratio") || - (vArgs[i] == "-score-lower-bound-ratio")) + else if (vArgs[i] == "-doggxy") { try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) + if ((i+5 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || + (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || + (vArgs[i+3] == "") || (vArgs[i+3][0] == '-')) throw invalid_argument(""); - score_upper_bound = stof(vArgs[i+1]); - //score_lower_bound = -std::numeric_limits::infinity(); - score_bounds_are_ratios = true; + width_a[0] = stof(vArgs[i+1]); + width_a[1] = width_a[0]; + width_b[0] = stof(vArgs[i+2]); + width_b[1] = width_b[0]; + //The "-doggxy" filter is a Difference-of-Generalized-Gaussians + //in the X,Y directions, multiplied by an ordinary Gaussian in + //the Z direction. + width_a[2] = stof(vArgs[i+3]); + width_b[2] = -1.0; //(This disables the second (negative) Gaussian) + //if (width_b[0] <= width_a[0]) + // throw InputErr("Error: The two arguments following " + vArgs[i] + " must be\n" + // " increasing. (Ie., the 2nd argument must be > 1st argument.)\n"); + filter_type = DOGGXY; } catch (invalid_argument& exc) { - throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number number.\n"); + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by 3 positive numbers:\n" + " a_xy b_xy a_z\n" + " (I.E., the \"A\" and \"-B\" Gaussian widths in the XY plane,\n" + " followed by the Gaussian width in the Z direction.)\n"); } - num_arguments_deleted = 2; - } //if (vArgs[i] == "-minima_ratio") + num_arguments_deleted = 4; + } //if (vArgs[i] == "-doggxy") + #endif //#ifndef DISABLE_DOGGXY - else if ((vArgs[i] == "-maxima-ratio") || - (vArgs[i] == "-score-upper-bound-ratio")) + else if ((vArgs[i] == "-log") || + (vArgs[i] == "-log-r") || + (vArgs[i] == "-log-d")) { try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) + if ((i+1 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) throw invalid_argument(""); - score_lower_bound = stof(vArgs[i+1]); - //score_upper_bound = std::numeric_limits::infinity(); - score_bounds_are_ratios = true; + float blob_sigma_multiplier = 1.0; + blob_sigma_multiplier = 1.0; + if (vArgs[i] == "-log") + blob_sigma_multiplier = 1.0; + if (vArgs[i] == "-log-r") + // (for a solid uniform 3-D sphere) + blob_sigma_multiplier = 1.0/sqrt(3.0); + if (vArgs[i] == "-log-d") + // (for a solid uniform 3-D sphere) + blob_sigma_multiplier = 1.0/(2.0*sqrt(3.0)); + log_width[0] = stof(vArgs[i+1]) * blob_sigma_multiplier; + log_width[1] = log_width[0]; + log_width[2] = log_width[0]; + m_exp = 2.0; + n_exp = 2.0; + filter_type = LOG_DOG; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number number.\n"); + " argument must be followed by a positive number:\n" + " the (approximate) width of the objects of interest.\n"); } num_arguments_deleted = 2; - } //if (vArgs[i] == "-maxima-ratio") - - - - else if ((vArgs[i] == "-spheres-nonmax-radii-range") || - (vArgs[i] == "-sphere-nonmax-radii-range")) { - try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "-") || (vArgs[i+1][0] == '-') || - (i+2 >= vArgs.size()) || (vArgs[i+2] == "-") || (vArgs[i+2][0] == '-')) - throw invalid_argument(""); - sphere_diameters_lower_bound = stof(vArgs[i+1]); - sphere_diameters_upper_bound = stof(vArgs[i+2]); - } - catch (invalid_argument& exc) { - throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number number.\n"); - } - num_arguments_deleted = 3; - } //if (vArgs[i] == "-spheres-nonmax-radii-range") + } //if (vArgs[i] == "-log") - else if ((vArgs[i] == "-spheres-nonmax-score-range") || - (vArgs[i] == "-sphere-nonmax-score-range")) { + else if (vArgs[i] == "-log-aniso") + { try { - if (i+2 >= vArgs.size()) + if ((i+3 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || + (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || + (vArgs[i+3] == "") || (vArgs[i+3][0] == '-')) throw invalid_argument(""); - score_lower_bound = stof(vArgs[i+1]); - score_upper_bound = stof(vArgs[i+2]); + log_width[0] = stof(vArgs[i+1]); + log_width[1] = stof(vArgs[i+2]); + log_width[2] = stof(vArgs[i+3]); + m_exp = 2.0; + n_exp = 2.0; + filter_type = LOG_DOG; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by two numbers.\n"); + " argument must be followed by 3 positive numbers:\n" + " the approximate widths of the objects of interest in the XYZ directions\n" + " (Due to the missing-wedge artifact they might differ.\n"); } - num_arguments_deleted = 3; - } //if (vArgs[i] == "-spheres-nonmax-score-range") + num_arguments_deleted = 4; + } //if (vArgs[i] == "-log-aniso") + + + @@ -1373,6 +1228,7 @@ Settings::ParseArgs(vector& vArgs) } //if (vArgs[i] == "-gdog-exponents") + else if ((vArgs[i] == "-gauss-exponent") || (vArgs[i] == "-exponent")) { @@ -1395,288 +1251,470 @@ Settings::ParseArgs(vector& vArgs) } //if (vArgs[i] == "-gauss-exponent") - //Commenting out: The user cannot set filter_truncate_halfwidth directly anymore: - //else if (vArgs[i] == "-window") { - // if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) - // throw InputErr("Error: The " + vArgs[i] + - // " argument must be followed by either 1 or 3 positive integers.\n"); - // filter_truncate_halfwidth[0] = stoi(vArgs[i+1]); - // num_arguments_deleted = 2; - // if ((i+4 >= vArgs.size()) && - // (vArgs[i+2][0] != '-') && (vArgs[i+2] != "") && - // (vArgs[i+3][0] != '-') && (vArgs[i+3] != "")) { - // filter_truncate_halfwidth[1] = stoi(vArgs[i+2]); - // filter_truncate_halfwidth[2] = stoi(vArgs[i+3]); - // num_arguments_deleted = 4; - // } - // else { - // filter_truncate_halfwidth[1] = filter_truncate_halfwidth[0]; - // filter_truncate_halfwidth[2] = filter_truncate_halfwidth[0]; - // } - //} //if (vArgs[i] == "-window") + + else if ((vArgs[i] == "-spheres-nonmax-overlap") || + (vArgs[i] == "-max-volume-overlap") || + (vArgs[i] == "-max-overlap")) + { + try { + if (i+1 >= vArgs.size()) + throw invalid_argument(""); + nonmax_max_volume_overlap_large = stof(vArgs[i+1]); + nonmax_min_radial_separation_ratio = 0.0; + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by a number:\n" + " the maximum overlap allowed between a pair of detected blobs\n" + " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-spheres-nonmax-overlap") - else if (vArgs[i] == "-truncate") + + else if ((vArgs[i] == "-spheres-nonmax-overlap-small") || + (vArgs[i] == "-max-volume-overlap-small") || + (vArgs[i] == "-max-overlap-small")) { try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) + if (i+1 >= vArgs.size()) throw invalid_argument(""); - filter_truncate_ratio = stof(vArgs[i+1]); - filter_truncate_threshold=-1.0; //(disables)override any filter_truncate_threshold settings - //filter_truncate_ratio_exp = exp(-pow(filter_truncate_ratio, - // n_exp)); + nonmax_max_volume_overlap_small = stof(vArgs[i+1]); + nonmax_min_radial_separation_ratio = 0.0; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number.\n"); + " argument must be followed by a number:\n" + " the maximum overlap allowed between a pair of detected blobs\n" + " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); } num_arguments_deleted = 2; - } // if (vArgs[i] == "-truncate") + } //if (vArgs[i] == "-spheres-nonmax-overlap-small") - else if (vArgs[i] == "-truncate-threshold") + + else if ((vArgs[i] == "-spheres-nonmax-overlap-radial") || + (vArgs[i] == "-max-overlap-radial")) { try { - if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || (vArgs[i+1][0] == '-')) + if (i+1 >= vArgs.size()) throw invalid_argument(""); - filter_truncate_threshold = stof(vArgs[i+1]); - filter_truncate_ratio = -1.0; //(disables) override any filter_truncate_ratio settings + float max_overlap = stof(vArgs[i+1]); + nonmax_min_radial_separation_ratio = (1.0 - max_overlap); + //REMOVE THIS CRUFT: + //nonmax_min_radial_separation_ratio = (1.0 - max_overlap) * sqrt(3.0); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number.\n"); - } + " argument must be followed by a number:\n" + " the maximum overlap allowed between a pair of detected blobs\n" + " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + } num_arguments_deleted = 2; - } // if (vArgs[i] == "-truncate-thresold") + } //if (vArgs[i] == "-spheres-nonmax-overlap-radial") + + else if ((vArgs[i] == "-radial-separation") || + (vArgs[i] == "-blob-separation") || + (vArgs[i] == "-blob-r-separation") || + (vArgs[i] == "-blobr-separation") || + (vArgs[i] == "-spheres-nonmax-separation-radius")) + { + try { + if (i+1 >= vArgs.size()) + throw invalid_argument(""); + nonmax_min_radial_separation_ratio = stof(vArgs[i+1]); + //REMOVE THIS CRUFT + //nonmax_min_radial_separation_ratio = stof(vArgs[i+1]) * sqrt(3.0); + //find_extrema_occlusion_ratio = stof(vArgs[i+1]); + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by a number:\n" + " the minimum allowed separation between a pair of detected blobs\n" + " (as a fraction of the sum of their radii, estimated using r≈σ√3)."); + } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-radial-separation") + - else if (vArgs[i] == "-rescale") + else if ((vArgs[i] == "-blob") || + (vArgs[i] == "-blob-sigma") || + (vArgs[i] == "-blob-s") || + (vArgs[i] == "-blobs") || + (vArgs[i] == "-blob-radii") || + (vArgs[i] == "-blob-r") || + (vArgs[i] == "-blobr") || + (vArgs[i] == "-blob-diameters") || + (vArgs[i] == "-blob-d") || + (vArgs[i] == "-blob")) { + try { + if ((i+5 >= vArgs.size()) || + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || + (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || + (vArgs[i+3] == "") || (vArgs[i+3][0] == '-') || + (vArgs[i+4] == "") || (vArgs[i+4][0] == '-') || + (vArgs[i+5] == "") || (vArgs[i+5][0] == '-')) + throw invalid_argument(""); + + string blob_file_name_base = vArgs[i+2]; + //if (EndsWith(blob_file_name_base, ".txt")) + // blob_file_name_base = + // blob_file_name_base.substr(0, + // blob_file_name_base.length()-4); + + if ((vArgs[i+1] == "minima") || (vArgs[i+1] == "min")) { + blob_minima_file_name = blob_file_name_base; + blob_maxima_file_name = ""; // disable searching for maxima + score_upper_bound = 0.0; // disable searching for maxima + } + else if ((vArgs[i+1] == "maxima") || (vArgs[i+1] == "max")) { + blob_maxima_file_name = blob_file_name_base; + blob_minima_file_name = ""; // disable searching for minima + score_lower_bound = 0.0; // disable searching for maxima + } + else if (vArgs[i+1] == "all") { + blob_minima_file_name = blob_file_name_base + string(".minima.txt"); + blob_maxima_file_name = blob_file_name_base + string(".maxima.txt"); + if (score_lower_bound == 0.0) + // if searching for minima was disabled, then enable it + score_lower_bound = -std::numeric_limits::infinity(); + if (score_upper_bound == 0.0) + // if searching for maxima was disabled, then enable it + score_upper_bound = std::numeric_limits::infinity(); + } + else { + throw InputErr("Error: The 1st parameter to the \"" + vArgs[i] + "\" argument must be one of:\n" + " \"minima\", \"maxima\", (or) \"all\"\n" + "\n" + " (It indicates whether you are try to detect\n" + " dark blobs, bright blobs, or both.)\n"); + } + + float blob_width_min = stof(vArgs[i+3]); + float blob_width_max = stof(vArgs[i+4]); + if ((blob_width_min <= 0.0) || + (blob_width_max <= 0.0) || + (blob_width_min >= blob_width_max)) + throw invalid_argument(""); + + // The old version of the code allows user to set "N" directly + //int N = stof(vArgs[i+5]); + //if (N < 3) + // throw InputErr("The 4th parameter to the \"" + vArgs[i] + "\" argument must be >= 3."); + //float growth_ratio = pow(blob_width_max / blob_width_min, 1.0/(N-1)); + + // In the new version of the code, the user specifies "growth_ratio" + float growth_ratio = stof(vArgs[i+5]); + if (growth_ratio <= 1.0) + throw invalid_argument(""); + int N = 1+ceil(log(blob_width_max/blob_width_min) / log(growth_ratio)); + growth_ratio = pow(blob_width_max/blob_width_min, 1.0/N); + blob_diameters.resize(N); + + // REMOVE THIS CRUFT: + // Optional: + // Make sure the difference between the width of the Gaussians used + // to approximate the Laplacian-of-Gaussian (delta_sigma_over_sigma), + // is no larger than the difference between the successive widths + // in the list of Gaussian blurs we will apply to our source image. + // (Actually, I make sure it is no larger than 1/3 this difference.) + //if (((growth_ratio-1.0)/3 < delta_sigma_over_sigma) + // && + // (! user_set_delta_manually)) + // delta_sigma_over_sigma = (growth_ratio-1.0)/3; + + blob_width_multiplier = 1.0; + if ((vArgs[i] == "-blob-sigma") || (vArgs[i] == "-blob-s")) + // (for a solid uniform 3-D sphere) + blob_width_multiplier = 2.0*sqrt(3.0); + //REMOVE THIS CRUFT + //blob_width_multiplier = 1.0/sqrt(3.0); + if ((vArgs[i] == "-blob-radii") || (vArgs[i] == "-blob-r") || (vArgs[i] == "-blobr")) + // (for a solid uniform 3-D sphere) + blob_width_multiplier = 2.0; + //REMOVE THIS CRUFT + //blob_width_multiplier = 1.0/sqrt(3.0); + else if ((vArgs[i] == "-blob-diameters") || (vArgs[i] == "-blob-d")) + // (for a solid uniform 3-D sphere) + blob_width_multiplier = 1.0; + //REMOVE THIS CRUFT + //blob_width_multiplier = 1.0/(2.0*sqrt(3.0)); + blob_diameters[0] = blob_width_min * blob_width_multiplier; + for (int n = 1; n < N; n++) { + blob_diameters[n] = blob_diameters[n-1] * growth_ratio; + } + m_exp = 2.0; // (<--not necessary, we ignore these parameters anyway) + n_exp = 2.0; // (<--not necessary, we ignore these parameters anyway) + filter_type = BLOB; + } + catch (invalid_argument& exc) { + throw InputErr("Error: The " + vArgs[i] + + " argument must be followed by a type (either\n" + " \"minima\", \"maxima\", or \"all\")\n" + " as well as a file name, followed by 3 positive numbers.\n" + "For example:\n" + " minima file_name.txt 200.0 280.0 1.02\n"); + } + num_arguments_deleted = 6; + } //if (vArgs[i] == "-blob") + + + + + else if ((vArgs[i] == "-discard-blobs") || + (vArgs[i] == "-blob-nonmax") || + (vArgs[i] == "-blobs-nonmax")) { try { if ((i+2 >= vArgs.size()) || - (vArgs[i+1] == "") || - (vArgs[i+2] == "")) + (vArgs[i+1] == "") || (vArgs[i+1][0] == '-') || + (vArgs[i+2] == "") || (vArgs[i+2][0] == '-') || + (vArgs[i+1] == vArgs[i+2])) throw invalid_argument(""); - use_intensity_map = true; - use_rescale_multiply = true; - out_rescale_multiply = stof(vArgs[i+1]); - out_rescale_offset = stof(vArgs[i+2]); + in_coords_file_name = vArgs[i+1]; + out_coords_file_name = vArgs[i+2]; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 2 numbers:\n" - " outA outB\n" - " (the desired minimum and maximum voxel intensity values for the final image)\n"); + " argument must be followed by two different file names\n"); } + filter_type = SPHERE_NONMAX_SUPPRESSION; num_arguments_deleted = 3; - } // if (vArgs[i] == "-rescale") + } //if (vArgs[i] == "-discard-blobs") - else if (vArgs[i] == "-fill") + + + else if (vArgs[i] == "-auto-thresh") { + auto_thresh_score = true; try { if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); - use_intensity_map = true; - use_rescale_multiply = true; - out_rescale_multiply = 0.0; - out_rescale_offset = stof(vArgs[i+1]); + string list_of_attributes = vArgs[i+1]; + if (list_of_attributes != "score") + throw invalid_argument(""); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by a number."); + " argument must be followed by \"score\"\n" + " (In the future, it may be followed by a list of attributes surrounded\n" + " in quotes. But currently only the \"score\" attribute is supported.\n" + " Thresholds for these attributes will be determined automatically.)\n"); } num_arguments_deleted = 2; - } // if (vArgs[i] == "-rescale") + } //if (vArgs[i] == "-auto-thresh") - else if ((vArgs[i] == "-thresh-range") || - (vArgs[i] == "-thresh-range-out")) + + else if (vArgs[i] == "-supervised") { try { - if ((i+2 >= vArgs.size()) || - (vArgs[i+1] == "") || - (vArgs[i+2] == "")) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "") || + (i+2 >= vArgs.size()) || (vArgs[i+2] == "")) throw invalid_argument(""); - out_thresh_a_value = stof(vArgs[i+1]); - out_thresh_b_value = stof(vArgs[i+2]); + training_pos_fname = vArgs[i+1]; + training_neg_fname = vArgs[i+2]; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 2 numbers:\n" - " outA outB\n" - " (the desired minimum and maximum voxel intensity values for the final image)\n"); + " argument must be followed by two file names:\n" + " FILE_POS.txt FILE_NEG.txt\n" + " containing positive and negative training data, respectively.\n"); } num_arguments_deleted = 3; - } // if (vArgs[i] == "-outab") + } //if (vArgs[i] == "-supervised") - else if (vArgs[i] == "-rescale-min-max") + else if (vArgs[i] == "-supervised-multi") { + filter_type = SPHERE_NONMAX_SUPERVISED_MULTI; try { - if ((i+2 >= vArgs.size()) || - (vArgs[i+1] == "") || - (vArgs[i+2] == "")) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); - rescale_min_max_out = true; - out_rescale_max = stof(vArgs[i+1]); - out_rescale_min = stof(vArgs[i+2]); + string fname = vArgs[i+1]; + + { //Read multi_training_meta_fname + fstream f; + f.open(fname, ios::in); + const char comment_char = '#'; + string strLine; + size_t i_line = 1; + while (getline(f, strLine)) + { + // ignore text after comments + size_t ic = strLine.find(comment_char); + if (ic != string::npos) + strLine = strLine.substr(0, ic); + stringstream ssLine(strLine); + vector tokens; + string x; + try { + while (ssLine >> x) { + tokens.push_back(x); + } + } // try { + catch ( ... ) { + stringstream err_msg; + err_msg << "Error: Read error (invalid entry?) in file:\n" + << " \"" << fname << "\", line: " << i_line << "\n"; + throw VisfdErr(err_msg.str()); + } + if (tokens.size() == 0) + continue; + else if (tokens.size() != 3) { + stringstream err_msg; + err_msg << "Error: Format error in file:\n" + << " \"" << fname << "\", line: " << i_line << "\n" + << "\n" + << " Expected 3 file names in every line in this file:\n" + << "\n" + << " training_pos_file training_neg_file blob_info_list.txt\n"; + throw VisfdErr(err_msg.str()); + } + multi_training_pos_fnames.push_back(tokens[0]); + multi_training_neg_fnames.push_back(tokens[1]); + multi_in_coords_file_names.push_back(tokens[2]); + } //while (getline(f, strLine)) + f.close(); + } //Read multi_training_meta_fname } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 2 numbers:\n" - " outA outB\n" - " (the desired minimum and maximum voxel intensity values for the final image)\n"); + " argument must be followed a file name. Example: FILE_TRAINING_SETS.txt\n" + "\n" + "File format details:\n" + " This should be a 3-column text file containing file names. Example:\n" + " training_pos_1.txt training_crds_neg_1.txt blob_info_1.txt\n" + " training_pos_2.txt training_crds_neg_2.txt blob_info_2.txt\n" + " training_pos_3.txt training_crds_neg_3.txt blob_info_3.txt\n" + " : : :\n" + " training_pos_N.txt training_crds_neg_N.txt blob_info_N.txt\n" + "\n" + " containing positive and negative training data, and a list of detected blobs\n" + " for each of the N images you want to simultaneously analyze.\n"); } - num_arguments_deleted = 3; - } // if (vArgs[i] == "-rescale-min-max") + num_arguments_deleted = 2; - else if ((vArgs[i] == "-no-rescale") || (vArgs[i] == "-norescale")) - { - rescale_min_max_out = false; - in_threshold_01_a = 1.0; - in_threshold_01_b = 1.0; - num_arguments_deleted = 1; - } // if (vArgs[i] == "-norescale") - + } //if (vArgs[i] == "-supervised-multi") - else if ((vArgs[i] == "-invert") || - (vArgs[i] == "-inv")) - { - invert_output = true; - num_arguments_deleted = 1; - } // if (vArgs[i] == "-invert") - else if ((vArgs[i] == "-thresh") || - (vArgs[i] == "-thresh-out")) { + else if ((vArgs[i] == "-minima-threshold") || + (vArgs[i] == "-score-upper-bound")) + { try { - if (i+1 >= vArgs.size()) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); - use_intensity_map = true; - use_dual_thresholds = false; - in_threshold_01_a = stof(vArgs[i+1]); - in_threshold_01_b = in_threshold_01_a; + score_upper_bound = stof(vArgs[i+1]); + //score_lower_bound = -std::numeric_limits::infinity(); + score_bounds_are_ratios = false; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 1 number.\n"); + " argument must be followed by a number number.\n"); } num_arguments_deleted = 2; - } + } //if (vArgs[i] == "-minima-threshold") - else if ((vArgs[i] == "-thresh2") || - (vArgs[i] == "-thresh2-out")) { + + else if ((vArgs[i] == "-maxima-threshold") || + (vArgs[i] == "-score-lower-bound")) + { try { - if (i+2 >= vArgs.size()) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); - use_intensity_map = true; - use_dual_thresholds = false; - in_threshold_01_a = stof(vArgs[i+1]); - in_threshold_01_b = stof(vArgs[i+2]); - out_thresh2_use_clipping = false; + score_lower_bound = stof(vArgs[i+1]); + //score_upper_bound = std::numeric_limits::infinity(); + score_bounds_are_ratios = false; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 2 numbers.\n"); + " argument must be followed by a number number.\n"); } - num_arguments_deleted = 3; - } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-maxima-threshold") - else if ((vArgs[i] == "-clip") || - (vArgs[i] == "-cl")) { + + else if ((vArgs[i] == "-minima-ratio") || + (vArgs[i] == "-score-lower-bound-ratio")) + { try { - if (i+2 >= vArgs.size()) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); - use_intensity_map = true; - use_dual_thresholds = false; - in_threshold_01_a = stof(vArgs[i+1]); - in_threshold_01_b = stof(vArgs[i+2]); - out_thresh2_use_clipping = true; - if (vArgs[i] == "-cl") - out_thresh2_use_clipping_sigma = true; - else - out_thresh2_use_clipping_sigma = false; + score_upper_bound = stof(vArgs[i+1]); + //score_lower_bound = -std::numeric_limits::infinity(); + score_bounds_are_ratios = true; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 2 numbers.\n"); + " argument must be followed by a number number.\n"); } - num_arguments_deleted = 3; - } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-minima_ratio") - else if ((vArgs[i] == "-thresh4") || - (vArgs[i] == "-thresh4-out")) { + + else if ((vArgs[i] == "-maxima-ratio") || + (vArgs[i] == "-score-upper-bound-ratio")) + { try { - if (i+4 >= vArgs.size()) - throw invalid_argument(""); - use_intensity_map = true; - use_dual_thresholds = true; - in_threshold_01_a = stof(vArgs[i+1]); - in_threshold_01_b = stof(vArgs[i+2]); - in_threshold_10_a = stof(vArgs[i+3]); - in_threshold_10_b = stof(vArgs[i+4]); - bool all_increasing = ((in_threshold_01_a <= in_threshold_01_b) && - (in_threshold_01_b <= in_threshold_10_a) && - (in_threshold_10_a <= in_threshold_10_b)); - bool all_decreasing = ((in_threshold_01_a >= in_threshold_01_b) && - (in_threshold_01_b >= in_threshold_10_a) && - (in_threshold_10_a >= in_threshold_10_b)); - if (! (all_increasing || all_decreasing)) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "")) throw invalid_argument(""); + score_lower_bound = stof(vArgs[i+1]); + //score_upper_bound = std::numeric_limits::infinity(); + score_bounds_are_ratios = true; } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 4 numbers\n" - " (These numbers must be either in increasing or decreasing order.)\n"); + " argument must be followed by a number number.\n"); } - num_arguments_deleted = 5; - } + num_arguments_deleted = 2; + } //if (vArgs[i] == "-maxima-ratio") - else if ((vArgs[i] == "-thresh-interval") || - (vArgs[i] == "-thresh-interval-out")) { + + else if ((vArgs[i] == "-spheres-nonmax-radii-range") || + (vArgs[i] == "-sphere-nonmax-radii-range")) { try { - if (i+2 >= vArgs.size()) + if ((i+1 >= vArgs.size()) || (vArgs[i+1] == "-") || (vArgs[i+1][0] == '-') || + (i+2 >= vArgs.size()) || (vArgs[i+2] == "-") || (vArgs[i+2][0] == '-')) throw invalid_argument(""); - use_intensity_map = true; - use_dual_thresholds = true; - in_threshold_01_a = stof(vArgs[i+1]); - in_threshold_01_b = stof(vArgs[i+1]); - in_threshold_10_a = stof(vArgs[i+2]); - in_threshold_10_b = stof(vArgs[i+2]); + sphere_diameters_lower_bound = stof(vArgs[i+1]); + sphere_diameters_upper_bound = stof(vArgs[i+2]); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 4 numbers.\n"); + " argument must be followed by a number number.\n"); } num_arguments_deleted = 3; - } + } //if (vArgs[i] == "-spheres-nonmax-radii-range") - else if ((vArgs[i] == "-thresh-gauss") || - (vArgs[i] == "-thresh-gauss-out")) { + + else if ((vArgs[i] == "-spheres-nonmax-score-range") || + (vArgs[i] == "-sphere-nonmax-score-range")) { try { if (i+2 >= vArgs.size()) throw invalid_argument(""); - use_intensity_map = true; - use_gauss_thresholds = true; - out_thresh_gauss_x0 = stof(vArgs[i+1]); - out_thresh_gauss_sigma = stof(vArgs[i+2]); + score_lower_bound = stof(vArgs[i+1]); + score_upper_bound = stof(vArgs[i+2]); } catch (invalid_argument& exc) { throw InputErr("Error: The " + vArgs[i] + - " argument must be followed by 4 numbers.\n"); + " argument must be followed by two numbers.\n"); } num_arguments_deleted = 3; - } + } //if (vArgs[i] == "-spheres-nonmax-score-range") + else if (vArgs[i] == "-bs") { diff --git a/bin/filter_mrc/settings.hpp b/bin/filter_mrc/settings.hpp index 7033c5b..0a91cc3 100644 --- a/bin/filter_mrc/settings.hpp +++ b/bin/filter_mrc/settings.hpp @@ -33,28 +33,30 @@ class Settings { // this program. (Some of the algorithms share the same parameters.) typedef enum eFilterType { - NONE, //(the user has not selected a filter, or is using thresholding) + NONE, //(The user has not selected a filter, or is using thresholding) + DILATION, // grayscale image dilation + EROSION, // grayscale image erosion + OPENING, // grayscale image opening + CLOSING, // grayscale image closing GAUSS, //3D Gaussian filter GGAUSS, //3D Generalized Gaussian filter - DOG, //3D Difference of Gaussians (arbitrary widths, - // fast, seperable Gaussians) - DOGG, //3D Generaliws Difference-of-Generalized-Gaussians with arbitrary + DOG, //3D Difference of Gaussians (arbitrary widths, fast) + DOGG, //3D Generalized Difference-of-Generalized-Gaussians with arbitrary // widths and arbitrary exponents, slow) DOGGXY, //2D Generalized Difference-of-Generalized-Gaussians with arbitrary - // widths andarbitrary exponents, slow) - // and a 1D Gaussian filter (in the Z direction) - LOG_DOG, // Approximation to Laplacian-of-Guassian using DoG + // widths and arbitrary exponents in XY, and Gaussian blur in Z. + LOG_DOG, // Approximation to Laplacian-of-Guassian using DOG (fast) TEMPLATE_GAUSS, // Perform a least-squares fit between a Gaussian // and the surrounding voxel brightnesses TEMPLATE_GGAUSS, // Perform a least-squares fit between a generalized // Gaussian and the surrounding voxel brightnesses BOOTSTRAP_DOGG, // DOGG filter with bootstrapping (significance testing) - BLOB, // detect (scale-free) blobs - SURFACE_EDGE, // detect the edge of a light-dark boundary (gradient) - SURFACE_RIDGE, // detect 2D surface-like ridges (a thin membrane detector) - CURVE, // detect 1D curve-like ridges (a filament detector) WATERSHED, // Watershed segmentation - CLUSTER_CONNECTED, // Similar to Watershed segmentation + BLOB, // Scale-free blob detection + SURFACE_EDGE, // Detect the edge of a light-dark boundary (gradient) + SURFACE_RIDGE, // Detect 2D surface-like ridges (a thin membrane detector) + CURVE, // Detect 1D curve-like ridges (a filament detector) + CLUSTER_CONNECTED, // Watershed applied to voxels with tensor attributes LOCAL_FLUCTUATIONS, // Report the fluctuation of nearby voxel intensities SPHERE_NONMAX_SUPPRESSION,//throw away overlapping objects(detected earlier) SPHERE_NONMAX_SUPERVISED_MULTI,//use training data to throw away bad scoring blobs in multiple different images @@ -94,6 +96,7 @@ class Settings { string save_intermediate_fname_base = ""; //save progress of slow calculations string load_intermediate_fname_base = ""; //load progress from a previous run? + float thickness_morphology; // ---- parameters for "difference of generalized gaussians" filters ---- // // A significant number of the filters used by this program diff --git a/doc/doc_filter_mrc.md b/doc/doc_filter_mrc.md index cec2eec..928b14c 100644 --- a/doc/doc_filter_mrc.md +++ b/doc/doc_filter_mrc.md @@ -6,7 +6,8 @@ filter_mrc *(using [3D tensor voting](https://www.ncbi.nlm.nih.gov/pubmed/24625523))*. It can also detect point-like (or sphere-like) **blobs**. *(* **WARNING:** *The detection of curves is experimental as of 2021-6-21.)* - +*(Note: In this documentation, I use the words "image" and "tomogram" +interchangeably. A tomogram is a 3-D image.)* **filter_mrc** can be used local minima-finding, clustering, and [classic watershed segmentation](https://imagej.net/Classic_Watershed). @@ -19,6 +20,7 @@ Several primitive filters are included including low-pass, high-pass, edge detectors, ridge detectors, +morphological filters, thresholding, brightness inversions, [generalized](https://en.wikipedia.org/wiki/Generalized_normal_distribution#Version_1) @@ -38,7 +40,7 @@ The contributions from remaining voxels are normalized, so that objects located within narrow confined spaces can be detected accurately and without penalty.) Masks can also be used to give some voxels more consideration than others during the blurring (filtering) process. (A.K.A. "weighting".) -You can use a mask to to apply a filter to an image +In this way, you can use a mask to to apply a filter to an image whose boundaries are smooth and gradual as opposed to jagged and rectangular, @@ -211,9 +213,9 @@ too close to the surface. The following command will contract the boundary of the surface by approximately 40 Angstroms inward, creating a new file "segmented_interior.rec". ``` -filter_mrc -i segmented.rec -o segmented_interior.rec -erode-gauss 40 +filter_mrc -i segmented.rec -o segmented_interior.rec -erosion-gauss 40 ``` -*(See [-erode-gauss](#-erode-gauss-thickness) for details.)* +*(See [-erosion-gauss](#-erosion-gauss-thickness) for details.)* **Note:** @@ -444,28 +446,36 @@ structures which are bright on a dark background. The *thickness* parameter should be a number indicating the approximate target *thickness* for the thin membrane-like feature you are interested in detecting. -(This parameter should be expressed in physical units, not voxels. - Membranes whose thickness within a factor of 2 of this target - are also likely to be detected. - Technically, the *thickness* parameter controls the width Gaussian that will - be initially convolved with the original image, according to σ=thickness/√3. - For details, see Steger IEEE Trans. Pattern Anal. Machine Int. 1998.) +(This parameter should be expressed in physical units, not voxels. +Membranes whose thickness within a factor of 2 of this target +are also likely to be detected. +*Technically, the "thickness" parameter controls the width Gaussian, σ, +that will be initially convolved with the original image, according to +σ=thickness/√3. For details, see +[Steger IEEE Trans. Pattern Anal. Machine Int. 1998](https://doi.org/10.1109/34.659930).*) The *output* of this filter will be bright wherever the derivative of the brightness varies much more rapidly in one direction than it does in the two orthogonal directions. -Because this filter depends on second derivatives, -it is prone to detect a large number of spurious fluctuations in brightness -(in addition to the membrane-like structures you are interested in). -Tensor-voting (using the [-tv](#-tv-σ_ratio) argument) +Because this filter depends on second derivatives, it is prone to also +detect a large number of unwanted spurious fluctuations in brightness. +Tensor-voting (using the [**-tv**](#-tv-σ_ratio) argument) can be used to remove these spurious detected peaks -and improve the signal-to-noise ratio. +and greatly improve the signal-to-noise ratio. + +Groups of connected voxels that belong to the same surface can be +identified using the [**-connect**](#-connect-threshold) argument, +and their locations can be saved to a simple text file for geometric +analysis using the +[**-surface-normals-file**](#-surface-normals-file-PLY_FILE) argument. **WARNING** *This filter requires a very large amount of memory (enough to store at least 9 copies of the original image in 32bit-float mode).* + + #### Detection and masking You may wish to ignore certain voxels which are not part of the region @@ -490,8 +500,9 @@ noise near this boundary instead of the surface or curve that you are seeking. Alternatively, to get around this problem you can dilate (expand) the size of the mask region to include these voxels near the periphery. *(One way to do this is using the -[-dilate-gauss](#-dilate-gauss-thickness) -argument with a large thickness parameter. This is discussed below. +[-dilation](#-dilation-thickness) or +[-dilation-gauss](#-dilation-gauss-thickness) +arguments with a large thickness parameter. This is discussed below. Alternatively, if your only goal is to reduce the computation time, you can use the [-mask-rect](#-mask-rect-xmin-xmax-ymin-ymax-zmin-zmax) and [-mask-sphere](#-mask-sphere-x0-y0-z0-r) arguments to create a mask @@ -502,8 +513,10 @@ It may also help to vary the *size of the region* in the mask file using [dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)). (If you are using the "-mask" argument to load a mask image file, you can do this by creating a new mask file using *filter_mrc* again with the -[-erode-gauss](#-erode-gauss-thickness) or -[-dilate-gauss](#-dilate-gauss-thickness) arguments.) +[-erosion](#-erosion-thickness), +[-dilation](#-dilation-thickness), +[-erosion-gauss](#-erosion-gauss-thickness), or +[-dilation-gauss](#-dilation-gauss-thickness) arguments.) Eroding the mask can elimate traces of any structure at the edge of the mask boundary that might otherwise confuse or distract the detector. But, as mentioned earlier, dilating the mask is frequently more effective @@ -536,7 +549,7 @@ that you want to detect. As with **-surface** argument, the edge detection fidelity can be improved using tensor voting (with the [-tv](#-tv-σ_ratio) argument). -(Unless the "-w 1" argument was specified, then +(Unless the ["-w 1"](#Voxel-Width) argument was specified, then the *thickness* parameter should be in units of physical distance, such as Angstroms, not in voxels.) @@ -576,7 +589,7 @@ of the curve (in physical units). *As with surface detection, the curve detection fidelity can be improved using tensor voting (with the [-tv](#-tv-σ_ratio) argument).* -(Unless the "-w 1" argument was specified, +(Unless the ["-w 1"](#Voxel-Width) argument was specified, the *thickness* parameter should be in units of physical distance, such as Angstroms, not in voxels.) Also see the documentation for the [-surface](#-surface-type-thickness) @@ -681,360 +694,321 @@ The *fraction* parameter should lie in the range from 0 to 1. with using lower values. If the resulting surfaces are incomplete, it may help to increase this number and try again.) -### -connect threshold - -The -connect argument is a simple voxel clustering tool -used to distinguish different bright objects in the same image. -This is a two-step process. -First, all the voxels whose brightness exceeds *threshold* are selected. -Then, this subset of voxels is partitioned into different "islands". -"Islands" are defined as sets of voxels which are all -physically adjacent to (touching) eachother. -(The "bucket-fill" tool in a paint program behaves - in a similar way to the way islands are chosen. - The definition of "adjacency" can be controlled using the - [-neighbor-connectivity](#-neighbor-connectivity-nc) - argument.) +## -save-progress FILENAME +## -load-progress FILENAME -Once membership of each island has been decided, -a new image is generated showing which -voxels belong to each island. -(Note: This behavior is identical to the behavior of the - ["*-watershed maxima*"](#-watershed-type) - argument when used together with the - ["*-watershed-threshold*"](#-watershed-threshold-threshold) - argument.) +These two optional arguments work together to minimize the amount of time users +must spend redundantly waiting for slow operations (tensor voting) to complete. +The "-save-progress" argument will create 6 temporary files +(whose name begins with FILENAME and +ends with "_tensor_1.rec", ..., "_tensor_6.rec"). +This is useful because the process of stitching together a closed +membrane surface (or a continuous filamentous curve) +typically involves many iterations of trial and error, +using different thresholds and constraints until the +resulting surface (or curve) is reasonable. +If we use "-save-progress", then we can skip the time consuming +process of running tensor-voting each time. +We can do this using the "-load-progress FILENAME" argument. +The *FILENAME* argument should match the argument +you used with "-save-progress" earlier. +Be warned that files generated by "-save-progress" will consume +a substantial amount of disk space. You should delete them eventually. -*If the -["-connect"](#-connect-threshold) -argument is used together with the -["-surface"](#Detecting-membranes) -argument,* -(which is typically used for membrane detection), then it means that additional, -*more stringent* criteria will be used to distinguish nearby thin, curved -membrane-like objects from each other. -In particular, surfaces are assumed to be moderately smooth. -This means that adjacent voxels with radically different orientations -will never be grouped together. -Only voxels with similar orientations will be grouped into connected surfaces. -(The degree of similarity can be set by the user.) -In this case, the "*threshold*" parameter determines how *membrane-like* -a voxel must be in order for it to be included. -If the *threshold* parameter is chosen carefully, then these -different islands will hopefully correspond to different membranes -in the original image. -This *threshold* parameter will vary from image to image -and -[must be chosen carefully](#determining-the--connect-threshold-parameter). -If the *threshold* parameter is too large, -then individual objects (eg. membranes) in the image -will be split into multiple pieces. -If too small, then separate objects in the image will be joined together. +Note that when using "-load-progress", you must again specify +all the same arguments that you used earlier when you used +"-save-progress" +(including the "-in" "-surface", "-edge", and "-curve", and "-tv" arguments). +When you do this, you are permitted to make some changes which are listed below. +When you use "-load-progress", you are permitted to change the "-out" argument. +You can also add an image mask (or change the mask, for example +using the "-mask", "-mask-rect", "-mask-sphere" arguments). +Adding a mask is useful if you want to restrict the detection +of surface features to certain regions inside the image, for example. +(It makes sense to start using a mask now because +[masks are not typically used during the initial stages of detection](#Surface-detection-and-masking).) -Because it often takes several iterations to choose the correct -thresholds, it is recommended that you run *filter_mrc* once in advance -to detect the membrane, saving your progress using the -["-save-progress"](#-save-progress-FILENAME) -argument. *Then* when you are ready to connect the surfaces (or curves) -together using the "-connect" argument, use the -["-load-progress"](#-load-progress-FILENAME) -argument to load the directional features of the image that you measured earlier -(to avoid having to recalculate them again). -(This was demonstrated in [example 1](#Example-1) and [example 3](#Example-3).) -*Note:* If you are unable to find thresholds which connect all of -the pieces together correctly, you can also use the -["**-must-link**"](#-must-link-FILE) -argument. -This will manually force different bright regions in the image -to belong to the same cluster (a.k.a. "island". See below.) -#### Determining the -connect threshold parameter +## Blob detection -To choose the *threshold* parameter, -run membrane-detection -(for example using -["-surface"](#Detecting-membranes) -and -["-tv"](#-tv-σ_ratio) -) -once in advance without the -["-connect"](#-connect-threshold) -argument -with the -["-save-progress"](#-save-progress-FILENAME) -argument -(as we did in the membrane-detection -[example](#Example-1)). -Open the file created during that step -(eg. "membranes.rec") in IMOD. -Find a place in this image where the saliency (brightness) -of the membrane you are interested in is weak. -Click on voxels located near the weakest point (a.k.a. "junction point", -or "saddle point") between two different bright blobs -corresponding to the *same* surface you are interested in. -These two islands will not be joined unless the *-connect* argument -is less than the weakest link connecting these two islands. -(and even then, they might not be joined - if the voxel orientations are dissimilar.) -Select "Edit"->"Point"->"Value" menu option in IMOD to -see the "saliency" (brightness) of that voxel. -Do this several times in different places near the junction -write down the largest "saliency" number. -Then reduce this number by 20% (ie. multiply it by 0.8). -This makes a good first guess for the "*-connect*" parameter. -After using the "*-connect*" argument you can can -open the REC/MRC file we created -(eg "*membranes_clusters.rec*") -in IMOD, and click on different voxels (using "Edit"->"Point"->"Value") -to see whether the voxels were clustered correctly into the same object. -The voxel brightness values in that image should be integers -indicating which cluster they belong to -(reverse-sorted by cluster size, starting at 1). +### -blob, -blob-r, -blob-s -If some clusters are too big, you can either increase the *threshold* -value, *or* you can alter increase angular sensitivity by decreasing -the [*-connect-angle*](#-connect-angle-theta) from 45 degrees to 25 degrees -(or, equivalently by increasing the -*-cts*,*-ctn*,*-cvn*, and *-cvs* parameters from 0.707 to, say 0.9). +The "**-blob**", +["**-blob-r**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest), +and +["**-blob-s**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) +arguments are used for +[Scale-Free Blob-Detection](https://en.wikipedia.org/wiki/Blob_detection). +When this is used, the program will apply a LoG filter to the image +multiple times using a range of Gaussian widths (σ) (specified by the user) +in an attempt to determine the correct size (scale) for the relevant objects +in the image automatically. (As such, this operation is comparatively slow.) +A list of local minima and maxima in *X,Y,Z* space (and scale-space) +will generated and saved in a file, using the method described in: +Lindeberg,T., Int. J. Comput. Vision., 30(2):77-116, (1998) -Because it will probably take several tries to choose these parameters, -you can use the -["-load-progress"](#-load-progress-FILENAME) -argument to avoid having to recalculate the directional features -of the image that you (hopefully) measured earlier. -This was demonstrated in [example 3](#Example-3). +The "**-blob**", +["**-blob-r**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest), +and +["**-blob-s**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) +filters are followed +by 5 arguments (whose meaning depends on the filter selected): +``` + -blob type file_name d_min d_max gratio + -blob-r type file_name r_min r_max gratio + -blob-s type file_name σ_min σ_max gratio +``` +The "**type**" argument must be either "*minima*", "*maxima*", or "*all*". +*(It's meaning is explained in detail below. + For nearly all CryoEM images, it is safe to use "minima".)* -Also: It is a bad idea to try this on the original full-sized tomogram image. -1) Instead try this on one or more small cropped versions of the image -near the junction points of interest. -(You can crop images either by using IMOD, - or by using the "crop_mrc" program distributed with *visfd*.) -2) Again, you can also reduce the size and resolution of the image -(during the detection process), using the ["-bin"](#-bin-binsize) argument. -For example, using "-bin 2" will often also produce reasonable results -and will reduce the computation time considerably. Any features -detected at reduced resolution will be scaled up in size accordingly -(along with the voxel width). +The "**file_name**" argument is the name of a text file which will store +a list of detected blobs. The format of this file is explained below. +If "**-blob**" is selected, then the **d_min** and **d_max** parameters +(3rd and 4th parameters), specify a range of diameters +of the objects that you wish to detect. +(Simlarly, +["-blob-r"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) +allows the user to specify blob sizes in terms of their radii.) +Either way, a LoG filter will be applied to the image using a series of +different Gaussians whose widths (σ) vary between "**σ_min**", and "**σ_max**". +In this series, each successive Gaussian width (σ) will be larger than the +previous one by a factor of (no larger than) "**gratio**", +a number which should be > 1. +(**1.01** is a safe choice for **gratio**, + but you can speed up the calculation by increasing this parameter. + Values as high as 1.1 are not uncommon.) +(Note that: *σ_min*, and *σ_max* are equal to *d_min/(2√3)* and *d_max/(2√3)*, + or *r_min/√3* and *r_max/√3*, respectively. + If you prefer, you can use the + ["*-blob-s*"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) + argument to directly specify the + range of Gaussian widths you wish to use"*σ_min*", and "*σ_max*".) +After applying the LoG filter, if a voxel in the image is a (scale-adjusted) +local maxima (or minima) in *x,y,z,* and *σ*, +*then* its location and size will be recorded in a file. +"**file_name**" is the name of a file which will store these +locations of all the blobs detected in the image +as well as their sizes and scores (see below). +These detected blobs are either local minima or maxima in +X,Y,Z,[scale-space](https://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach). +Each file is a 5-column ascii text file with the following format: +``` +x1 y1 z1 diameter1 score1 +x2 y2 z2 diameter2 score2 + : : : : : +xM yM zM diameterM scoreM +``` +...where +"**M**" is the number of blobs (maxima or minima) which were detected, +On each line of that file, +**x,y,z** are the coordinates for the blob's center, +**diameter** is the diameter of the blob (which equals (2√3)σ), +and **score** is the intensity of that voxel after +a (scale-adjusted) LoG filter of that size was applied to it. +For *minima* type blobs, the list is ordered +from low score (most negative) to high score score +For *maxima* type blobs, the list is ordered from high score to low score. +These blobs can be visualized using the "**-draw-spheres**" argument +(see below). -Make sure clustering was successful *before* attempting to -close holes in the surface using programs like *meshlab* or *PoissonRecon*. +#### Blob types +The first argument indicates the **type** blob that the user wants to detect. +(The **type** is the 1st argument following the "-blob" argument.) +It must be one of the following choices: -### -must-link FILE +| type | meaning | +| -------- | ------- | +| *minima* | detect dark blobs on a bright background | +| *maxima* | detect bright blobs on a darker background | +| *all* | detect both dark and bright blobs | -*WARNING: This is an experimental feature. Please report bugs.* +Note: If **type** is set to *all*, then two different files will be generated +whose names will end in ".minima.txt" and ".maxima.txt", respectively. -If the "**-connect**" argument fails, you can manually force -different regions in the image to the same object -by specifying a text file containing the locations of pairs of voxels -that you wish to link together. You must prepare this text file in advance. -The coordinates of each voxel in the same connected group are listed -on separate lines in the file (which is in ASCII format). -Blank lines are used to separate voxels belonging to different objects. -There are two different file formats supported. +#### Blob detection and masking +It is recommended that you refrain from specifying a mask +*during* the process of blob detection +(for example, by using the +["**-mask**"](#-mask-MRC_FILE) +argument). +This is because the detector is more sensitive to noise in the image +near the mask boundary (or image boundary). +Instead, it is recommended that you use the +"-mask" argument *later on, after you have finished blob detection*. +(This was shown in [example 2](#Example-2) above.) -You can define links for many different objects in the image -in a single file by using blank-lines between different objects. -The following examples describe two *different* connected objects. -(For example, two different membranes.) - -#### Example 1: *must-link* coordinates in units of voxels (from IMOD) - -*If* the coordinates are enclosed in round-parenthesis () -and separated by spaces, (as shown below) -*then* it is assumed that these coordinates are in units of voxels, -and should range from 1 to Nx, Ny, or Nz (which define the image size). - -Considering the following text file: -``` -(16.5 77 73) -(134 75 73) - -(141.083 83.0833 62) -(123.833 133.333 64) -``` -*(Note: You* ***must include spaces*** between each coordinate, - even if commas are also present.)* - -This example describes two *different* connected objects -(separated by a blank line). -For each of these objects, all of the voxels connected to the -first voxel (for example, near position 16.5, 77, 73) -will be joined with all of the voxels connected to the second voxel -(near position 134, 75, 73). -The voxels in the last two lines of the file -will be connected together as well. - -*(Note: When using this format, coordinates are assumed to begin at 1, - and are rounded down to the next lowest integer. - These coordinates do not have to lie exactly on the object you - are trying to segment. Nearby voxels should work adequately well.)* - -**WARNING:** *Do not forget to put spaces between each integer (not commas).* - -#### Suggestion: *Use IMOD (3dmod)* - -The text above happens to be the text format printed by IMOD. -If you view a tomogram using "*3dmod -S*", -left-click anywhere on the image and press the "**F**" key, -the coordinates where you clicked will be printed to the -*3dmod* window in this format. -(The yellow "+" cursor, if visible, denotes this location.) -You can click on two places in the image that you want to link together. -Then copy the two lines printed by IMOD into your text file. -You can use this text file with the "**-must-link**" argument. - -For convenience, you can copy the entire line of text printed by IMOD -into your text file. (This includes the commas and the -extra text at the beginning and end of each line: -such as "Pixel" and "=...". This text will be ignored.) - -**WARNING:** *Remember to put blank lines between voxels that you *don't* - want to join together. - (And be warned that, depending on your thresholds, - they may get joined anyway.)* - - -#### Example 2: *must-link* coordinates in units of *physical distance* - -If no parenthesis are used, then it is assumed -that the coordinates are in physical distance units (eg Angstroms). -``` -379.68 1923.71 1822.46 -3366.5 1873.09 1822.46 - -3543.68 2075.58 1544.03 -3088.06 3341.18 1594.66 -``` -As with the previous example, -this example also describes two *different* connected objects. -For each of these objects, all of the voxels connected to the -first voxel (for example, at position 379.68,1923.71,1822.46) -will be joined with all of the voxels connected to the second voxel -(at position 3366.5,1873.09,1822.46). - -*(Note: When using this format, voxel positions begin at 0, not 1.)* - +Sometimes you might want to use the "-mask" argument during blob detection +because it can accelerate the search for blobs (by ignoring voxels outside +the mask). However don't do this if you care about blobs near the periphery. +Alternatively, to get around this problem you can dilate (expand) the size +of the mask region to include these voxels near the periphery. +*(One way to do this is using the +[-dilation](#-dilation-thickness) or +[-dilation-gauss](#-dilation-gauss-thickness) arguments.)* -### -select-cluster cluster-ID -### -surface-normals-file PLY_FILE +#### Non-max suppression: automatic disposal of blobs -Once clustering is working, you can select one of the clusters using -the "**-select-cluster**" argument. -(Cluster-IDs are assigned in reverse order according to their size, - beginning with the largest cluster, which has ID *1*.) -You can create a file which contains a list of voxels locations -(for the voxels belonging to that cluster), -as well as their surface orientations -using the "**-surface-normals-file**". -The resulting file will be in .PLY format. -This oriented point-cloud file can be used for further processing -(for example for hole-repair using the "PoissonRecon" program). +***By default, all minima and maxima are +reported during blob detection +(no matter how insignificant).*** +Unfortunately, the blob-detector will typically +an enormous number of blobs in a typical image. +(Running a blob detector on some Electron-Cryotomography + images can result in tens of millions of spurious blobs.) +The *vast majority* of these detected blobs are either +not relevant due to noise in the source image, or are overlapping. +You can discard these blobs *during* the blob detection process +using the arguments below. +However, because blob-detection is slow, +it is recommended that you save all of these blobs to a file. +*Then* you can decide which of these blobs are meaningful +by running non-max suppression later on (as shown below). +Then you can run "filter_mrc" a third time using the +["**-draw-spheres**"](#-draw-spheres-filename) +argument to visualize the remaining blobs you did not throw away. +(This is what we did in [example 2](#Example-2).) +Sometimes several iterations of non-max suppression followed by visualization +are needed before you achieve visually pleasing results. -Note: Voxel coordinates in the point-cloud file are expressed -in physical units (ie Angstroms), not voxels, -Consequently they are not integers (unless the -["-w 1"](#Voxel-Width) -argument was used). +Details are provided below explaining how visualize these blobs +(using the +[**-draw-spheres**](#-draw-spheres-filename) + argument), +as well as how to discard low-quality and overlapping blobs +(using the +[**-discard-blobs**](#-discard-blobs), +[**-radial-separation**](#Automatic-disposal-of-overlapping-blobs), +[**-max-volume-overlap**](#Automatic-disposal-of-overlapping-blobs), +and +[**-max-volume-overlap-small**](#Automatic-disposal-of-overlapping-blobs) +arguments). +#### Recommendation: -### -connect-angle theta +Blob-detection is computationally expensive, +but it only has to be performed once. +Typically, users will perform blob detection with the default (permissive) +score thresholds, so that they keep all of the blobs detected +(no matter how small or insignificant). +Then later, they will run **filter_mrc** again using *either* the +[**-discard-blobs**](#-discard-blobs), +["*-auto-thresh*"](#-auto-thresh-score--supervised-file_accepttxt-file_rejecttxt) +argument, *or* the +[**-minima-threshold**](#Automatic-disposal-of-poor-scoring-blobs) +(or +[**-maxima-threshold**](#Automatic-disposal-of-poor-scoring-blobs) +[**-draw-spheres**](#-draw-spheres-filename) +arguments, perhaps several times with different thresholds) +...to decide which of these blobs are worth keeping. -The *theta* argument determines how similar the -orientations of each voxel must be in order for a pair of neighboring voxels -to be merged into the same cluster (ie. the same membrane or same filament). -Neighboring voxels whose orientations differ by more than *theta* will -not be grouped into the same cluster (ie. the same surface or same curve). -If not specified, the default value is assumed to be 45 degrees. -However *theta* values from 15-60 degrees are common. -Note: The *-connect-angle* argument has no effect unless the -["-surface" and "-connect"](#Detecting-membranes) arguments are also supplied. -#### -cts, -ctn, -cvs, -cvn +# Morphology -*(Warning: This is an advanced experimental feature. -These arguments can probably be ignored by most users.)* +### -erosion thickness +### -dilation thickness -Alternatively, if you need more detailed control over the criteria -used to decide whether to add voxels to an existing cluster, -then you can specify 4 different thresholds using the -"**-cts threshold**", "**-ctn threshold**", -"**-cvs threshold**", and "**-cvn *threshold**" arguments. -This is an *alternative* to (specifying the "**-connect-angle**" argument). +This will apply a *grayscale* +[image dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)) and +[image erosion](https://en.wikipedia.org/wiki/Erosion_(morphology)) +filter to your image. This will enlarge the +[bright](https://en.wikipedia.org/wiki/Dilation_(morphology)#Grayscale_dilation) +or +[dark](https://en.wikipedia.org/wiki/Erosion_(morphology)#Grayscale_erosion) +regions in the image appear by a distance of *thickness*, respectively +(which is in physical units, not voxels, unless the +[-w 1](#Voxel-Width) argument is used). +*(Details: A spherical structure-factor with radius=thickness is used.)* + +Image dilation and erosion are useful for erasing or emphasizing +features (or noise) in the image that are in a given size range. +Dilation and erosion are also very useful to modify an image that +you eventually intend to use with the "-mask" argument +(in a future invocation of filter_mrc). +Recall that a *mask image* is a 3-D image used with the +["**-mask**"](#-mask-MRC_FILE) argument. +The brightness value of each voxel in a mask image is typically either 0 or 1 +*(depending on whether the voxel is supposed to be +ignored or considered during subsequent calculations)*. -Motivation and explanation: -Ideally, if tensor-voting worked, then the results of tensor voting should -agree with the actual direction of the features located at each voxel location. -If not, then you can exclude these voxels from consideration. -To do that, the "*-cts*" and "*-cvs*" threshold arguments specify -the minimum similarity between the hessian of the image brightness with the -tensor and vector features that resulted from tensor voting, respectively. -Additionally, the tensors from *neighboring* voxels which are part of the -same surface or curve should be poinging in the same direction. -So the "*-ctn*" and "*-cvn*" threshold arguments specify is the -minimum similarity between the orientations of the tensors and vectors -from neighboring voxels, respectively. -All of these threshold parameters should be in the range from -1 to 1. -A *-cvn -threshold* value of 0.707 ≈ cos(45°) and corresponds to a -45 degree difference between orientations of neighboring voxels. -In that case, neighboring voxels pointing in directions which -differ by more than 45° will be assigned to different clusters (membranes). -(This works reasonably well for phospholipid membranes.) +***WARNING:*** Dilation and erosion are slow for large *thickness* parameters. +If you have a binary image (images whose voxels have brightnesses of +either 0 or 1), you can use "-dilation-binary" or "-erosion-binary" instead. +*(<--NOTE: THIS FEATURE IS NOT IMPLEMENTED YET. -A 2021-7-11.)* +For binary images, you can also use the Gaussian dilation and erosion +method, which is even faster, using the +[-dilation-gauss](#-dilation-gauss-thickness) and +[-erosion-gauss](#-erosion-gauss-thickness) +arguments. +*(For additional morphological filter operations, see +[dilate](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.dilation) +[erode](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.erosion) +[binary_dilate](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.binary_dilation) +[binary_erode](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.binary_erosion) +from +[scikit-image](https://scikit-image.org), +along with the +[mrcfile](https://mrcfile.readthedocs.io/en/latest/readme.html#basic-usage) +module.)* -Note that the *-cts*, *-ctn*, *-cvs*, and *-cvn* arguments -have no effect unless the ["-surface" and "-connect"](#Detecting-membranes) -arguments are also in use. +### -opening thickness +### -closing thickness +This will apply a *grayscale* +[image opening](https://en.wikipedia.org/wiki/Opening_(morphology)) and +[image closing](https://en.wikipedia.org/wiki/Closing_(morphology)) +filter to your image. +The morphological opening of an image is defined as an +[erosion followed by a dilation](https://scikit-image.org/docs/dev/api/skimage.morphology.html#skimage.morphology.opening). +Opening removes small objects, while closing removes small holes. +The size of the objects removed is determined by the *thickness* parameter +(which is in physical units, not voxels, unless the +[-w 1](#Voxel-Width) argument is used). +*(Details: A spherical structure-factor with radius=thickness is used.)* -## -save-progress FILENAME -## -load-progress FILENAME +### -erosion-gauss thickness +### -dilation-gauss thickness -These two optional arguments work together to minimize the amount of time users -must spend redundantly waiting for slow operations (tensor voting) to complete. -The "-save-progress" argument will create 6 temporary files -(whose name begins with FILENAME and -ends with "_tensor_1.rec", ..., "_tensor_6.rec"). -This is useful because the process of stitching together a closed -membrane surface (or a continuous filamentous curve) -typically involves many iterations of trial and error, -using different thresholds and constraints until the -resulting surface (or curve) is reasonable. -If we use "-save-progress", then we can skip the time consuming -process of running tensor-voting each time. -We can do this using the "-load-progress FILENAME" argument. -The *FILENAME* argument should match the argument -you used with "-save-progress" earlier. -Be warned that files generated by "-save-progress" will consume -a substantial amount of disk space. You should delete them eventually. +This is a fast and crude implementation of +[image dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)) and +[image erosion](https://en.wikipedia.org/wiki/Erosion_(morphology)) +that uses fast Gaussian filters to blur an image, followed by thresholding. +The **-dilation-gauss** and **-erosion-gauss** arguments will enlarge +or shrink the size of the region stored in a binary image by +a distance of roughly *thickness* (which has units of physical distance +unless the ["-w 1"](#Voxel-Width) argument is used). +Note that, due to blurring, features in the image which are smaller or +narrower than *thickness* will be lost after this operation is completed. +#### Details +This is implemented by blurring the image, and then applying a threshold +so that all voxels whose brightness are below a threshold are discarded. +Consequently, you cannot combine these arguments with either the *-gauss* +or *-thresh* arguments, because: -Note that when using "-load-progress", you must again specify -all the same arguments that you used earlier when you used -"-save-progress" -(including the "-in" "-surface", "-edge", and "-curve", and "-tv" arguments). -When you do this, you are permitted to make some changes which are listed below. -When you use "-load-progress", you are permitted to change the "-out" argument. -You can also add an image mask (or change the mask, for example -using the "-mask", "-mask-rect", "-mask-sphere" arguments). -Adding a mask is useful if you want to restrict the detection -of surface features to certain regions inside the image, for example. -(It makes sense to start using a mask now because -[masks are not typically used during the initial stages of detection](#Surface-detection-and-masking).) +- "-dilation-gauss *thickness*" is equivalent to +"-gauss *thickness* -thresh 0.157299", *(where 0.157299 ≈ 1-erf(1))* +- "-erosion-gauss *thickness*" is equivalent to +"-gauss *thickness* -thresh 0.842701", *(where 0.842701 ≈ erf(1))* +These thresholds were chosen so that the selected region +expands or shrinks by a distance of *thickness*. +(But, technically, this is only true near a flat boundary.) @@ -1076,7 +1050,7 @@ You can use the to create an image showing which voxels belong to each minima or maxima. (In that image, the voxel brightness will be an integer indicating - the minima or maxima to which the voxel belongs, if any, starting at 1. + the minima or maxima to which the voxel belongs, if any, starting at 1. When both minima and maxima are displayed, the minima are represented by negative integers.) @@ -1085,16 +1059,6 @@ of the image (or on the boundary of the mask region), are considered. To ignore these extrema (ie, to consider only voxels which are surrounded by other voxels) use the "**-ignore-boundary-extrema**" argument. -##### -neighbor-connectivity nc - -*Note:* When searching neighboring voxels, all 26 neighbors (3x3x3-1) -surrounding the voxel are considered by default. -The "*nc*" parameter represents the maximum squared distance a voxel can -be from the center voxel in order for it to be considered a "neighbor". -(You can use the **-neighbor-connectivity nc** argument to skip - 3D corner voxels by setting nc=2, and also 2D corners by setting nc=1. - This may increase the number of spurious local minima and maxima discovered.) - #### Non-max suppression: automatic disposal of minima or maxima @@ -1128,190 +1092,7 @@ then the poorer scoring minima will be discarded. -## Blob detection - - -### -blob, -blob-r, -blob-s - -The "**-blob**", -["**-blob-r**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest), -and -["**-blob-s**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) -arguments are used for -[Scale-Free Blob-Detection](https://en.wikipedia.org/wiki/Blob_detection). -When this is used, the program will apply a LoG filter to the image -multiple times using a range of Gaussian widths (σ) (specified by the user) -in an attempt to determine the correct size (scale) for the relevant objects -in the image automatically. (As such, this operation is comparatively slow.) -A list of local minima and maxima in *X,Y,Z* space (and scale-space) -will generated and saved in a file, using the method described in: -Lindeberg,T., Int. J. Comput. Vision., 30(2):77-116, (1998) - -The "**-blob**", -["**-blob-r**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest), -and -["**-blob-s**"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) -filters are followed -by 5 arguments (whose meaning depends on the filter selected): -``` - -blob type file_name d_min d_max gratio - -blob-r type file_name r_min r_max gratio - -blob-s type file_name σ_min σ_max gratio -``` -The "**type**" argument must be either "*minima*", "*maxima*", or "*all*". -*(It's meaning is explained in detail below. - For nearly all CryoEM images, it is safe to use "minima".)* - -The "**file_name**" argument is the name of a text file which will store -a list of detected blobs. The format of this file is explained below. - -If "**-blob**" is selected, then the **d_min** and **d_max** parameters -(3rd and 4th parameters), specify a range of diameters -of the objects that you wish to detect. -(Simlarly, -["-blob-r"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) -allows the user to specify blob sizes in terms of their radii.) -Either way, a LoG filter will be applied to the image using a series of -different Gaussians whose widths (σ) vary between "**σ_min**", and "**σ_max**". -In this series, each successive Gaussian width (σ) will be larger than the -previous one by a factor of (no larger than) "**gratio**", -a number which should be > 1. -(**1.01** is a safe choice for **gratio**, - but you can speed up the calculation by increasing this parameter. - Values as high as 1.1 are not uncommon.) -(Note that: *σ_min*, and *σ_max* are equal to *d_min/(2√3)* and *d_max/(2√3)*, - or *r_min/√3* and *r_max/√3*, respectively. - If you prefer, you can use the - ["*-blob-s*"](#Specifying-the-radius-or-Gaussian-sigma-parameters-for-the-objects-of-interest) - argument to directly specify the - range of Gaussian widths you wish to use"*σ_min*", and "*σ_max*".) -After applying the LoG filter, if a voxel in the image is a (scale-adjusted) -local maxima (or minima) in *x,y,z,* and *σ*, -*then* its location and size will be recorded in a file. -"**file_name**" is the name of a file which will store these -locations of all the blobs detected in the image -as well as their sizes and scores (see below). -These detected blobs are either local minima or maxima in -X,Y,Z,[scale-space](https://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach). -Each file is a 5-column ascii text file with the following format: -``` -x1 y1 z1 diameter1 score1 -x2 y2 z2 diameter2 score2 - : : : : : -xM yM zM diameterM scoreM -``` -...where -"**M**" is the number of blobs (maxima or minima) which were detected, -On each line of that file, -**x,y,z** are the coordinates for the blob's center, -**diameter** is the diameter of the blob (which equals (2√3)σ), -and **score** is the intensity of that voxel after -a (scale-adjusted) LoG filter of that size was applied to it. -For *minima* type blobs, the list is ordered -from low score (most negative) to high score score -For *maxima* type blobs, the list is ordered from high score to low score. -These blobs can be visualized using the "**-draw-spheres**" argument -(see below). - -#### Blob types - -The first argument indicates the **type** blob that the user wants to detect. -(The **type** is the 1st argument following the "-blob" argument.) -It must be one of the following choices: - -| type | meaning | -| -------- | ------- | -| *minima* | detect dark blobs on a bright background | -| *maxima* | detect bright blobs on a darker background | -| *all* | detect both dark and bright blobs | - -Note: If **type** is set to *all*, then two different files will be generated -whose names will end in ".minima.txt" and ".maxima.txt", respectively. - - -#### Blob detection and masking - -It is recommended that you refrain from specifying a mask -*during* the process of blob detection -(for example, by using the -["**-mask**"](#-mask-MRC_FILE) -argument). -This is because the detector is more sensitive to noise in the image -near the mask boundary (or image boundary). -Instead, it is recommended that you use the -"-mask" argument *later on, after you have finished blob detection*. -(This was shown in [example 2](#Example-2) above.) - -Sometimes you might want to use the "-mask" argument during blob detection -because it can accelerate the search for blobs (by ignoring voxels outside -the mask). However don't do this if you care about blobs near the periphery. -Alternatively, to get around this problem you can dilate (expand) the size -of the mask region to include these voxels near the periphery. -*(One way to do this is using the -[-dilate-gauss](#-dilate-gauss-thickness) argument.)* - - -#### Non-max suppression: automatic disposal of blobs - -***By default, all minima and maxima are -reported during blob detection -(no matter how insignificant).*** -Unfortunately, the blob-detector will typically -an enormous number of blobs in a typical image. -(Running a blob detector on some Electron-Cryotomography - images can result in tens of millions of spurious blobs.) -The *vast majority* of these detected blobs are either -not relevant due to noise in the source image, or are overlapping. -You can discard these blobs *during* the blob detection process -using the arguments below. -However, because blob-detection is slow, -it is recommended that you save all of these blobs to a file. -*Then* you can decide which of these blobs are meaningful -by running non-max suppression later on (as shown below). -Then you can run "filter_mrc" a third time using the -["**-draw-spheres**"](#-draw-spheres-filename) -argument to visualize the remaining blobs you did not throw away. -(This is what we did in [example 2](#Example-2).) -Sometimes several iterations of non-max suppression followed by visualization -are needed before you achieve visually pleasing results. - -Details are provided below explaining how visualize these blobs -(using the -[**-draw-spheres**](#-draw-spheres-filename) - argument), -as well as how to discard low-quality and overlapping blobs -(using the -[**-discard-blobs**](#-discard-blobs), -[**-radial-separation**](#Automatic-disposal-of-overlapping-blobs), -[**-max-volume-overlap**](#Automatic-disposal-of-overlapping-blobs), -and -[**-max-volume-overlap-small**](#Automatic-disposal-of-overlapping-blobs) -arguments). - - - -#### Recommendation: - -Blob-detection is computationally expensive, -but it only has to be performed once. -Typically, users will perform blob detection with the default (permissive) -score thresholds, so that they keep all of the blobs detected -(no matter how small or insignificant). -Then later, they will run **filter_mrc** again using *either* the -[**-discard-blobs**](#-discard-blobs), -["*-auto-thresh*"](#-auto-thresh-score--supervised-file_accepttxt-file_rejecttxt) -argument, *or* the -[**-minima-threshold**](#Automatic-disposal-of-poor-scoring-blobs) -(or -[**-maxima-threshold**](#Automatic-disposal-of-poor-scoring-blobs) -[**-draw-spheres**](#-draw-spheres-filename) -arguments, perhaps several times with different thresholds) -...to decide which of these blobs are worth keeping. - - - - -## Segmentation +# Segmentation ### -watershed type @@ -1399,226 +1180,411 @@ When this argument is included, these voxels will have intensities which are assigned to *label* instead. (This parameter is a number.) +### -neighbor-connectivity nc +*Note:* When searching neighboring voxels, all 26 neighbors (3x3x3-1) +surrounding the voxel are considered by default. +The "*nc*" parameter represents the maximum squared distance a voxel can +be from the center voxel in order for it to be considered a "neighbor". +(You can use the **-neighbor-connectivity nc** argument to skip + 3D corner voxels by setting nc=2, and also 2D corners by setting nc=1. + This may increase the number of spurious local minima and maxima discovered.) -## General filters: -### -gauss and -ggauss -The **-gauss** and **-gauss-aniso** arguments must be followed by one or more numbers specifying the width of the Gaussian filter to apply to your image: -``` - -gauss σ -``` -or -``` - -gauss-aniso σ_x σ_y σ_z -``` -In either case, the original image is convolved with: -``` - h(x,y,z) = A*exp(-0.5 * ((x/σ_x)^2 + (y/σ_y)^2 + (z/σ_z)^2)) -``` -When *-gauss* is used, *σ_x = σ_y = σ_z = σ* -If **-ggauss** or **-ggauss-aniso** is used, -then the image is instead convolved with -[a generalized gaussian function](https://en.wikipedia.org/wiki/Generalized_normal_distribution#Version_1): -``` - h(x,y,z) = A*exp(-r^m) - where r = √((x/σ_x)^2 + (y/σ_y)^2 + (z/σ_z)^2)) -``` -The width of the Gaussian (ie, the σ_x, σ_y, σ_z arguments) should be specified in units of physical distance, not in voxels. -(The A coefficient will be determined automatically by normalization.) -By default **m** *= 2*, -however this can be overridden using the optional "**-exponent m**" -*Note:* that these *σ*, *σ_x*, *σ_y*, *σ_z*, parameters are a factor -of *√2 times larger* than the corresponding *σ* ("sigma") parameter -traditionally used in the -[Gaussian function](https://en.wikipedia.org/wiki/Normal_distribution). -(For your convenience, you can plot these functions for various parameters using -the "draw_filter1D.py -ggauss A s m" python script which is located in the same directory. -In this case, the "s" parameter is specified in voxels, not physical distance.) -*Note:* The calculation is -[fast](https://en.wikipedia.org/wiki/Separable_filter) -if you use the default exponent of 2. -*Changing the exponent will slow down the filter considerably.* +### -connect threshold -The filter is truncated far away from the central peak at a point which is chosen automatically according the σ, σ_x, σ_y, σ_z parameters selected by the user. However this can be customized using the -["-truncate-threshold"](#Filter-Size) -and -["-truncate"](#Filter-Size) -arguments if necessary (see below). +The -connect argument is a voxel clustering tool +used to distinguish different bright objects in the same image. +This argument will cluster nearby (adjacent) bright voxels +*(or voxels with high "saliency")* into separate islands and then +generate a new image showing *which* island each voxel belongs to (if any). -### -dog -When the "**-dog**" or "**-dog-aniso**" filter is selected, a -[Difference-of-Gaussians](https://en.wikipedia.org/wiki/Difference_of_Gaussians) -filter will be used. -The **-dog** and **-dog-aniso** arguments must be followed by numbers indicating the width of the Gaussians you wish to convolve with your image. -(Equivalently, the numbers indicate the band of frequencies you wish to keep from your original image.) -``` - -dog a b -``` -or: -``` - -dog-aniso a_x a_y a_z b_x b_y b_z -``` -In either case, -The original image will be convolved with the following function: -``` - h(x,y,z) = A*exp(-0.5 * ((x/a_x)^2 + (y/a_y)^2 + (z/a_z)^2)) - - B*exp(-0.5 * ((x/b_x)^2 + (y/b_y)^2 + (z/b_z)^2)) -``` -When *-dog* is used, -*a_x = a_y = a_z = a* -and -*b_x = b_y = b_z = b* +This is a two-step process. +First, all the voxels whose brightness exceeds *threshold* are selected. +Then, this subset of voxels is partitioned into different "islands". +"Islands" are defined as sets of voxels which are all +physically adjacent to (touching) eachother. +(The definition of "adjacency" can be controlled using the +[-neighbor-connectivity](#-neighbor-connectivity-nc) +argument.) -### LoG filters -### -log, -log-r, -log-d, and -log-aniso +Once membership of each island has been decided, +a new image is generated showing which +voxels belong to each island. +(Note: This behavior is identical to the behavior of the +["*-watershed maxima*"](#-watershed-type) +argument when used together with the +["*-watershed-threshold*"](#-watershed-threshold-threshold) +argument. However the "-connect" argument has additional options.) +*If the +["-connect"](#-connect-threshold) +argument is used together with the +["-surface"](#Detecting-membranes) +argument,* +(which is typically used for membrane detection), then it means that additional, +*more stringent* criteria will be used to distinguish nearby thin, curved +membrane-like objects from each other. +In that case, it is not the brightness of the voxel that matters, but its +*saliency* (a number indicating how closely the region near the voxel resembles +the feature you are trying to detect, like a curve, surface, or edge). +In order for two neighboring voxels to be in the same island, they must +both have high saliency above the *threshold* parameter specified by the user. +In this case, the "*threshold*" parameter determines how *membrane-like* +a voxel must be in order for it to be included. +If the *threshold* parameter is chosen carefully, then these +different islands will hopefully correspond to different membranes +in the original image. -When the "**-log**", "**-log-r**", and "**-log-d**", or "**-log-aniso**" -arguments are selected, a -[Laplacian-of-a-Gaussian (LoG)](https://en.wikipedia.org/wiki/Blob_detection) -filter is applied on the source image. -(See implementation details below.) -The Laplacian-of-a-Gaussian filter applies a Gaussian blur to the image -(whose width is specified by the user), -and then calculates the Laplacian of the result. -Features in the image of various sizes can be emphasized -using this kind of filter, and the result can be saved as a new image using -["**-out** filename.rec"](#Input-and-Output-files). -The "**-log**", "**-log-aniso**", "**-log-r**", and "**-log-d**" -arguments described here are typically only chosen when the user already -knows the approximate size of the features they want to emphasize -in the image. -``` - -log σ -``` -When "**-log σ**" is used, a LoG filter will be applied -to the image using a Gaussian width (σ). -Alternatively, if the user wishes to specify the actual size -(the radius or diameter) of the objects they want to emphasize, -then they can (equivalently) use -``` - -log-r radius -``` -or -``` - -log-d diameter -``` -instead. -There is also an anisotropic version of this filter as well: -``` - -log-aniso σ_x σ_y σ_z** -``` -The new image which is generated will *tend* to have have bright and dark spots -wherever bright and dark objects of the corresponding size can be found. -This image can be searched for (local) minima and maxima -by running this program again on the new image using the -"**-find-minima**", "**-find-maxima**" and "-radial-separation" arguments. -These locations can be saved to a file and then a new anotated image can be -created using the -["**-draw-spheres**"](#-draw-spheres-filename) -argument. This provides a fast, sloppy, -and unselective way to search for features in an image. -However scale-free -[**blob detection**](https://en.wikipedia.org/wiki/Blob_detection) -is a more robust (and computationally expensive) way to detect objects -of a given size within an image. +This *threshold* parameter will vary from image to image +and +[must be chosen carefully](#determining-the--connect-threshold-parameter). +If the *threshold* parameter is too large, +then individual objects (eg. membranes) in the image +will be split into multiple pieces. +If too small, then separate objects in the image will be joined together. -*(Implementation note: The LoG filter described here is approximated internally - using a DoG filter. You can control the accuracy of this approximation using - the - ["-dog-delta δ"](#Implementation-of-LoG-filters-and-blob-detectors) - argument.) +Because it often takes several iterations to choose the correct +thresholds, it is recommended that you run *filter_mrc* once in advance +to detect the membrane, saving your progress using the +["-save-progress"](#-save-progress-FILENAME) +argument. *Then* when you are ready to connect the surfaces (or curves) +together using the "-connect" argument, use the +["-load-progress"](#-load-progress-FILENAME) +argument to load the directional features of the image that you measured earlier +(to avoid having to recalculate them again). +(This was demonstrated in [example 1](#Example-1) and [example 3](#Example-3).) +*Note:* If you are unable to find thresholds which connect all of +the pieces together correctly, you can also use the +["**-must-link**"](#-must-link-FILE) +argument. +This will manually force different bright regions in the image +to belong to the same cluster (a.k.a. "island". See below.) +#### Clustering based on directional similarity + +Some detection methods (including surface, edge, and curve detectors) +also report *directional* features in the image (such as the curve +direction or the surface normal direction at each location in the image). +*(Tensor attributes, such as the second derivative (Hessian), +are also reported.)* +The clustering algorithm used here is aware of any directional +*(as well as tensor)* attributes of the voxels (if applicable). + +These directional attributes can be used to increase clustering sensitivity. +Directional information can be used distinguish different objects from +each other that are otherwise touching and might be grouped together. +For example, curves and surfaces are assumed to be moderately smooth. +You can prevent adjacent voxels with radically different +orientations from being grouped together (even if they both have high saliency). +This way, only voxels with similar orientations will be grouped into +connected surfaces. (The degree of similarity can be set by the user +using the [*-connect-angle*](#-connect-angle-theta) argument.) +This strategy was used in [example 3](#Example-3). -### -fluct and -fluct-aniso -Usage: +#### Determining the -connect threshold parameter + +Let's assume you are using the "-surface" argument to detect thin +membrane-like surfaces in an image. +To choose the *threshold* parameter, +run membrane-detection +(for example using +["-surface"](#Detecting-membranes) +and +["-tv"](#-tv-σ_ratio) +) +once in advance without the +["-connect"](#-connect-threshold) +argument +with the +["-save-progress"](#-save-progress-FILENAME) +argument +(as we did in the membrane-detection +[example](#Example-1)). +Open the file created during that step +(eg. "membranes.rec") in IMOD. +Find a place in this image where the saliency (brightness) +of the membrane you are interested in is weak. +Click on voxels located near the weakest point (a.k.a. "junction point", +or "saddle point") between two different bright blobs +corresponding to the *same* surface you are interested in. +These two islands will not be joined unless the *-connect* argument +is less than the weakest link connecting these two islands. +(and even then, they might not be joined + if the voxel orientations are dissimilar.) +Select "Edit"->"Point"->"Value" menu option in IMOD to +see the "saliency" (brightness) of that voxel. +Do this several times in different places near the junction +write down the largest "saliency" number. +Then reduce this number by 20% (ie. multiply it by 0.8). +This makes a good first guess for the "*-connect*" parameter. + +After using the "*-connect*" argument you can can +open the REC/MRC file we created +(eg "*membranes_clusters.rec*") +in IMOD, and click on different voxels (using "Edit"->"Point"->"Value") +to see whether the voxels were clustered correctly into the same object. +The voxel brightness values in that image should be integers +indicating which cluster they belong to +(reverse-sorted by cluster size, starting at 1). + +If some clusters are too big, you can either increase the *threshold* +value, *or* you can alter increase angular sensitivity by reducing the +[*-connect-angle*](#-connect-angle-theta) from 45 degrees to 15 degrees +(or, equivalently by increasing the +*-cts*,*-ctn*,*-cvn*, and *-cvs* thresholds from 0.707 to 0.966). + +Because it will probably take several tries to choose these parameters, +you can use the +["-load-progress"](#-load-progress-FILENAME) +argument to avoid having to recalculate the directional features +of the image that you (hopefully) measured earlier. +This was demonstrated in [example 3](#Example-3). + +Also: It is a bad idea to try this on the original full-sized tomogram image. +1) Instead try this on one or more small cropped versions of the image +near the junction points of interest. +(You can crop images either by using IMOD, + or by using the "crop_mrc" program distributed with *visfd*.) +2) Again, you can also reduce the size and resolution of the image +(during the detection process), using the ["-bin"](#-bin-binsize) argument. +For example, using "-bin 2" will often also produce reasonable results +and will reduce the computation time considerably. Any features +detected at reduced resolution will be scaled up in size accordingly +(along with the voxel width). + +Make sure clustering was successful *before* attempting to +close holes in the surface using programs like *meshlab* or *PoissonRecon*. + + +### -must-link FILE + +*WARNING: This is an experimental feature. Please report bugs.* + +If the "**-connect**" argument fails, you can manually force +different regions in the image to the same object +by specifying a text file containing the locations of pairs of voxels +that you wish to link together. You must prepare this text file in advance. + +The coordinates of each voxel in the same connected group are listed +on separate lines in the file (which is in ASCII format). +Blank lines are used to separate voxels belonging to different objects. +There are two different file formats supported. + + +You can define links for many different objects in the image +in a single file by using blank-lines between different objects. +The following examples describe two *different* connected objects. +(For example, two different membranes.) + +#### Example 1: *must-link* coordinates in units of voxels (from IMOD) + +*If* the coordinates are enclosed in round-parenthesis () +and separated by spaces, (as shown below) +*then* it is assumed that these coordinates are in units of voxels, +and should range from 1 to Nx, Ny, or Nz (which define the image size). + +Considering the following text file: ``` - fluct r +(16.5 77 73) +(134 75 73) + +(141.083 83.0833 62) +(123.833 133.333 64) ``` -The "**-fluct**" filter calculates the magnitude of -the voxel intensity fluctuations relative to their local surroundings. -This is useful for finding regions in an image with nearly uniform brightness -and regions whose brightness fluctuates wildly. -*r* is the effective search radius. +*(Note: You* ***must include spaces*** between each coordinate, +even if commas are also present.)* -To compute this, near every voxel, a Gaussian-weighted average -intensity is computed with the same effective radius, *r*. -A Gaussian-weighted squared-difference between all the nearby voxel intensities -and this local average intensity is also computed. -A new image is generated whose voxel intensities equal this -local average squared difference intensity of nearby voxels. +This example describes two *different* connected objects +(separated by a blank line). +For each of these objects, all of the voxels connected to the +first voxel (for example, near position 16.5, 77, 73) +will be joined with all of the voxels connected to the second voxel +(near position 134, 75, 73). +The voxels in the last two lines of the file +will be connected together as well. -The "**-fluct-aniso**" variant allows the user to control the width -of the Gaussian independently in the x,y,z directions: +*(Note: When using this format, coordinates are assumed to begin at 1, + and are rounded down to the next lowest integer. + These coordinates do not have to lie exactly on the object you + are trying to segment. Nearby voxels should work adequately well.)* + +**WARNING:** *Do not forget to put spaces between each integer (not commas).* + +#### Suggestion: *Use IMOD (3dmod)* + +The text above happens to be the text format printed by IMOD. +If you view a tomogram using "*3dmod -S*", +left-click anywhere on the image and press the "**F**" key, +the coordinates where you clicked will be printed to the +*3dmod* window in this format. +(The yellow "+" cursor, if visible, denotes this location.) +You can click on two places in the image that you want to link together. +Then copy the two lines printed by IMOD into your text file. +You can use this text file with the "**-must-link**" argument. + +For convenience, you can copy the entire line of text printed by IMOD +into your text file. (This includes the commas and the +extra text at the beginning and end of each line: +such as "Pixel" and "=...". This text will be ignored.) + +**WARNING:** +*Remember to put blank lines between voxels that you *don't* want to join +together. (And be warned that, depending on your thresholds, +they may get joined anyway.)* + + +#### Example 2: *must-link* coordinates in units of *physical distance* + +If no parenthesis are used, then it is assumed +that the coordinates are in physical distance units (eg Angstroms). ``` - fluct-aniso r_x r_y r_z +379.68 1923.71 1822.46 +3366.5 1873.09 1822.46 + +3543.68 2075.58 1544.03 +3088.06 3341.18 1594.66 ``` -#### Implementation detail: - The width of the Gaussian, σ, is chosen so that the effective - number of voxels under the 3D Gaussian, - (which can be interpreted as the 3D integral of the unnormalized 3D Gaussian, - (2\*π\*σ^2)^(3/2)) - equals the volume of a sphere of radius r - (4/3)\*π\*r^3. - This yields σ=(9π/2)^(-1/6)r. +As with the previous example, +this example also describes two *different* connected objects. +For each of these objects, all of the voxels connected to the +first voxel (for example, at position 379.68,1923.71,1822.46) +will be joined with all of the voxels connected to the second voxel +(at position 3366.5,1873.09,1822.46). - It might seem more straightforward to simply consider the fluctuations in - brightnesses of all of the voxels which lie within a fixed radius from - each target voxel. However using Gaussian weighted averages instead - accomplishes essentially the same goal and is much more - [computationally efficient](https://en.wikipedia.org/wiki/Separable_filter). +*(Note: When using this format, voxel positions begin at 0, not 1.)* - If you prefer, instead of using a Gaussian, you can instead perform the - fluctuation calculation over a hard spherical region by using the - "-exponent n" argument. - (This parameter is 2 by default. Changing it will slow down the calculation.) - As n is increased, the region over which the fluctuations in brightness - are calculated becomes more and more like a uniform sphere with radius σ. - (So if you do this, be sure to multiply your radius argument, r, - by (9π/2)^(1/6)≈1.5549880806696572 to compensate.) -## Resizing the image +### -select-cluster cluster-ID +### -surface-normals-file PLY_FILE -### -bin binsize +Once clustering is working, you can select one of the clusters using +the "**-select-cluster**" argument. +(Cluster-IDs are assigned in reverse order according to their size, + beginning with the largest cluster, which has ID *1*.) +You can create a file which contains a list of voxels locations +(for the voxels belonging to that cluster), +as well as their surface orientations +using the "**-surface-normals-file**". +The resulting file will be in .PLY format. +This oriented point-cloud file can be used for further processing +(for example for hole-repair using the "PoissonRecon" program). -Reduce the resolution of the image in each direction by a factor of *binsize*. -(*binsize* must be a positive integer.) -The new image will be smaller in each direction by a factof of *binsize*. +Note: Voxel coordinates in the point-cloud file are expressed +in physical units (ie Angstroms), not voxels, +Consequently they are not integers (unless the +["-w 1"](#Voxel-Width) argument was used). -Motivation: Reducing the resolution of an image can often dramatically -reduce the computation time necessary to detect objects in the image. -If you are detecting features in the image (such as membranes or blobs) -that are significantly larger (or thicker) than the voxel width, -then a reduction in resolution by a factor of 2 or so -should not effect your ability to detect it accurately, -(and could make tensor voting up to 64 times faster (as of 2021-6-21)). -Note that you can use this argument together with other arguments. -If you do that, the reduction of resolution occurs before all -other operations are performed (such as "-blob" or "-surface" detection). -This allows you to reduce the image size before -these other expensive calculations are performed. +### -connect-angle theta -Note that the width of each voxel will be increased accordingly, -so that the positions of any features in the image that you detect -should not be significantly effected by the change in resolution. -(This includes the coordinates of the surface mesh.) -Similarly, when using "-bin", the user does not need to alter -any of the other parameters or arguments passed to the program -(such as the thickness parameters or sigma parameters -that characterize the size of the objects they want to detect). -The software will automatically adjust these parameters -to compensate for the change in image resolution. +The *theta* argument determines how similar the +orientations of each voxel must be in order for a pair of neighboring voxels +to be merged into the same cluster (ie. the same membrane or same filament). +Neighboring voxels whose orientations differ by more than *theta* will +not be grouped into the same cluster (ie. the same surface or same curve). +If not specified, the default value is assumed to be 45 degrees. +However *theta* values from 15-60 degrees are common. + +Note: The *-connect-angle* argument has no effect unless the +["-surface" and "-connect"](#Detecting-membranes) arguments are also supplied. + +#### -cts, -ctn, -cvs, -cvn + +*(Warning: This is an advanced experimental feature. +These arguments can probably be ignored by most users.)* + +Alternatively, if you need more detailed control over the criteria +used to decide whether to add voxels to an existing cluster, +then you can specify 4 different thresholds using the +"**-cts threshold**", "**-ctn threshold**", +"**-cvs threshold**", and "**-cvn *threshold**" arguments. +This is an *alternative* to (specifying the "**-connect-angle**" argument). + + +Motivation and explanation: +Ideally, if tensor-voting worked, then the results of tensor voting should +agree with the actual direction of the features located at each voxel location. +If not, then you can exclude these voxels from consideration. +To do that, the "*-cts*" and "*-cvs*" threshold arguments specify +the minimum similarity between the hessian of the image brightness with the +tensor and vector features that resulted from tensor voting, respectively. +Additionally, the tensors from *neighboring* voxels which are part of the +same surface or curve should be poinging in the same direction. +So the "*-ctn*" and "*-cvn*" threshold arguments specify is the +minimum similarity between the orientations of the tensors and vectors +from neighboring voxels, respectively. +All of these threshold parameters should be in the range from -1 to 1. +A *-cvn -threshold* value of 0.707 ≈ cos(45°) and corresponds to a +45 degree difference between orientations of neighboring voxels. +In that case, neighboring voxels pointing in directions which +differ by more than 45° will be assigned to different clusters (membranes). +(This works reasonably well for phospholipid membranes.) +Note that the *-cts*, *-ctn*, *-cvs*, and *-cvn* arguments +have no effect unless the ["-surface" and "-connect"](#Detecting-membranes) +arguments are also in use. + +# Filters + +### -gauss and -ggauss +The **-gauss** and **-gauss-aniso** arguments must be followed by one or more numbers specifying the width of the Gaussian filter to apply to your image: +``` + -gauss σ +``` +or +``` + -gauss-aniso σ_x σ_y σ_z +``` +In either case, the original image is convolved with: +``` + h(x,y,z) = A*exp(-0.5 * ((x/σ_x)^2 + (y/σ_y)^2 + (z/σ_z)^2)) +``` +When *-gauss* is used, *σ_x = σ_y = σ_z = σ* + +If **-ggauss** or **-ggauss-aniso** is used, +then the image is instead convolved with +[a generalized gaussian function](https://en.wikipedia.org/wiki/Generalized_normal_distribution#Version_1): +``` + h(x,y,z) = A*exp(-r^m) + where r = √((x/σ_x)^2 + (y/σ_y)^2 + (z/σ_z)^2)) +``` +The width of the Gaussian (ie, the σ_x, σ_y, σ_z arguments) should be specified in units of physical distance, not in voxels. +(The A coefficient will be determined automatically by normalization.) +By default **m** *= 2*, +however this can be overridden using the optional "**-exponent m**" +*Note:* that these *σ*, *σ_x*, *σ_y*, *σ_z*, parameters are a factor +of *√2 times larger* than the corresponding *σ* ("sigma") parameter +traditionally used in the +[Gaussian function](https://en.wikipedia.org/wiki/Normal_distribution). +(For your convenience, you can plot these functions for various parameters using +the "draw_filter1D.py -ggauss A s m" python script which is located in the same directory. +In this case, the "s" parameter is specified in voxels, not physical distance.) +*Note:* The calculation is +[fast](https://en.wikipedia.org/wiki/Separable_filter) +if you use the default exponent of 2. +*Changing the exponent will slow down the filter considerably.* + +The filter is truncated far away from the central peak at a point which is chosen automatically according the σ, σ_x, σ_y, σ_z parameters selected by the user. However this can be customized using the +["-truncate-threshold"](#Filter-Size) +and +["-truncate"](#Filter-Size) +arguments if necessary (see below). + +Note: By default, the image blurring that results from applying a +Gaussian filter is adjusted near the image boundary (or *mask* boundary) +so that the average brightness there is not effected. +This can be disabled using the +[-normalize-filters no](#-normalize-filters-yes/no) argument. + ## Thresholding, Clipping and Rescaling (Intensity Map Filters): Thresholding filters provide a way to rescale (map) the @@ -1632,17 +1598,6 @@ a narrow range of interest. *Note:* All threshold operations are performed *after* normal filtering operations have been applied (such as -gauss, -dog, -dogg, -fluct, or -log filters). -### -invert - -This filter replaces bright voxels with dark ones, and visa versa, -while keeping the average voxel brightness the same. -(To calculate the new voxel intensity, - the *difference* between the original voxel intensity - and the average intensity is *subtracted* from the average intensity.) -Inverting an image can be useful to prepare files which you plan -to read with other programs (such as ChimeraX and IMOD). - - ### -rescale m b @@ -1680,15 +1635,32 @@ image voxel intensities in the final image will be shifted and rescaled so that the minimum and maximum voxel intensities will be *outA* and *outB*. -*(Note: As of 2018-6-30, the "Black" and "White" controls in *IMOD* - report values between 0 and 255, even if the actual voxel brightnesses - in the file are floating point numbers (between 0 and 1, for example). - Ignore those numbers. - You can find the brightness of a particular voxel by left-clicking - on that voxel (to move the yellow-cursor) - and then by selecting the "Edit"->"Point"->"Value" menu option within IMOD. - You can also use the "histogram_mrc.py" script to - see if your intensity values in the image lie in the range you expect.)* +*(Note: You can use the "histogram_mrc.py" script to +see if your intensity values in the image lie in the range you expect. +You can find the brightness of a particular voxel in IMOD (3dmod) +by left-clicking on that voxel (to move the yellow-cursor), and then by +selecting the "Edit"->"Point"->"Value" menu option (or pressing the "F" key). +However, for 8-bit MRC files, IMOD/3dmod will report the brightness values +as unsigned integers between 0-255, shifting the brightness values upwards +to avoid negative values. So if you have an 8-bit MRC or REC +file, you should use IMOD's "newstack" program to convert the image into 32-bit +floats using "newstack -input mask_file.rec -output new_file.rec -mode 2". +Then use 3dmod to view the 32-bit version of the image, the brightness values +should be reported correctly.) + + + +### -invert + +This filter replaces bright voxels with dark ones, and visa versa, +while keeping the average voxel brightness the same. +(To calculate the new voxel intensity, + the *difference* between the original voxel intensity + and the average intensity is *subtracted* from the average intensity.) +Inverting an image can be useful to prepare files which you plan +to read with other programs (such as ChimeraX and IMOD). + + ### -clip a b @@ -1916,9 +1888,181 @@ centered around *x0* with standard deviation σ. resulting voxel intensities from outA to outB (instead of from 0 to 1).* +### -dog +When the "**-dog**" or "**-dog-aniso**" filter is selected, a +[Difference-of-Gaussians](https://en.wikipedia.org/wiki/Difference_of_Gaussians) +filter will be used. +The **-dog** and **-dog-aniso** arguments must be followed by numbers indicating the width of the Gaussians you wish to convolve with your image. +(Equivalently, the numbers indicate the band of frequencies you wish to keep from your original image.) +``` + -dog a b +``` +or: +``` + -dog-aniso a_x a_y a_z b_x b_y b_z +``` +In either case, +The original image will be convolved with the following function: +``` + h(x,y,z) = A*exp(-0.5 * ((x/a_x)^2 + (y/a_y)^2 + (z/a_z)^2)) + - B*exp(-0.5 * ((x/b_x)^2 + (y/b_y)^2 + (z/b_z)^2)) +``` +When *-dog* is used, +*a_x = a_y = a_z = a* +and +*b_x = b_y = b_z = b* + + +### LoG filters +### -log, -log-r, -log-d, and -log-aniso + + +When the "**-log**", "**-log-r**", and "**-log-d**", or "**-log-aniso**" +arguments are selected, a +[Laplacian-of-a-Gaussian (LoG)](https://en.wikipedia.org/wiki/Blob_detection) +filter is applied on the source image. +(See implementation details below.) +The Laplacian-of-a-Gaussian filter applies a Gaussian blur to the image +(whose width is specified by the user), +and then calculates the Laplacian of the result. +Features in the image of various sizes can be emphasized +using this kind of filter, and the result can be saved as a new image using +["**-out** filename.rec"](#Input-and-Output-files). +The "**-log**", "**-log-aniso**", "**-log-r**", and "**-log-d**" +arguments described here are typically only chosen when the user already +knows the approximate size of the features they want to emphasize +in the image. +``` + -log σ +``` +When "**-log σ**" is used, a LoG filter will be applied +to the image using a Gaussian width (σ). +Alternatively, if the user wishes to specify the actual size +(the radius or diameter) of the objects they want to emphasize, +then they can (equivalently) use +``` + -log-r radius +``` +or +``` + -log-d diameter +``` +instead. +There is also an anisotropic version of this filter as well: +``` + -log-aniso σ_x σ_y σ_z** +``` +The new image which is generated will *tend* to have have bright and dark spots +wherever bright and dark objects of the corresponding size can be found. +This image can be searched for (local) minima and maxima +by running this program again on the new image using the +"**-find-minima**", "**-find-maxima**" and "-radial-separation" arguments. +These locations can be saved to a file and then a new anotated image can be +created using the +["**-draw-spheres**"](#-draw-spheres-filename) +argument. This provides a fast, sloppy, +and unselective way to search for features in an image. +However scale-free +[**blob detection**](https://en.wikipedia.org/wiki/Blob_detection) +is a more robust (and computationally expensive) way to detect objects +of a given size within an image. + +*(Implementation note: The LoG filter described here is approximated internally + using a DoG filter. You can control the accuracy of this approximation using + the + ["-dog-delta δ"](#Implementation-of-LoG-filters-and-blob-detectors) + argument.) + + + +### -fluct and -fluct-aniso +Usage: +``` + -fluct r +``` +The "**-fluct**" filter calculates the magnitude of +the voxel intensity fluctuations relative to their local surroundings. +This is useful for finding regions in an image with nearly uniform brightness +and regions whose brightness fluctuates wildly. +*r* is the effective search radius. + +To compute this, near every voxel, a Gaussian-weighted average +intensity is computed with the same effective radius, *r*. +A Gaussian-weighted squared-difference between all the nearby voxel intensities +and this local average intensity is also computed. +A new image is generated whose voxel intensities equal this +local average squared difference intensity of nearby voxels. + +The "**-fluct-aniso**" variant allows the user to control the width +of the Gaussian independently in the x,y,z directions: +``` + -fluct-aniso r_x r_y r_z +``` +#### Implementation detail: + The width of the Gaussian, σ, is chosen so that the effective + number of voxels under the 3D Gaussian, + (which can be interpreted as the 3D integral of the unnormalized 3D Gaussian, + (2\*π\*σ^2)^(3/2)) + equals the volume of a sphere of radius r + (4/3)\*π\*r^3. + This yields σ=(9π/2)^(-1/6)r. + + It might seem more straightforward to simply consider the fluctuations in + brightnesses of all of the voxels which lie within a fixed radius from + each target voxel. However using Gaussian weighted averages instead + accomplishes essentially the same goal and is much more + [computationally efficient](https://en.wikipedia.org/wiki/Separable_filter). + + If you prefer, instead of using a Gaussian, you can instead perform the + fluctuation calculation over a hard spherical region by using the + "-exponent n" argument. + (This parameter is 2 by default. Changing it will slow down the calculation.) + As n is increased, the region over which the fluctuations in brightness + are calculated becomes more and more like a uniform sphere with radius σ. + (So if you do this, be sure to multiply your radius argument, r, + by (9π/2)^(1/6)≈1.5549880806696572 to compensate.) + + + + + +## Resizing the image + +### -bin binsize + +Reduce the resolution of the image in each direction by a factor of *binsize*. +(*binsize* must be a positive integer.) +The new image will be smaller in each direction by a factof of *binsize*. + +Motivation: Reducing the resolution of an image can often dramatically +reduce the computation time necessary to detect objects in the image. +If you are detecting features in the image (such as membranes or blobs) +that are significantly larger (or thicker) than the voxel width, +then a reduction in resolution by a factor of 2 or so +should not effect your ability to detect it accurately, +(and could make tensor voting up to 64 times faster (as of 2021-6-21)). +Note that you can use this argument together with other arguments. +If you do that, the reduction of resolution occurs before all +other operations are performed (such as "-blob" or "-surface" detection). +This allows you to reduce the image size before +these other expensive calculations are performed. + +Note that the width of each voxel will be increased accordingly, +so that the positions of any features in the image that you detect +should not be significantly effected by the change in resolution. +(This includes the coordinates of the surface mesh.) +Similarly, when using "-bin", the user does not need to alter +any of the other parameters or arguments passed to the program +(such as the thickness parameters or sigma parameters +that characterize the size of the objects they want to detect). +The software will automatically adjust these parameters +to compensate for the change in image resolution. -## Annotation of detected objects in an image + + + +# Annotation of detected objects in an image ### (-draw-spheres filename) ### (-draw-hollow-spheres filename) @@ -1952,8 +2096,8 @@ and it superimposes them upon the original image. Each line in the file contains information describing a different sphere. The first three numbers on each line are assumed to store the x,y, and z coordinates of the center of that sphere. -*(Unless the "-w 1" argument is used, the coordinates and diameters -are assumed to be in physical units (eg. Angstroms), not voxels.)* +*(Unless the ["-w 1"](#Voxel-Width) argument is used, the coordinates and +diameters are assumed to be in physical units (eg. Angstroms), not voxels.)* If the file contains a 4th column, then it is assumed to store the diameter of each sphere. (Otherwise the "spheres" will be only 1 voxel wide by default.) @@ -1974,7 +2118,6 @@ or by pressing the "F" key. The brightness of voxel at that location will be printed to the IMOD control window. That brightness is the score of the corresponding blob.)* - #### Background voxels Again, by default, the background voxels will be the same as the voxels @@ -2034,7 +2177,6 @@ the *ratio* parameter between 0 and 1 until the image looks reasonable. (Typical values of the *ratio* parameter are between 0.1 and 0.3.) - #### Thin hollow spheres If "**-draw-hollow-spheres**" is used, then hollow spheres are drawn instead. @@ -2044,7 +2186,7 @@ The thickness of each hollow spherical shell can be controlled using the In the first example, the *thickness* of the sphere is the same for all spheres, regardless of sphere size. In that case, the *thickness* argument should be specified in physical distance units -(*or* in voxels, if you are using the "*-w 1*" argument). +(*or* in voxels, if you are using the ["-w 1"](#Voxel-Width) argument). On the other hand, if you use the "**-spheres-shell-ratio ratio**" argument, the thickness is expressed as a ratio, relative to the radius of the sphere. (Larger spheres will have thicker shells.) @@ -2147,8 +2289,7 @@ sphere. If the small sphere's volume is less than half of the larger sphere To prevent the small sphere from being discarded, you can use [**-max-volume-overlap 0.5**](#-max-volume-overlap-fraction) together with -[**-max-volume-overlap-small 1**](#-max-volume-overlap-small-fraction) -. +[**-max-volume-overlap-small 1**](#-max-volume-overlap-small-fraction). #### Automatic disposal of poor scoring blobs @@ -2370,7 +2511,6 @@ or the Gaussian width (σ) of the objects that you wish to detect within the image (instead of the object's diameter). - #### Visualizing blobs using the "*-draw-spheres*" argument Again, after @@ -2387,9 +2527,7 @@ arguments. - - -## Masking +# Masking The optional "-mask", "-mask-select", and "-mask-out" arguments allow you to ignore certain voxels from the source image (tomogram). @@ -2569,89 +2707,18 @@ The resulting mask region can be visualized using the ### -fill brightness -Replace all of the voxel brightnesses with *brightness*. -If masks are used, then voxels outside the mask will be omitted. -When using the *-mask-rect* and *-mask-sphere* arguments, -using "**-fill 1**" will generate a 3-D image -(which you can view in IMOD, for example) -This provides a way to see where the mask region is located. - - - -### -erode-gauss thickness -### -dilate-gauss thickness - - -*(Note: This is a crude implementation of -[image dilation](https://en.wikipedia.org/wiki/Dilation_(morphology)) and -[image erosion](https://en.wikipedia.org/wiki/Erosion_(morphology)) -For a more flexible and sophisticated implementation in python, try using -[dilate](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.dilation) -[erode](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.erosion) -[binary_dilate](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.binary_dilation) -[binary_erode](https://scikit-image.org/docs/dev/api/skimage.morphology.html?highlight=skeletonize#skimage.morphology.binary_erosion) -from -[scikit-image](https://scikit-image.org), -along with the -[mrcfile](https://mrcfile.readthedocs.io/en/latest/readme.html#basic-usage) -module.)* - -The **-erode-gauss** and **-dilate-gauss** arguments are useful -for modifying the size of an existing binary image (eg. mask), -making the region appear bigger or smaller. -(This is useful to modify an image that you eventually intend -to use with the "-mask" argument in a future invocation of filter_mrc. -But these operations work on any binary image.) - -Recall that a *mask image* is a 3-D image used with the -["**-mask**"](#-mask-MRC_FILE) argument. -The brightness value of each voxel in a mask image is typically either 0 or 1 -*(depending on whether the voxel is supposed to be -ignored or considered during subsequent calculations)*. -I will refer to an image (of 1s and 0s) as a "binary image". -The **-dilate-gauss** and **-erode-gauss** arguments will enlarge -or shrink the size of the region stored in a binary image by -a distance of roughly *thickness* (which has units of physical distance -unless the "*-w 1*" argument is used). - -Note that, due to blurring, features in the image which are smaller or -narrower than *thickness* will be lost after this operation is completed. -Consquently, for large distances, it might work better to -repeat this operation multiple time with the same small *thickness* argument, -instead of using a single large *thickness* argument. -You will have to experiment to see what works best. +Replace all of the voxel brightnesses in the image with *brightness*. +If a mask is in use, then voxels outside the mask will be omitted. +*(You can specify the brightness of these voxels outside the mask using the +[-mask-out](#-mask-out-intensity_value) argument, which is 0 by default.)* -#### Details -This is implemented by blurring the image, and then applying a threshold -so that all voxels whose brightness are below a threshold are discarded. -Consequently, you cannot combine these arguments with either the *-gauss* -or *-thresh* arguments, because: +If you use "**-fill 1**", this will generate a 3-D image showing where the +mask region is located. (This can be useful if you are using the +[-mask-sphere](#-mask-sphere-x0-y0-z0-r) and +[-mask-rect](#-mask-rect-xmin-xmax-ymin-ymax-zmin-zmax) arguments. +You can view this image in IMOD, for example.) -- "-dilate-gauss *thickness*" is equivalent to -"-gauss *thickness* -thresh 0.157299", *(where 0.157299 ≈ 1-erf(1))* -- "-erode-gauss *thickness*" is equivalent to -"-gauss *thickness* -thresh 0.842701", *(where 0.842701 ≈ erf(1))* -The motivation for this strategy is outlined below: -If you start with a binary image (consisting of 1s and 0s), -and blur it, and discard all voxels whose brightness is less than 0.5, -(setting their brightness to 0, and the other voxel brightnesses to 1), -then you will get back a new binary image which strongly resembles -the original binary image -(with some high frequencies removed due to the blurring). -It will have roughly the same number of bright and dark voxels -as the original image. -But if, after blurring, you discard all voxels with threshold brightness -less than, say, 0.75 (instead of 0.5), then this is a more stringent -requirement. As a result, *fewer* voxels will remain, compared to if -you had chosen 0.5. -This effectively shrinks the size of the selected region (consisting of 1s) -from the binary image. Alternatively, if you only discard voxels whose -brightness less than, say, 0.25, then this is a more forgiving criteria. -This will select *more* voxels, effectively expanding the size of the -selected region from the binary image. -The thresholds above were chosen so that the selected region -expands or shrinks by roughly a distance of *thickness*. ## Miscellaneous @@ -2691,6 +2758,29 @@ the running time is proportional to the sum of these numbers, Wx+Wy+Wz. Keep this in mind when specifying filter window widths.)* +### -normalize-filters yes/no +When a Gaussian blur filter is applied to an image, each resulting voxel has +a brightness which is a (weighted) average of the brightnesses of nearby voxels. +By default, voxels which are *outside* the boundaries of the image +(or outside the *mask*), will not be included in this average. +Averaging is only done over the voxels which *are* within +the boundaries. *(In other words, the weights given to the remaining voxels +are increased when computing the average nearby brightness.)* +Unfortunately, by averaging fewer voxels near the image boundaries, +the resulting image is likely to be noiser at these locations. +If this is not what you want, you can turn this feature off using +``` +-normalize-filters no +``` +If you do that, the image may appear artificially dark or bright +near the image or mask boundaries. + +This setting effects the "-gauss" and "-ggauss" arguments. +It also effects blob detection, surface detection, curve detection, +and any other filters which apply a Gaussian blur to the image +at some point during the calculation. + + ### Distance Units: Angstroms or Nanometers ``` -a2nm diff --git a/lib/visfd/clustering.hpp b/lib/visfd/clustering.hpp index afea9a6..048f1ad 100644 --- a/lib/visfd/clustering.hpp +++ b/lib/visfd/clustering.hpp @@ -66,6 +66,10 @@ typedef enum eClusterSortCriteria { /// voxels which lie near a maxima (or a minima) of saliency, /// and groups voxels of similar saliency together (as well as voxels /// whose vector and/or tensor directions are compatible, if applicable) +/// Unlike the watershed algorithm, this function can cluster voxels +/// according to their vector and tensor attributes, clustering +/// adjacent voxels only if their directions and tensors are +/// sufficiently similar. /// /// @return The function does not have a return value. /// After the function is finished, the aaaiDest[][][] array will diff --git a/lib/visfd/morphology.hpp b/lib/visfd/morphology.hpp index 354b51b..f3fc376 100644 --- a/lib/visfd/morphology.hpp +++ b/lib/visfd/morphology.hpp @@ -221,89 +221,109 @@ FindExtrema(int const image_size[3], //!< size of input image array -/// @brief Computes the grayscale dilation of an image as defined here: +/// @brief Computes the grayscale dilation of an image with an arbitrary +/// structure factor, as defined here: /// https://en.wikipedia.org/wiki/Dilation_(morphology)#Grayscale_dilation /// The structure factor is implemented as a vector of tuples. /// Each tuple contains the voxel coordinates and brightness value ("b") /// for the structuer factor at that location. -template +/// (A version of this function using a spherical structure factors +/// is defined elsewhere.) +template void -Dilate(vector > structure_factor, // a list of (ix,iy,iz,b) values, one for each voxel - Integer const image_size[3], //!< size of the image in x,y,z directions +Dilate(vector > structure_factor, // a list of (ix,iy,iz,b) values, one for each voxel + const int image_size[3], //!< size of the image in x,y,z directions Scalar const *const *const *aaafSource, //!< image array aaafI[iz][iy][ix] Scalar ***aaafDest, //!< image array aaafI[iz][iy][ix] Scalar const *const *const *aaafMask = nullptr, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 ostream *pReportProgress = nullptr) { - for (Integer iz=0; iz < image_size[2]; iz++) { + for (int iz=0; iz < image_size[2]; iz++) { if (pReportProgress) *pReportProgress << " z = " << iz+1 << " of " << image_size[2] << endl; #pragma omp parallel for collapse(2) - for (Integer iy=0; iy < image_size[1]; iy++) { - for (Integer ix=0; ix < image_size[0]; ix++) { + for (int iy=0; iy < image_size[1]; iy++) { + for (int ix=0; ix < image_size[0]; ix++) { Scalar max_f_plus_b = -std::numeric_limits::infinity(); for (auto pVoxel = structure_factor.begin(); pVoxel != structure_factor.end(); pVoxel++) { - Integer jx = (*pVoxel)[0]; - Integer jy = (*pVoxel)[1]; - Integer jz = (*pVoxel)[2]; - if ((aaafMask) && (aaafMask[iz+jz][iy+jy][ix+jx] == 0.0)) + int jx = std::get<0>(*pVoxel); + int jy = std::get<1>(*pVoxel); + int jz = std::get<2>(*pVoxel); + int Ix = ix + jx; + int Iy = iy + jy; + int Iz = iz + jz; + if ((Ix < 0) || (Ix >= image_size[0]) || + (Iy < 0) || (Iy >= image_size[1]) || + (Iz < 0) || (Iz >= image_size[2])) continue; - Scalar b = (*pVoxel)[3]; - Scalar f = aaafSource[iz+jz][iy+jy][ix+jx]; + if ((aaafMask) && (aaafMask[Iz][Iy][Ix] == 0.0)) + continue; + Scalar b = std::get<3>(*pVoxel); + Scalar f = aaafSource[Iz][Iy][Ix]; max_f_plus_b = std::max(max_f_plus_b, f + b); } aaafDest[iz][iy][ix] = max_f_plus_b; } } } -} // Dilate(vector > structure_factor,...) +} // Dilate(vector > structure_factor,...) -/// @brief Computes the grayscale erosion of an image as defined here: +/// @brief Computes the grayscale erosion of an image with an arbitrary +/// structure factor, as defined here: /// https://en.wikipedia.org/wiki/Erosion_(morphology)#Grayscale_erosion /// The structure factor is implemented as a vector of tuples. /// Each tuple contains the voxel coordinates and brightness value ("b") /// for the structuer factor at that location. -template +/// (A version of this function using a spherical structure factors +/// is defined elsewhere.) +template void -Erode(vector > structure_factor, // a list of (ix,iy,iz,b) values, one for each voxel - Integer const image_size[3], //!< size of the image in x,y,z directions +Erode(vector > structure_factor, // a list of (ix,iy,iz,b) values, one for each voxel + const int image_size[3], //!< size of the image in x,y,z directions Scalar const *const *const *aaafSource, //!< image array aaafI[iz][iy][ix] Scalar ***aaafDest, //!< image array aaafI[iz][iy][ix] Scalar const *const *const *aaafMask = nullptr, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 ostream *pReportProgress = nullptr) { - for (Integer iz=0; iz < image_size[2]; iz++) { + for (int iz=0; iz < image_size[2]; iz++) { if (pReportProgress) *pReportProgress << " z = " << iz+1 << " of " << image_size[2] << endl; #pragma omp parallel for collapse(2) - for (Integer iy=0; iy < image_size[1]; iy++) { - for (Integer ix=0; ix < image_size[0]; ix++) { - Scalar min_f_minus_b = -std::numeric_limits::infinity(); + for (int iy=0; iy < image_size[1]; iy++) { + for (int ix=0; ix < image_size[0]; ix++) { + Scalar min_f_minus_b = std::numeric_limits::infinity(); for (auto pVoxel = structure_factor.begin(); pVoxel != structure_factor.end(); pVoxel++) { - Integer jx = (*pVoxel)[0]; - Integer jy = (*pVoxel)[1]; - Integer jz = (*pVoxel)[2]; - if ((aaafMask) && (aaafMask[iz+jz][iy+jy][ix+jx] == 0.0)) + int jx = std::get<0>(*pVoxel); + int jy = std::get<1>(*pVoxel); + int jz = std::get<2>(*pVoxel); + int Ix = ix + jx; + int Iy = iy + jy; + int Iz = iz + jz; + if ((Ix < 0) || (Ix >= image_size[0]) || + (Iy < 0) || (Iy >= image_size[1]) || + (Iz < 0) || (Iz >= image_size[2])) + continue; + if ((aaafMask) && (aaafMask[Iz][Iy][Ix] == 0.0)) continue; - Scalar b = (*pVoxel)[3]; - Scalar f = aaafSource[iz+jz][iy+jy][ix+jx]; + Scalar b = std::get<3>(*pVoxel); + Scalar f = aaafSource[Iz][Iy][Ix]; min_f_minus_b = std::min(min_f_minus_b, f - b); } aaafDest[iz][iy][ix] = min_f_minus_b; } } } -} // Erode(vector > structure_factor,...) +} // Erode(vector > structure_factor,...) @@ -313,21 +333,24 @@ Erode(vector > structure_factor, // a list /// https://en.wikipedia.org/wiki/Dilation_(morphology)#Grayscale_dilation /// The radius paramater can be a floating point number, but its /// units are in voxels. -template +template void DilateSphere(Scalar radius, //!< radius of the sphere int const image_size[3], //!< size of the image in x,y,z directions Scalar const *const *const *aaafSource, //!< image array aaafI[iz][iy][ix] Scalar ***aaafDest, //!< image array aaafI[iz][iy][ix] Scalar const *const *const *aaafMask = nullptr, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 - bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels + bool smooth_boundary = false, //!< attempt to produce a rounder smoother structure factor that varies between 0 and -1 at its boundary voxels ostream *pReportProgress = nullptr //!< report progress to the user ) { - vector > structure_factor; // a list of (ix,iy,iz,b) values, one for each voxel - for (Integer iz=0; iz < image_size[2]; iz++) { - for (Integer iy=0; iy < image_size[1]; iy++) { - for (Integer ix=0; ix < image_size[0]; ix++) { + // structure_factor = a vector of (ix,iy,iz,b) values, one for each voxel + vector > structure_factor; + + int Ri = ceil(radius); + for (int iz = -Ri; iz <= Ri; iz++) { + for (int iy = -Ri; iy <= Ri; iy++) { + for (int ix = -Ri; ix <= Ri; ix++) { bool add_this_voxel = false; Scalar b = 0.0; if (! smooth_boundary) { @@ -352,9 +375,9 @@ DilateSphere(Scalar radius, //!< radius of the sphere // Here, I am doing it a sloppy way, but it should be good enough.) Scalar r_max = -std::numeric_limits::infinity(); Scalar r_min = std::numeric_limits::infinity(); - for (Integer jz=0; jz <= 1; jz++) { - for (Integer jy=0; jy <= 1; jy++) { - for (Integer jx=0; jx <= 1 ; jx++) { + for (int jz=0; jz <= 1; jz++) { + for (int jy=0; jy <= 1; jy++) { + for (int jx=0; jx <= 1 ; jx++) { Scalar r = sqrt(SQR(ix+jx-0.5)+SQR(iy+jy-0.5)+SQR(iz+jz-0.5)); if (r < r_min) r_min = r; @@ -378,10 +401,10 @@ DilateSphere(Scalar radius, //!< radius of the sphere } } // if (smooth_boundary) if (add_this_voxel) - structure_factor.push_back(tuple(ix,iy,iz,b)); - } // for (Integer ix=0; ix < image_size[0]; ix++) - } // for (Integer iy=0; iy < image_size[1]; iy++) - } // for (Integer iz=0; iz < image_size[2]; iz++) + structure_factor.push_back(tuple(ix,iy,iz,b)); + } // for (int ix=0; ix < image_size[0]; ix++) + } // for (int iy=0; iy < image_size[1]; iy++) + } // for (int iz=0; iz < image_size[2]; iz++) if (pReportProgress) *pReportProgress << "DilateSphere() progress:" << endl; @@ -402,21 +425,24 @@ DilateSphere(Scalar radius, //!< radius of the sphere /// https://en.wikipedia.org/wiki/Dilation_(morphology)#Grayscale_dilation /// The radius paramater can be a floating point number, but its /// units are in voxels. -template +template void ErodeSphere(Scalar radius, //!< radius of the sphere int const image_size[3], //!< size of the image in x,y,z directions Scalar const *const *const *aaafSource, //!< image array aaafI[iz][iy][ix] Scalar ***aaafDest, //!< image array aaafI[iz][iy][ix] Scalar const *const *const *aaafMask = nullptr, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 - bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels + bool smooth_boundary = false, //!< attempt to produce a rounder smoother structure factor that varies between 0 and -1 at its boundary voxels ostream *pReportProgress = nullptr //!< report progress to the user ) { - vector > structure_factor; // a list of (ix,iy,iz,b) values, one for each voxel - for (Integer iz=0; iz < image_size[2]; iz++) { - for (Integer iy=0; iy < image_size[1]; iy++) { - for (Integer ix=0; ix < image_size[0]; ix++) { + // structure_factor = a vector of (ix,iy,iz,b) values, one for each voxel + vector > structure_factor; + + int Ri = ceil(radius); + for (int iz = -Ri; iz <= Ri; iz++) { + for (int iy = -Ri; iy <= Ri; iy++) { + for (int ix = -Ri; ix <= Ri; ix++) { bool add_this_voxel = false; Scalar b = 0.0; if (! smooth_boundary) { @@ -441,9 +467,9 @@ ErodeSphere(Scalar radius, //!< radius of the sphere // Here, I am doing it a sloppy way, but it should be good enough.) Scalar r_max = -std::numeric_limits::infinity(); Scalar r_min = std::numeric_limits::infinity(); - for (Integer jz=0; jz <= 1; jz++) { - for (Integer jy=0; jy <= 1; jy++) { - for (Integer jx=0; jx <= 1 ; jx++) { + for (int jz=0; jz <= 1; jz++) { + for (int jy=0; jy <= 1; jy++) { + for (int jx=0; jx <= 1 ; jx++) { Scalar r = sqrt(SQR(ix+jx-0.5)+SQR(iy+jy-0.5)+SQR(iz+jz-0.5)); if (r < r_min) r_min = r; @@ -466,10 +492,10 @@ ErodeSphere(Scalar radius, //!< radius of the sphere } } // if (smooth_boundary) if (add_this_voxel) - structure_factor.push_back(tuple(ix,iy,iz,b)); - } // for (Integer ix=0; ix < image_size[0]; ix++) - } // for (Integer iy=0; iy < image_size[1]; iy++) - } // for (Integer iz=0; iz < image_size[2]; iz++) + structure_factor.push_back(tuple(ix,iy,iz,b)); + } // for (int ix=0; ix < image_size[0]; ix++) + } // for (int iy=0; iy < image_size[1]; iy++) + } // for (int iz=0; iz < image_size[2]; iz++) if (pReportProgress) *pReportProgress << "ErodeSphere() progress:" << endl; @@ -485,6 +511,100 @@ ErodeSphere(Scalar radius, //!< radius of the sphere + +/// @brief Computes the grayscale opening of an image with a flat sphere. +/// https://en.wikipedia.org/wiki/Opening_(morphology) +/// The radius paramater can be a floating point number, but its +/// units are in voxels. +template +void +OpenSphere(Scalar radius, //!< radius of the sphere + const int image_size[3], //!< size of the image in x,y,z directions + Scalar const *const *const *aaafSource, //!< image array aaafI[iz][iy][ix] + Scalar ***aaafDest, //!< image array aaafI[iz][iy][ix] + Scalar const *const *const *aaafMask = nullptr, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels + ostream *pReportProgress = nullptr //!< report progress to the user + ) +{ + float *afTmp; //temporary space to store image after erosion + float ***aaafTmp; //temporary space to store image after erosion + + Alloc3D(image_size, + &afTmp, + &aaafTmp); + + ErodeSphere(radius, + image_size, + aaafSource, + aaafTmp, //<--save the resulting image in aaafTmp + aaafMask, + smooth_boundary, + pReportProgress); + + DilateSphere(radius, + image_size, + aaafTmp, //<--use aaafTmp as the source image + aaafDest, //<--save the resulting image in aaafDest + aaafMask, + smooth_boundary, + pReportProgress); + + Dealloc3D(image_size, + &afTmp, + &aaafTmp); + +} // OpenSphere() + + + + +/// @brief Computes the grayscale closing of an image with a flat sphere. +/// https://en.wikipedia.org/wiki/Closing_(morphology) +/// The radius paramater can be a floating point number, but its +/// units are in voxels. +template +void +CloseSphere(Scalar radius, //!< radius of the sphere + const int image_size[3], //!< size of the image in x,y,z directions + Scalar const *const *const *aaafSource, //!< image array aaafI[iz][iy][ix] + Scalar ***aaafDest, //!< image array aaafI[iz][iy][ix] + Scalar const *const *const *aaafMask = nullptr, //!< optional: ignore voxel ix,iy,iz if aaafMask!=nullptr and aaafMask[iz][iy][ix]==0 + bool smooth_boundary = false, //!< attempt to produce a round and smooth structure factor that varies between 0 and -1 at its boundary voxels + ostream *pReportProgress = nullptr //!< report progress to the user + ) +{ + float *afTmp; //temporary space to store image after dilation + float ***aaafTmp; //temporary space to store image after dilation + + Alloc3D(image_size, + &afTmp, + &aaafTmp); + + DilateSphere(radius, + image_size, + aaafSource, + aaafTmp, //<--save the resulting image in aaafTmp + aaafMask, + smooth_boundary, + pReportProgress); + + ErodeSphere(radius, + image_size, + aaafTmp, //<--use aaafTmp as the source image + aaafDest, //<--save the resulting image in aaafDest + aaafMask, + smooth_boundary, + pReportProgress); + + Dealloc3D(image_size, + &afTmp, + &aaafTmp); + +} // CloseSphere() + + + } //namespace visfd