LOOS  v2.3.2
gridgauss.cpp
1 /*
2  Gridgauss
3 
4  Apply a gaussian kernel to a grid
5 */
6 
7 
8 /*
9  This file is part of LOOS.
10 
11  LOOS (Lightweight Object-Oriented Structure library)
12  Copyright (c) 2009, Tod D. Romo, Alan Grossfield
13  Department of Biochemistry and Biophysics
14  School of Medicine & Dentistry, University of Rochester
15 
16  This package (LOOS) is free software: you can redistribute it and/or modify
17  it under the terms of the GNU General Public License as published by
18  the Free Software Foundation under version 3 of the License.
19 
20  This package is distributed in the hope that it will be useful,
21  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23  GNU General Public License for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with this program. If not, see <http://www.gnu.org/licenses/>.
27 */
28 
29 
30 #include <loos.hpp>
31 #include <boost/format.hpp>
32 #include <DensityGrid.hpp>
33 #include <GridUtils.hpp>
34 
35 using namespace std;
36 using namespace loos;
37 using namespace loos::DensityTools;
38 
39 
40 int main(int argc, char *argv[]) {
41 
42  if (argc != 5) {
43  cerr <<
44  "DESCRIPTION\n\tApply a gaussian kernel convolution with a grid\n"
45  "\nUSAGE\n\tgridgauss width size scaling sigma <grid >output\n"
46  "Width controls the size (in grid units) of the kernel. Size\n"
47  "determines how the gaussian is mapped onto the kernel, i.e.\n"
48  "-size <= x < size. The gaussian is f(x) = exp(-0.5*(x/sigma)^2)\n"
49  "and is normalized so the sum of f(x) is one, then multiplied by\n"
50  "the scaling factor.\n"
51  "\nEXAMPLES\n\tgridgauss 10 3 1 1 <foo.grid >foo_smoothed.grid\n"
52  "This convolves the grid with a 10x10 kernel with sigma=1, and is a good\n"
53  "starting point for smoothing out water density grid.\n";
54  exit(0);
55  }
56 
57  string hdr = invocationHeader(argc, argv);
58 
59  int k = 1;
60  uint width = strtol(argv[k++], 0, 10);
61  double scaling = strtod(argv[k++], 0);
62  double normalization = strtod(argv[k++], 0);
63  double sigma = strtod(argv[k++], 0);
64 
65 
66  vector<double> kernel;
67  double scaling2 = 2.0 * scaling;
68 
69  double sum = 0.0;
70  for (uint i=0; i<width; ++i) {
71  double x = scaling2 * i / width - scaling;
72  double f = exp(-0.5*(x/sigma)*(x/sigma));
73  sum += f;
74  kernel.push_back(f);
75  }
76 
77  double a = normalization / sum;
78  for (uint i=0; i<kernel.size(); ++i)
79  kernel[i] *= a;
80 
81  cerr << "Kernel (" << kernel.size() << "): ";
82  copy(kernel.begin(), kernel.end(), ostream_iterator<double>(cerr, ","));
83  cerr << endl;
84  cerr << "normalization = " << sum << endl;
85 
87  cin >> grid;
88  gridConvolve(grid, kernel);
89 
90  grid.addMetadata(hdr);
91  cout << grid;
92 }
93 
94 
95 
STL namespace.
std::string invocationHeader(int argc, char *argv[])
Create an invocation header.
Definition: utils.cpp:124
Namespace for Density package.
Definition: DensityGrid.hpp:48
void gridConvolve(DensityGrid< T > &grid, DensityGrid< T > &kernel)
Convolve a grid with another grid (kernel)
Definition: GridUtils.hpp:202
Namespace for most things not already encapsulated within a class.