Line data Source code
1 : //# dImageHistograms.cc: This program generates histograms from images
2 : //#
3 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id: dImageHistograms.cc 20567 2009-04-09 23:12:39Z gervandiepen $
28 : //
29 : // IMHIST iterates through an image and creates and displays histograms from
30 : // data chunks specified by the axes keyword.
31 : //
32 : //
33 : // in The name on the input aips++ image. Currently must be of type <Float>
34 : // and can be of any dimension.
35 : //
36 : // There is no default.
37 : //
38 : // axes An array of integers (1 relative) specifying which axes are to be
39 : // histogrammed and which ones the histograms will be displayed as
40 : // a function of. For example, axes=3 would cause imhist to work
41 : // out histograms of the third (z) axis and display them as a function of
42 : // location along the first (x) and second (y) axes. axes=1,3 would cause
43 : // histograms of 1-3 (x-z) planes to be made and then to be displayed as
44 : // a function of the second axis (y) location.
45 : //
46 : // The default is to make a histogram of the entire image.
47 : //
48 : //
49 : // blc,trc Region (1 relative)
50 : // inc Increment to step through image
51 : //
52 : // nbins This specifies the number of bins, which is the same for each
53 : // histogram. Note that the bin width is worked out for each histogram
54 : // separately from the data minimum and maximum for that data chunk
55 : // (unless you set the range keyword).
56 : //
57 : // The default is 25
58 : //
59 : // include This specifies a range of pixel values to *include* in generating
60 : // the histograms. If two values are given then include pixels in
61 : // the range include(1) to include(2). If one value is given, then
62 : // include pixels in the range -abs(include) to abs(include)
63 : // This range applies to all histograms.
64 : //
65 : // The default is to include all data.
66 : //
67 : // gauss If true, a Gaussian overlay with the same mean, sigma
68 : // (of the pixels that were binned) and integral (of the histogram)
69 : // will be drawn on each histogram.
70 : //
71 : // The default is true.
72 : //
73 : // cumu If true, a cumulative histogram is drawn.
74 : //
75 : // The defaults is false.
76 : //
77 : // log If true, the log of the histogram values is drawn.
78 : //
79 : // The default is false.
80 : //
81 : // list If true list some statistical information about each histogram
82 : //
83 : // The default is false.
84 : //
85 : // plotter The PGPLOT device.
86 : //
87 : // The default is /xs
88 : //
89 : // nxy The number of subplots to put on each page in x and y. Each
90 : // histogram takes one subplot.
91 : //
92 : // The default is 1,1
93 : //
94 : // To do:
95 : // add keyword exclude ? Bit hard as requires restucturing of min/max
96 : // storage image
97 : //
98 : #include <casacore/casa/aips.h>
99 : #include <casacore/casa/Arrays.h>
100 : #include <casacore/casa/Exceptions/Error.h>
101 : #include <casacore/casa/Inputs/Input.h>
102 : #include <casacore/casa/Logging.h>
103 : #include <casacore/casa/BasicSL/String.h>
104 :
105 : #include <casacore/images/Images/PagedImage.h>
106 : #include <casacore/images/Images/SubImage.h>
107 : #include <casacore/images/Regions/ImageRegion.h>
108 : #include <imageanalysis/ImageAnalysis/ImageHistograms.h>
109 : #include <casacore/lattices/LRegions/LCSlicer.h>
110 : #include <casacore/lattices/LRegions/LCBox.h>
111 : #include <casacore/casa/System/PGPlotter.h>
112 : #include <casacore/casa/OS/EnvVar.h>
113 :
114 : #include <casacore/casa/iostream.h>
115 :
116 :
117 : #include <casacore/casa/namespace.h>
118 : enum defaults {AXES, REGION, RANGE, NDEFAULTS};
119 :
120 :
121 1 : int main (int argc, const char* argv[])
122 : {
123 : try {
124 :
125 1 : Input inputs(1);
126 1 : inputs.version ("$Revision: 20567 $");
127 :
128 :
129 : // Get inputs
130 1 : String casadata = EnvironmentVariable::get("CASADATA");
131 1 : if (casadata.empty()) {
132 0 : cerr << "CASADATA env variable not defined. Can't find fixtures. Did you source the casainit.(c)sh file?" << endl;
133 0 : return 1;
134 : }
135 3 : String *parts = new String[2];
136 1 : split(casadata, parts, 2, String(" "));
137 1 : String datadir = parts[0] + "/unittest/image/";
138 4 : delete [] parts;
139 :
140 1 : String name = datadir + "test_image.im";
141 1 : inputs.create("in", name, "Input image name");
142 1 : inputs.create("axes", "-10", "Cursor axes");
143 1 : inputs.create("blc", "-10", "blc");
144 1 : inputs.create("trc", "-10", "trc");
145 1 : inputs.create("inc", "-10", "inc");
146 1 : inputs.create("nbins", "25", "Number of bins");
147 1 : inputs.create("include", "0.0", "Pixel range to include");
148 1 : inputs.create("gauss", "true", "Plot Gaussian equivalent ?");
149 1 : inputs.create("cumu", "false", "Plot cumulative histogram ?");
150 1 : inputs.create("log", "false", "Take log of y axis ?");
151 1 : inputs.create("list", "false", "List statistics for each histogram");
152 1 : inputs.create("plotter", "", "Plot device");
153 1 : inputs.create("nxy", "1,1", "Number of subplots in x & y");
154 1 : inputs.create("disk", "false", "Force storage image to be disk based");
155 1 : inputs.readArguments(argc, argv);
156 :
157 1 : const String in = inputs.getString("in");
158 1 : const Block<Int> cursorAxesB(inputs.getIntArray("axes"));
159 1 : const Block<Int> blcB(inputs.getIntArray("blc"));
160 1 : const Block<Int> trcB(inputs.getIntArray("trc"));
161 1 : const Block<Int> incB(inputs.getIntArray("inc"));
162 1 : const Int nBins = inputs.getInt("nbins");
163 1 : const Bool doGauss = inputs.getBool("gauss");
164 1 : const Bool doCumu = inputs.getBool("cumu");
165 1 : const Bool doLog = inputs.getBool("log");
166 1 : const Block<Double> includeB = inputs.getDoubleArray("include");
167 1 : const Bool doList = inputs.getBool("list");
168 1 : const Block<Int> nxyB(inputs.getIntArray("nxy"));
169 1 : const String device = inputs.getString("plotter");
170 1 : const Bool force = inputs.getBool("disk");
171 :
172 :
173 : // Create defaults array
174 :
175 1 : Vector<Bool> validInputs(NDEFAULTS);
176 1 : validInputs = false;
177 2 : LogOrigin lor("imhist", "main()", WHERE);
178 1 : LogIO os(lor);
179 :
180 : // Check inputs
181 :
182 1 : if (in.empty()) {
183 0 : os << LogIO::SEVERE << "You must give an input image" << LogIO::POST;
184 0 : return 1;
185 : }
186 :
187 : // Convert cursor axes array to a vector (0 relative)
188 :
189 1 : Vector<Int> cursorAxes;
190 1 : if (cursorAxesB.nelements() == 1 && cursorAxesB[0] == -10) {
191 1 : cursorAxes.resize(0);
192 : } else {
193 0 : cursorAxes.resize(cursorAxesB.nelements());
194 0 : for (uInt i=0; i<cursorAxes.nelements(); i++) cursorAxes(i) = cursorAxesB[i] - 1;
195 0 : validInputs(AXES) = true;
196 : }
197 :
198 : // Convert region things to IPositions (0 relative)
199 :
200 1 : IPosition blc;
201 1 : IPosition trc;
202 1 : IPosition inc;
203 1 : if (blcB.nelements() == 1 && blcB[0] == -10) {
204 1 : blc.resize(0);
205 : } else {
206 0 : blc.resize(blcB.nelements());
207 0 : for (uInt i=0; i<blcB.nelements(); i++) blc(i) = blcB[i] - 1;
208 0 : validInputs(REGION) = true;
209 : }
210 1 : if (trcB.nelements() == 1 && trcB[0] == -10) {
211 1 : trc.resize(0);
212 : } else {
213 0 : trc.resize(trcB.nelements());
214 0 : for (uInt i=0; i<trcB.nelements(); i++) trc(i) = trcB[i] - 1;
215 0 : validInputs(REGION) = true;
216 : }
217 1 : if (incB.nelements() == 1 && incB[0] == -10) {
218 1 : inc.resize(0);
219 : } else {
220 0 : inc.resize(incB.nelements());
221 0 : for (uInt i=0; i<incB.nelements(); i++) inc(i) = incB[i];
222 0 : validInputs(REGION) = true;
223 : }
224 :
225 : // Convert inclusion range to vector
226 :
227 1 : Vector<Float> include(includeB.nelements());
228 : uInt i;
229 2 : for (i=0;i<include.nelements(); i++) {
230 1 : include(i) = includeB[i];
231 : }
232 1 : if (include.nelements() == 1 && include(0)==0) {
233 1 : include.resize(0);
234 : } else {
235 0 : validInputs(RANGE) = true;
236 : }
237 :
238 :
239 : // Plotting things
240 :
241 1 : Vector<Int> nxy;
242 1 : if (nxyB.nelements() == 1 && nxyB[0] == -1) {
243 0 : nxy.resize(0);
244 : } else {
245 1 : nxy.resize(nxyB.nelements());
246 3 : for (i=0;i<nxy.nelements(); i++) nxy(i) = nxyB[i];
247 : }
248 :
249 :
250 : // Do the work
251 :
252 1 : os << LogIO::NORMAL << "Using image " << in << LogIO::POST;
253 1 : DataType imageType = imagePixelType(in);
254 1 : if (imageType==TpFloat) {
255 :
256 : // Construct image
257 :
258 1 : PagedImage<Float> inImage(in);
259 1 : SubImage<Float>* pSubImage2 = 0;
260 :
261 1 : if (validInputs(REGION)) {
262 0 : LCBox::verify(blc, trc, inc, inImage.shape());
263 0 : cout << "Selected region : " << blc+1<< " to "
264 0 : << trc+1 << endl;
265 0 : const LCSlicer region(blc, trc);
266 : //
267 0 : SubImage<Float>* pSubImage = 0;
268 0 : if (inImage.isMasked()) {
269 : ImageRegion mask =
270 0 : inImage.getRegion(inImage.getDefaultMask(),
271 0 : RegionHandler::Masks);
272 0 : pSubImage = new SubImage<Float>(inImage, mask);
273 : } else {
274 0 : pSubImage = new SubImage<Float>(inImage);
275 : }
276 0 : pSubImage2 = new SubImage<Float>(*pSubImage, ImageRegion(region));
277 0 : delete pSubImage;
278 : } else {
279 1 : if (inImage.isMasked()) {
280 : ImageRegion mask =
281 0 : inImage.getRegion(inImage.getDefaultMask(),
282 0 : RegionHandler::Masks);
283 0 : pSubImage2 = new SubImage<Float>(inImage, mask);
284 : } else {
285 1 : pSubImage2 = new SubImage<Float>(inImage);
286 : }
287 : }
288 :
289 : // Construct histogram object
290 :
291 1 : casa::ImageHistograms<Float> histo(*pSubImage2, os, true, force);
292 1 : if (pSubImage2!=0) delete pSubImage2;
293 :
294 :
295 : // Set state
296 :
297 1 : if (validInputs(AXES)) {
298 0 : if (!histo.setAxes(cursorAxes)) {
299 0 : os << histo.errorMessage() << LogIO::POST;
300 0 : return 1;
301 : }
302 : }
303 1 : if (!histo.setNBins(nBins)) {
304 0 : os << histo.errorMessage() << LogIO::POST;
305 0 : return 1;
306 : }
307 1 : if (validInputs(RANGE)) {
308 0 : if (!histo.setIncludeRange(include)) {
309 0 : os << histo.errorMessage() << LogIO::POST;
310 0 : return 1;
311 : }
312 : }
313 1 : if (!histo.setGaussian (doGauss)) {
314 0 : os << histo.errorMessage() << LogIO::POST;
315 0 : return 1;
316 : }
317 1 : if (!histo.setForm(doLog, doCumu)) {
318 0 : os << histo.errorMessage() << LogIO::POST;
319 0 : return 1;
320 : }
321 1 : if (!histo.setStatsList(doList)) {
322 0 : os << histo.errorMessage() << LogIO::POST;
323 0 : return 1;
324 : }
325 1 : if (!device.empty()) {
326 0 : PGPlotter plotter(device);
327 0 : if (!histo.setPlotting(plotter, nxy)) {
328 0 : os << histo.errorMessage() << LogIO::POST;
329 0 : return 1;
330 : }
331 : }
332 :
333 : // Display histograms
334 :
335 1 : if (!histo.display()) {
336 0 : os << histo.errorMessage() << LogIO::POST;
337 0 : return 1;
338 : }
339 :
340 : // Get histograms
341 :
342 1 : Array<Float> values, counts;
343 1 : os << LogIO::NORMAL << "Recovering histogram arrays" << LogIO::POST;
344 1 : if (!histo.getHistograms(values,counts)) {
345 0 : os << histo.errorMessage() << LogIO::POST;
346 0 : return 1;
347 : }
348 : // cout << "values=" << values << endl;
349 : // cout << "counts=" << counts << endl;
350 :
351 1 : os << LogIO::NORMAL << "Recovering individual histogram arrays" << LogIO::POST;
352 1 : Vector<Float> valuesV, countsV;
353 1 : IPosition pos(histo.displayAxes().nelements(),0);
354 1 : if (!histo.getHistogram(valuesV,countsV,pos,false)) {
355 0 : os << histo.errorMessage() << LogIO::POST;
356 0 : return 1;
357 : }
358 :
359 : // cout << "values=" << valuesV << endl;
360 : // cout << "counts=" << countsV << endl;
361 :
362 :
363 : // Test copy constructor
364 :
365 1 : os << LogIO::NORMAL << "Applying copy constructor" << endl;
366 1 : casa::ImageHistograms<Float> histo2(histo);
367 :
368 : // Test assignment operator
369 :
370 1 : os << "Applying assignment operator" << LogIO::POST;
371 1 : histo = histo2;
372 :
373 : // Test setNewImage
374 :
375 1 : os << "Test setNewImage" << LogIO::POST;
376 1 : if (!histo.setNewImage(inImage)) {
377 0 : os << histo.errorMessage() << LogIO::POST;
378 0 : return 1;
379 : }
380 1 : if (!histo.display()) {
381 0 : os << histo.errorMessage() << LogIO::POST;
382 0 : return 1;
383 : }
384 :
385 : } else {
386 : os << "images of type " << Int(imageType)
387 0 : << " not yet supported" << endl;
388 0 : return 1;
389 : }
390 : }
391 0 : catch (AipsError x) {
392 0 : cerr << "aipserror: error " << x.getMesg() << endl;
393 0 : return 1;
394 : }
395 :
396 1 : return 0;
397 :
398 : }
|
