LOOS  v2.3.2
blobid.cpp
1 /*
2  blobid.cpp
3 
4  Flood-fills a grid to identify blobs...
5 */
6 
7 
8 
9 
10 /*
11  This file is part of LOOS.
12 
13  LOOS (Lightweight Object-Oriented Structure library)
14  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
15  Department of Biochemistry and Biophysics
16  School of Medicine & Dentistry, University of Rochester
17 
18  This package (LOOS) is free software: you can redistribute it and/or modify
19  it under the terms of the GNU General Public License as published by
20  the Free Software Foundation under version 3 of the License.
21 
22  This package is distributed in the hope that it will be useful,
23  but WITHOUT ANY WARRANTY; without even the implied warranty of
24  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  GNU General Public License for more details.
26 
27  You should have received a copy of the GNU General Public License
28  along with this program. If not, see <http://www.gnu.org/licenses/>.
29 */
30 
31 
32 
33 #include <loos.hpp>
34 #include <limits>
35 
36 #include <DensityTools.hpp>
37 #include <DensityGrid.hpp>
38 #include <GridUtils.hpp>
39 
40 namespace opts = loos::OptionsFramework;
41 namespace po = loos::OptionsFramework::po;
42 
43 using namespace std;
44 using namespace loos;
45 using namespace loos::DensityTools;
46 
47 double lower, upper;
48 
49 // @cond TOOLS_INTERNAL
50 
51 
52 string fullHelpMessage(void) {
53  string msg =
54  "\n"
55  "SYNOPSIS\n"
56  "\n"
57  "\tIdentify blobs in a density grid.\n"
58  "\n"
59  "DESCRIPTION\n"
60  "\n"
61  "\tblobid identifies blobs by density values either in a range or above a threshold.\n"
62  "An edm grid (see for example water-hist) is expected for input.\n"
63  "Blobid then uses a flood-fill to determine how many separate blobs\n"
64  "meet the threshold/range criteria. A new grid is then written out\n"
65  "which identifies the separate blobs.\n"
66  "\nEXAMPLES\n"
67  "\tblobid --threshold 1 <foo.grid >foo_id.grid\n"
68  "Here we include all blobs above the threshold 1. foo_grid is a density\n"
69  "grid that has been created previously. For example a smoothed water \n"
70  "histogram grid may be used: \n"
71  "\twater-hist --radius=15 --bulk=25 --scale=1 b2ar.pdb b2ar.dcd |\\\n"
72  "\t grid2gauss 4 2 > foo_grid\n"
73  "The resulting blobs are then written to the grid \"foo_id\"\n"
74  "\n\n";
75 
76  return(msg);
77 }
78 
79 class ToolOptions : public opts::OptionsPackage {
80 
81  void addGeneric(po::options_description& o) {
82  o.add_options()
83  ("lower", po::value<double>(), "Sets the lower threshold for segmenting the grid")
84  ("upper", po::value<double>(), "Sets the upper threshold for segmenting the grid")
85  ("threshold", po::value<double>(), "Sets the threshold for segmenting the grid.");
86  }
87 
88  bool postConditions(po::variables_map& vm) {
89  if (vm.count("threshold")) {
90  lower = vm["threshold"].as<double>();
91  upper = numeric_limits<double>::max();
92  } else {
93  if (!(vm.count("lower-threshold") && vm.count("upper-threshold"))) {
94  cerr << "Error- you must specify either a treshold or a threshold range.\n";
95  return(false);
96  }
97 
98  lower = vm["lower-threshold"].as<double>();
99  upper = vm["upper-threshold"].as<double>();
100  }
101  return(true);
102  }
103 
104  string print() const {
105  ostringstream oss;
106 
107  oss << boost::format("lower=%f, upper=%f") % lower % upper;
108  return(oss.str());
109  }
110 
111 };
112 // @endcond
113 
114 
115 
116 
117 int fill(const DensityGridpoint seed, const int val, DensityGrid<double>& data_grid, DensityGrid<int>& blob_grid, const double low, const double high) {
118  vector<DensityGridpoint> blob = floodFill(seed, data_grid, val, blob_grid, ThresholdRange<double>(low, high));
119  return(blob.size());
120 }
121 
122 
123 boost::tuple<int, int, int, double> findBlobs(DensityGrid<double>& data_grid, DensityGrid<int>& blob_grid, const double low, const double high) {
124  DensityGridpoint dims = data_grid.gridDims();
125 
126  int id = 1;
127  int min = numeric_limits<int>::max();
128  int max = numeric_limits<int>::min();
129  double avg = 0.0;
130 
131  for (int k=0; k<dims[2]; k++) {
132  for (int j=0; j<dims[1]; j++)
133  for (int i=0; i<dims[0]; i++) {
134  DensityGridpoint point(i,j,k);
135  if (blob_grid(point) == 0 && (data_grid(point) >= low && data_grid(point) <= high)) {
136  int n = fill(point, id++, data_grid, blob_grid, low, high);
137  if (n < min)
138  min = n;
139  if (n > max)
140  max = n;
141  avg += n;
142  }
143  }
144  }
145 
146  avg /= (id-1);
147  boost::tuple<int, int, int, double> res(id-1, min, max, avg);
148  return(res);
149 }
150 
151 
152 int main(int argc, char *argv[]) {
153  string header = invocationHeader(argc, argv);
154 
155  opts::BasicOptions* bopts = new opts::BasicOptions(fullHelpMessage());
156  ToolOptions* topts = new ToolOptions;
157 
158  opts::AggregateOptions options;
159  options.add(bopts).add(topts);
160  if (!options.parse(argc, argv))
161  exit(-1);
162 
163 
164 
165  DensityGrid<double> data;
166  cin >> data;
167 
168  cerr << "Read in grid with size " << data.gridDims() << endl;
169 
170  DensityGrid<int> blobs(data.minCoord(), data.maxCoord(), data.gridDims());
171  boost::tuple<int, int, int, double> stats = findBlobs(data, blobs, lower, upper);
172  cerr << boost::format("Found %d blobs in range %6.4g to %6.4g\n") % boost::get<0>(stats) % lower % upper;
173  cerr << boost::format("Min blob size = %d, max blob size = %d, avg blob size = %6.4f\n")
174  % boost::get<1>(stats)
175  % boost::get<2>(stats)
176  % boost::get<3>(stats);
177 
178 
179  cout << blobs;
180 }
Namespace for encapsulating options processing.
STL namespace.
bool parse(int argc, char *argv[])
Parses a command line, returning true if parsing was ok.
Options common to all tools (including –fullhelp)
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Namespace for Density package.
Definition: DensityGrid.hpp:48
Functor that is true if value is between low and hi, inclusive.
Definition: GridUtils.hpp:54
Namespace for most things not already encapsulated within a class.
std::vector< DensityGridpoint > floodFill(const DensityGridpoint seed, const DensityGrid< T > &data_grid, const int id, DensityGrid< int > &blob_grid, const Functor &op)
Flood-fill a grid.
Definition: GridUtils.hpp:81
Combines a set of OptionsPackages.
AggregateOptions & add(OptionsPackage *pack)
Add a pointer to an OptionsPackage that will be used for options.