LOOS  v2.3.2
gridavg.cpp
1 /*
2  gridavg
3 
4  Average grids together. Requires that grids have the same dimensions.
5 */
6 
7 
8 /*
9  This file is part of LOOS.
10 
11  LOOS (Lightweight Object-Oriented Structure library)
12  Copyright (c) 2013, 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 <3) {
43  cerr <<
44  "DESCRIPTION\n\tAverage together multiple grids\n"
45  "\nUSAGE\n\tgridavg grid1 grid2 [grid3 ...] >averaged.grid\n"
46  "Requires that the grids have the same dimensions.\n"
47  "\nEXAMPLES\n\tgridavg water1.grid water2.grid water3.grid >water.grid\n";
48  exit(0);
49  }
50 
51  string hdr = invocationHeader(argc, argv);
52 
53 
55  int k = 1;
56  ifstream ifs(argv[k]);
57  if (ifs.fail()) {
58  cerr << "Error- cannot open " << argv[k] << endl;
59  exit(-1);
60  }
61 
62  ifs >> avg;
63  avg.addMetadata(hdr);
64  ifs.close();
65 
66  uint n = 0;
67  for (++k;k < argc; ++k) {
68  ifs.open(argv[k]);
69  if (ifs.fail()) {
70  cerr << "Error- cannot open " << argv[k] << endl;
71  exit(-2);
72  }
74  ifs >> grid;
75  if (grid.gridDims() != avg.gridDims()) {
76  cerr << "Error- grid in " << argv[k] << " has dimensions " << grid.gridDims() << ",\n"
77  << "but was expecting it to be " << avg.gridDims() << endl;
78  exit(-3);
79  }
80  if (grid.minCoord() != avg.minCoord() || grid.maxCoord() != avg.maxCoord())
81  cerr << "Warning- real world bounds for grid in " << argv[k] << " do not match. Proceeding anyway...\n";
82 
83  double mean = 0.0;
84 
85  for (long i=0; i<grid.size(); ++i) {
86  mean += grid(i);
87  avg(i) += grid(i);
88  }
89 
90  ++n;
91  }
92 
93  for (long i=0; i<avg.size(); ++i)
94  avg(i) /= n;
95 
96  cout << avg;
97 }
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
Namespace for most things not already encapsulated within a class.