LOOS  v2.3.2
blob_stats.cpp
1 /*
2  blob_stats.cpp
3 
4  Gather statistics on blobs...
5 */
6 
7 /*
8  This file is part of LOOS.
9 
10  LOOS (Lightweight Object-Oriented Structure library)
11  Copyright (c) 2008, Tod D. Romo, Alan Grossfield
12  Department of Biochemistry and Biophysics
13  School of Medicine & Dentistry, University of Rochester
14 
15  This package (LOOS) is free software: you can redistribute it and/or modify
16  it under the terms of the GNU General Public License as published by
17  the Free Software Foundation under version 3 of the License.
18 
19  This package is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  GNU General Public License for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with this program. If not, see <http://www.gnu.org/licenses/>.
26 */
27 
28 
29 
30 #include <loos.hpp>
31 #include <boost/format.hpp>
32 #include <boost/tuple/tuple.hpp>
33 #include <algorithm>
34 #include <sstream>
35 #include <limits>
36 
37 #include <DensityGrid.hpp>
38 
39 using namespace std;
40 using namespace loos;
41 using namespace loos::DensityTools;
42 
43 
44 
45 
46 int countBlobs(const DensityGrid<int>& grid) {
47  DensityGridpoint dims = grid.gridDims();
48  long k = dims[0] * dims[1] * dims[2];
49  int n = 0;
50 
51  for (long i=0; i<k; i++)
52  if (grid(i) > n)
53  n = grid(i);
54 
55  return(n);
56 }
57 
58 
59 vector<int> sizeBlobs(const DensityGrid<int>& grid) {
60 
61  int n = countBlobs(grid);
62  vector<int> sizes(n+1, 0);
63 
64  DensityGridpoint dims = grid.gridDims();
65  long k = dims[0] * dims[1] * dims[2];
66 
67  for (long i=0; i<k; i++)
68  sizes[grid(i)] += 1;
69 
70  return(sizes);
71 }
72 
73 
74 vector<GCoord> blobCentroids(const int n, const DensityGrid<int>& grid) {
75  vector<GCoord> centers(n+1, GCoord(0,0,0));
76  vector<int> counts(n+1, 0);
77 
78  DensityGridpoint dims = grid.gridDims();
79  for (int k=0; k<dims[2]; k++)
80  for (int j=0; j<dims[1]; j++)
81  for (int i=0; i<dims[0]; i++) {
82  DensityGridpoint point(i,j,k);
83  GCoord u = grid.gridToWorld(point);
84  int id = grid(point);
85 
86  centers[id] += u;
87  counts[id] += 1;
88  }
89 
90  for (int i=0; i<= n; i++)
91  centers[i] /= counts[i];
92 
93  return(centers);
94 }
95 
96 
97 int main(int argc, char *argv[]) {
98  DensityGrid<int> grid;
99 
100  if (argc != 1) {
101  cerr <<
102  "Usage- blob_stats <foo.grid\n"
103  "\n"
104  "Print out basic information about blobs in a grid (requires an integer grid)\n"
105  "See also blobid, and pick_blob\n";
106  exit(-1);
107  }
108 
109  cin >> grid;
110  cout << "Read in grid with dimensions " << grid.gridDims() << endl;
111  cout << "Grid extents (real-space) is " << grid.minCoord() << " x " << grid.maxCoord() << endl;
112  GCoord range = grid.maxCoord() - grid.minCoord();
113  cout << "Grid range is " << range << endl;
114 
115  vector<int> sizes = sizeBlobs(grid);
116  int n = sizes.size();
117  vector<GCoord> centers = blobCentroids(n, grid);
118 
119  GCoord delta = grid.gridDelta();
120  double voxel_volume = 1.0 / delta[0];
121  for (int i=1; i<3; i++)
122  voxel_volume *= (1.0 / delta[i]);
123 
124  cout << boost::format("Voxel volume = %8.6g\n") % voxel_volume;
125  cout << boost::format("%6s %12s %12s\t%s\n")
126  % "Id"
127  % "Voxels"
128  % "Size (in A^3)"
129  % "Centroid (in A)";
130  cout << boost::format("%-6s %-12s %-12s\t%s\n")
131  % "------"
132  % "------------"
133  % "------------"
134  % "------------------------------";
135 
136 
137 
138  for (uint i = 0; i < sizes.size(); ++i) {
139  cout << boost::format("%6d %12d %12.6g\t") % i % sizes[i] % (sizes[i] * voxel_volume);
140  cout << centers[i] << endl;
141  }
142 
143 }
STL namespace.
Namespace for Density package.
Definition: DensityGrid.hpp:48
loos::GCoord gridToWorld(const DensityGridpoint &v) const
Converts grid coords to real-space (world) coords.
Namespace for most things not already encapsulated within a class.