1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2015-2017 Imperial College London
5 * Copyright 2015-2017 Andreas Schuh
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19
20 #include "mirtk/Common.h"
21 #include "mirtk/Options.h"
22
23 #include "mirtk/ImageConfig.h"
24 #include "mirtk/IOConfig.h"
25
26 #include "mirtk/DataOp.h"
27 #include "mirtk/DataStatistics.h"
28 #include "mirtk/DataFunctions.h"
29
30 #if MIRTK_Image_WITH_VTK
31 #include "vtkDataSet.h"
32 #include "vtkSmartPointer.h"
33 #include "vtkPointData.h"
34 #include "vtkCellData.h"
35 #include "vtkDataArray.h"
36 #endif
37
38 using namespace mirtk;
39 using namespace mirtk::data;
40 using namespace mirtk::data::op;
41 using namespace mirtk::data::statistic;
42
43
44 // =============================================================================
45 // Help
46 // =============================================================================
47
48 // -----------------------------------------------------------------------------
PrintHelp(const char * name)49 void PrintHelp(const char *name)
50 {
51 cout << "\n";
52 cout << "Usage: " << name << " <input> [options]\n";
53 cout << "\n";
54 cout << "Description:\n";
55 cout << " This tool can be used for basic calculations from a sequence of data values read\n";
56 cout << " either from an image or a VTK pointset. It can be used, for example, to add two\n";
57 cout << " data sequences and to divide the result by a constant. The current sequence can\n";
58 cout << " be written to an output file again using :option:`-out`. Additionally, statistics\n";
59 cout << " of the current data sequence can be computed such as the mean or variance.\n";
60 cout << " The order of the data transformations and calculation of statistics is determined\n";
61 cout << " by the order of the command-line arguments.\n";
62 cout << "\n";
63 cout << " The data mask is used to include/exclude values from subsequent operations.\n";
64 cout << " Initially, all NaN values in the input data sequence are excluded.\n";
65 cout << " Further values can be excluded using one or more of the masking operations.\n";
66 cout << " Using the mask, operations can be performed on only a subset of the data,\n";
67 cout << " and the mask then reset using :option:`-reset-mask`.\n";
68 cout << "\n";
69 cout << " By default, data statistics are printed to STDOUT in a human readable format.\n";
70 cout << " This output can be appended to a text file using :option:`-append` instead.\n";
71 cout << " For a more machine readable output, e.g., as comma separated values (CSV),\n";
72 cout << " specify a delimiting string using :option:`-delimiter`. In this case, a header\n";
73 cout << " line is also printed when :option:`-header` is given with optional user\n";
74 cout << " specified column names for the individual output values.\n";
75 cout << "\n";
76 cout << "Input options:\n";
77 cout << " -pd, -point-data, -scalars <name> Name of input point data array. (default: active SCALARS array)\n";
78 cout << " -cd, -cell-data <name> Name of input cell data array. Overrides :option:`-pd`.\n";
79 cout << "\n";
80 cout << "Data masking options:\n";
81 cout << " -even\n";
82 cout << " Exclude values which are not an even number when cast to an integer.\n";
83 cout << " -odd\n";
84 cout << " Exclude values which are not an odd number when cast to an integer.\n";
85 cout << " -label <value|lower..upper>...\n";
86 cout << " Include data points with a value equal to either one of the given values.\n";
87 cout << " Closed intervals of values can be specified as \"lower..upper\".\n";
88 cout << " For example, \"-label 1 3 5..6 10 20..50\". This option is a shorthand for\n";
89 cout << " :option:`-mask-all` :option:`-threshold-inside` <lower> <upper> :option:`-invert-mask`\n";
90 cout << " where one :option:`-threshold-inside` operation is performed for each argument.\n";
91 cout << " -mask <value>... | <file> [<scalars>] [<value>]\n";
92 cout << " Exclude values equal a given threshold or with specified input mask <value>.\n";
93 cout << " The default mask value of values to be excluded is zero. When the input file\n";
94 cout << " is a point set file (e.g., .vtk, .vtp), the optional <scalars> argument can be\n";
95 cout << " used to specify the name of the point/cell data array to use as mask.\n";
96 cout << " Note that this operation does not modify the data values, but only marks them\n";
97 cout << " to be ignored from now on. Use :option:`-pad` following this operation to\n";
98 cout << " replace these values by a constant background value.\n";
99 cout << " -mask-all\n";
100 cout << " Exclude all values.\n";
101 cout << " -reset-mask\n";
102 cout << " Reset mask to include all values again.\n";
103 cout << " -invert-mask\n";
104 cout << " Invert mask to include all values that where excluded before and\n";
105 cout << " exclude all values that were included before.\n";
106 cout << " -set, -inside <value>\n";
107 cout << " Set new value for all currently included data values.\n";
108 cout << " -pad, -outside <value>\n";
109 cout << " Set new value for all currently excluded data values.\n";
110 cout << "\n";
111 cout << "Data thresholding options:\n";
112 cout << " -threshold <lower> [<upper>]\n";
113 cout << " This masking operation is equivalent to :option:`-threshold-outside`.\n";
114 cout << " When no upper threshold is specified, it defaults to +inf. Therefore,\n";
115 cout << " \"-threshold 0\" will exclude all negative values.\n";
116 cout << " -percentile-threshold, -pct-threshold <lower>\n";
117 cout << " This masking operation is equivalent to :option:`-threshold-outside-percentiles`.\n";
118 cout << " with an upper threshold of +inf. Therefore, \"-threshold 0\" excludes all negative values.\n";
119 cout << " -threshold-percentiles, -threshold-pcts <lower> <upper>\n";
120 cout << " This masking operation is equivalent to :option:`-threshold-outside-percentiles`.\n";
121 cout << " -threshold-inside, -mask-inside <lower> <upper>\n";
122 cout << " Exclude values which are inside a given closed interval.\n";
123 cout << " When the lower threshold is greater than the upper threshold,\n";
124 cout << " values less than or equal to the upper threshold and values greater\n";
125 cout << " than or equal to the lower threshold are excluded.\n";
126 cout << " -threshold-inside-percentiles, -threshold-inside-pcts, -mask-inside-percentiles, -mask-inside-pct <lower> <upper>\n";
127 cout << " Exclude values which are inside a given closed interval of percentiles.\n";
128 cout << " When the lower percentile is greater than the upper percentile,\n";
129 cout << " values less than or equal to the upper percentile and values greater\n";
130 cout << " than or equal to the lower percentile are excluded.\n";
131 cout << " -threshold-outside, -mask-outside <lower> <upper>\n";
132 cout << " Exclude values which are outside a given open interval.\n";
133 cout << " When the lower threshold is greater than the upper threshold,\n";
134 cout << " values inside the closed interval <upper>..<lower> are excluded.\n";
135 cout << " -threshold-outside-percentiles, -threshold-outside-pcts, -mask-outside-percentiles, -mask-outside-pcts <lower> <upper>\n";
136 cout << " Exclude values which are outside a given open interval of percentiles.\n";
137 cout << " When the lower percentile is greater than the upper percentile,\n";
138 cout << " values inside the closed interval <upper>..<lower> are excluded.\n";
139 cout << " -threshold-lt, -lower-threshold, -mask-lt <value>\n";
140 cout << " Exclude values less than a given threshold.\n";
141 cout << " -threshold-lt-percentile, -threshold-lt-pct, -lower-percentile-threshold, -lower-pct-threshold, -mask-lt-percentile, -mask-lt-pct <value>\n";
142 cout << " Exclude values less than a given precentile.\n";
143 cout << " -threshold-le, -mask-le, -mask-below <value>\n";
144 cout << " Exclude values less than or equal to a given threshold.\n";
145 cout << " -threshold-le-percentile, -threshold-le-pct, -mask-le-percentile, -mask-le-pct, -mask-below-percentile, -mask-below-pct <value>\n";
146 cout << " Exclude values less than or equal to a given percentile.\n";
147 cout << " -threshold-ge, -mask-ge, -mask-above <value>\n";
148 cout << " Exclude values greater than or equal to a given threshold.\n";
149 cout << " -threshold-ge-percentile, -threshold-ge-pct, -mask-ge-percentile, -mask-ge-pct, -mask-above-percentile, -mask-above-pct <value>\n";
150 cout << " Exclude values greater than or equal to a given percentile.\n";
151 cout << " -threshold-gt, -upper-threshold, -mask-gt <value>\n";
152 cout << " Exclude values greater than a given threshold.\n";
153 cout << " -threshold-gt-percentile, -threshold-gt-pct, -upper-percentile-threshold, -upper-pct-threshold, -mask-gt-percentile, -mask-gt-pct <value>\n";
154 cout << " Exclude values greater than a given percentile.\n";
155 cout << "\n";
156 cout << "Data rescaling options:\n";
157 cout << " -binarize <lower> [<upper>]\n";
158 cout << " Set values inside the closed interval <lower>..<upper> to one,\n";
159 cout << " and all other values to zero. The default upper threshold is +inf.\n";
160 cout << " When the lower threshold is greater than the upper threshold,\n";
161 cout << " values inside the closed interval <upper>..<lower> are set to zero\n";
162 cout << " and all other values to one instead. This operation is short for:\n";
163 cout << " :option:`-threshold-inside` <lower> <upper> :option:`-set` 1 :option:`-pad` 0\n";
164 cout << " -clamp <lower> <upper>\n";
165 cout << " Clamp values which are less than a lower or greater than an upper threshold.\n";
166 cout << " -clamp-percentiles, -clamp-pcts <lower> <upper>\n";
167 cout << " Clamp values which are less than a lower percentile or greater than an upper percentile.\n";
168 cout << " -clamp-below, -clamp-lt <value>\n";
169 cout << " Clamp values less than a given threshold.\n";
170 cout << " -clamp-below-percentile, -clamp-below-pct, -clamp-lt-percentile, -clamp-lt-pct <value>\n";
171 cout << " Clamp values less than a given percentile.\n";
172 cout << " -clamp-above, -clamp-gt <value>\n";
173 cout << " Clamp values greater than a given threshold.\n";
174 cout << " -clamp-above-percentile, -clamp-above-pct, -clamp-gt-percentile, -clamp-gt-pct <value>\n";
175 cout << " Clamp values greater than a given percentile.\n";
176 cout << " -rescale <min> <max>\n";
177 cout << " Linearly rescale values to the interval [min, max].\n";
178 cout << " -map <from> <to>...\n";
179 cout << " Replaces values equal to <from> by the specified <to> value. Multiple pairs of <from>\n";
180 cout << " and <to> value replacements can be specified in order to perform the substitutions in\n";
181 cout << " one step. For example, to swap the two values 1 and 2, use ``-map 1 2 2 1``.\n";
182 cout << "\n";
183 cout << "Arithmetic operation options:\n";
184 cout << " -add, -plus <value> | <file> [<scalars>]\n";
185 cout << " Add constant value or data sequence read from specified file.\n";
186 cout << " Another name for this option is the '+' sign, see Examples.\n";
187 cout << " -sub, -subtract, -minus <value> | <file> [<scalars>]\n";
188 cout << " Subtract constant value or data sequence read from specified file.\n";
189 cout << " Another name for this option is the '-' sign, see Examples.\n";
190 cout << " -mul, -multiply-with, -times <value> | <file> [<scalars>]\n";
191 cout << " Multiply by constant value or data sequence read from specified file.\n";
192 cout << " Another name for this option is the '*' sign, see Examples.\n";
193 cout << " -div, -divide-by, -over <value> | sum | <file> [<scalars>]\n";
194 cout << " Divide by constant value or data sequence read from specified file.\n";
195 cout << " When the argument is \"sum\", the divisor is the sum of the values.\n";
196 cout << " When dividing by zero values in the input file, the result is NaN.\n";
197 cout << " Use :option:`-mask` with argument NaN and :option:`-pad` to replace\n";
198 cout << " these undefined values by a constant such as zero.\n";
199 cout << " Another name for this option is the '/' sign, see Examples.\n";
200 cout << " -div-with-zero <value> | sum | <file> [<scalars>]\n";
201 cout << " Same as :option:`-div`, but set result to zero in case of division by zero.\n";
202 cout << " -abs\n";
203 cout << " Replace values by their respective absolute value.\n";
204 cout << " -pow, -power <exponent>\n";
205 cout << " Raise values to the power of the given exponent.\n";
206 cout << " -sq, -square\n";
207 cout << " Raise values to the power of 2 (i.e, -pow 2).\n";
208 cout << " -sqrt\n";
209 cout << " Calculate square root of each value (i.e, -pow .5).\n";
210 cout << " -exp\n";
211 cout << " Calculate exponential of data sequence.\n";
212 cout << " -log [<threshold>] [<base>]\n";
213 cout << " Compute logarithm after applying an optional threshold.\n";
214 cout << " (default threshold: min double, default base: e)\n";
215 cout << " -lb, -log2 [<threshold>]\n";
216 cout << " Compute binary logarithm, alias for :option:`-log` with base 2.\n";
217 cout << " -ln, -loge [<threshold>]\n";
218 cout << " Compute natural logarithm, alias for :option:`-log` with base e.\n";
219 cout << " -lg, -log10 [<threshold>]\n";
220 cout << " Compute logarithm to base 10, alias for :option:`-log` with base 10.\n";
221 cout << " -mod, -fmod <denominator>\n";
222 cout << " Compute modulo division of each value with specified denominator.\n";
223 cout << " -floor\n";
224 cout << " Round floating point values to largest integer value that is not greater.\n";
225 cout << " -ceil\n";
226 cout << " Round floating point values to smallest integer value that is greater.\n";
227 cout << " -round\n";
228 cout << " Round floating point values to the nearest integer value, away from zero for halfway cases.\n";
229 cout << "\n";
230 cout << "Data output options:\n";
231 cout << " -out, -o, -output <file> [<type>] [<name>]\n";
232 cout << " Write current data sequence to file in the format of the input file.\n";
233 cout << " Output data type can be: uchar, short, ushort, int, uint, float, double.\n";
234 cout << " The optional <name> argument can be used to save the modified data\n";
235 cout << " of an input point set data array with a different name along with the\n";
236 cout << " input data. Otherwise, the input data values are replaced by the modified\n";
237 cout << " values and stored with point data array name is unchanged.\n";
238 cout << " Another name for this option is the '=' sign, but the optional arguments are\n";
239 cout << " are not supported by this alternative notation. See Examples for usage.\n";
240 cout << "\n";
241 cout << "Data statistics options:\n";
242 cout << " -append <file>\n";
243 cout << " Append output to a file. (default: STDOUT)\n";
244 cout << " -delimiter, -delim, -d, -sep\n";
245 cout << " Delimiting character(s). (default: '')\n";
246 cout << " -header [<name>...]\n";
247 cout << " Request output of header line if delimiter was specified as well.\n";
248 cout << " If the output is appended to a text file, the header is only printed\n";
249 cout << " if it does not exist. If no or fewer custom column names are given,\n";
250 cout << " the default names for each statistic are printed. (default: none)\n";
251 cout << " -prefix <str>...\n";
252 cout << " One or more prefix strings to print. If no delimiter is specified,\n";
253 cout << " the concatenated strings are printed before each line of the output.\n";
254 cout << " Otherwise, each prefix string is printed as entry for the first columns\n";
255 cout << " in the delimited output row, separated by the specified delimiter. (default: none)\n";
256 cout << " -precision, -digits <int>\n";
257 cout << " Number of significant digits. (default: 5)\n";
258 cout << " -median\n";
259 cout << " Print median value, i.e., 50th percentile. (default: off)\n";
260 cout << " -mean, -avg, -average\n";
261 cout << " Print mean value. (default: on)\n";
262 cout << " -variance, -var\n";
263 cout << " Print variance of values. (default: off)\n";
264 cout << " -sigma, -std, -stddev, -stdev, -sd\n";
265 cout << " Print standard deviation of values. (default: on)\n";
266 cout << " -normal-distribution\n";
267 cout << " Print mean and standard deviation of values.\n";
268 cout << " Other option names: -mean+sigma, -mean+sd, -avg+std,... (default: off)\n";
269 cout << " -mad, -mean-absolute-difference, -mean-absolute-deviation\n";
270 cout << " Print mean absolute difference/deviation around the mean. (default: off)\n";
271 cout << " -mad-median, -median-absolute-difference, -median-absolute-deviation\n";
272 cout << " Print mean absolute difference/deviation around the median. (default: off)\n";
273 cout << " -minimum, -min\n";
274 cout << " Print minimum value. (default: off)\n";
275 cout << " -maximum, -max\n";
276 cout << " Print maximum value. (default: off)\n";
277 cout << " -extrema, -minmax\n";
278 cout << " Print minimum and maximum value. (default: on)\n";
279 cout << " -range\n";
280 cout << " Print range of values (i.e., max - min). (default: off)\n";
281 cout << " -percentile, -pct, -p <n>...\n";
282 cout << " Print n-th percentile. (default: none)\n";
283 cout << " -lower-percentile-mean, -lpctavg <n>\n";
284 cout << " Print mean intensity of values less than or equal to the n-th percentile. (default: off)\n";
285 cout << " -upper-percentile-mean, -upctavg <n>\n";
286 cout << " Print mean intensity of values greater than or equal to the n-th percentile. (default: off)\n";
287 cout << " -sum\n";
288 cout << " Print sum of values. Can be used to count values within a certain range using a thresholding\n";
289 cout << " followed by :option:`-set` 1 before summing these values. (default: off)\n";
290 cout << " -count\n";
291 cout << " Print number of values inside the mask, i.e., values not currently excluded. (default: off)\n";
292 PrintCommonOptions(cout);
293 cout << "\n";
294 cout << "Examples:\n";
295 cout << "\n";
296 cout << " " << name << " mni305.nii.gz\n";
297 cout << " Mean = 26.9753\n";
298 cout << " Standard deviation = 50.3525\n";
299 cout << " Extrema = [0, 254]\n";
300 cout << " Range = 254\n";
301 cout << "\n";
302 cout << " " << name << " mni305.nii.gz -pct 77\n";
303 cout << " 77th percentile = 25\n";
304 cout << "\n";
305 cout << " " << name << " mni305.nii.gz -padding 25 -range -percentile 25 50 75 -prefix MNI305 '[>25]'\n";
306 cout << " MNI305 [>25] range = 254\n";
307 cout << " MNI305 [>25] 25th percentile = 69\n";
308 cout << " MNI305 [>25] 50th percentile = 113\n";
309 cout << " MNI305 [>25] 75th percentile = 150\n";
310 cout << "\n";
311 cout << " " << name << " mni305.nii.gz -d , -prefix MNI305\n";
312 cout << " MNI305,26.9753,50.3525,0,254,254 [no newline at end of line]\n";
313 cout << "\n";
314 cout << " " << name << " mni305.nii.gz -d , -prefix MNI305 -header\n";
315 cout << " ,Mean,Sigma,Min,Max,Range\n";
316 cout << " MNI305,26.9753,50.3525,0,254,254\n";
317 cout << "\n";
318 cout << " " << name << " mni305.nii.gz -d , -prefix MNI305 -header ID Mean SD\n";
319 cout << " ID,Mean,SD,Min,Max,Range\n";
320 cout << " MNI305,26.9753,50.3525,0,254,254\n";
321 cout << "\n";
322 cout << " " << name << " a.nii.gz + b.nii.gz = c.nii.gz\n";
323 cout << "\n";
324 cout << " " << name << " a.vtk + b.nii.gz - 10 / c.nii = d.vtk\n";
325 cout << " Adds data values at identical sequential memory indices in a and b,\n";
326 cout << " subtracts the constant 10, and then divides by the values in image c.\n";
327 cout << "\n";
328 cout << " Note: Operations are always executed from left to right,\n";
329 cout << " i.e., no mathematical operator precedence is considered!\n";
330 cout << "\n";
331 }
332
333 // =============================================================================
334 // Main
335 // =============================================================================
336
337 // -----------------------------------------------------------------------------
338 // Some special options do not start with a '-' as otherwise required
339 #undef HAS_ARGUMENT
340 #define HAS_ARGUMENT \
341 _IsArgument(ARGIDX, argc, argv) && \
342 strcmp(argv[ARGIDX+1], "+") != 0 && \
343 strcmp(argv[ARGIDX+1], "/") != 0 && \
344 strcmp(argv[ARGIDX+1], "=") != 0
345
346 // -----------------------------------------------------------------------------
main(int argc,char ** argv)347 int main(int argc, char **argv)
348 {
349 InitializeIOLibrary();
350
351 // Initial data values
352 REQUIRES_POSARGS(1);
353
354 const char *input_name = POSARG(1);
355
356 UniquePtr<double[]> data;
357 int datatype = MIRTK_VOXEL_DOUBLE;
358 ImageAttributes attr;
359
360 #if MIRTK_Image_WITH_VTK
361 const char *scalars_name = nullptr;
362 bool cell_data = false;
363 for (ARGUMENTS_AFTER(1)) {
364 if (OPTION("-point-data") || OPTION("-pointdata") || OPTION("-pd") || OPTION("-scalars")) {
365 scalars_name = ARGUMENT;
366 cell_data = false;
367 }
368 else if (OPTION("-cell-data") || OPTION("-celldata") || OPTION("-cd")) {
369 scalars_name = ARGUMENT;
370 cell_data = true;
371 }
372 }
373 vtkSmartPointer<vtkDataSet> dataset;
374 vtkSmartPointer<vtkDataSetAttributes> arrays;
375 int n = Read(input_name, data, &datatype, &attr, &dataset, scalars_name, cell_data);
376 if (dataset) {
377 if (cell_data) {
378 arrays = dataset->GetCellData();
379 } else {
380 arrays = dataset->GetPointData();
381 }
382 }
383 #else // MIRTK_Image_WITH_VTK
384 int n = Read(input_name, data, &datatype, &attr);
385 #endif // MIRTK_Image_WITH_VTK
386
387 // Optional arguments
388 const double inf = numeric_limits<double>::infinity();
389 const double nan = numeric_limits<double>::quiet_NaN();
390
391 double a, b;
392 int p;
393 const char *append_name = NULL;
394 const char *delimiter = NULL;
395 bool print_header = false;
396 int digits = 5;
397
398 Array<string> header;
399 Array<string> prefix;
400
401 Array<UniquePtr<Op> > ops;
402
403 for (ARGUMENTS_AFTER(1)) {
404 if (OPTION("-append")) {
405 append_name = ARGUMENT;
406 } else if (OPTION("-point-data") || OPTION("-pointdata") || OPTION("-pd") || OPTION("-scalars")) {
407 #if MIRTK_Image_WITH_VTK
408 // Parsed before Read above
409 scalars_name = ARGUMENT;
410 cell_data = false;
411 #else
412 FatalError("Cannot process -point-data of VTK file because MIRTK Image library was built without VTK!");
413 #endif // MIRTK_Image_WITH_VTK
414 } else if (OPTION("-cell-data") || OPTION("-celldata") || OPTION("-cd")) {
415 #if MIRTK_Image_WITH_VTK
416 // Parsed before Read above
417 scalars_name = ARGUMENT;
418 cell_data = true;
419 #else
420 FatalError("Cannot process -cell-data of VTK file because MIRTK Image library was built without VTK!");
421 #endif // MIRTK_Image_WITH_VTK
422 } else if (OPTION("-prefix")) {
423 do {
424 prefix.push_back(ARGUMENT);
425 } while (HAS_ARGUMENT);
426 } else if (OPTION("-header")) {
427 print_header = true;
428 while (HAS_ARGUMENT) header.push_back(ARGUMENT);
429 // Masking
430 } else if (OPTION("-label")) {
431 ops.push_back(UniquePtr<Op>(new ResetMask(true)));
432 do {
433 const char *arg = ARGUMENT;
434 const Array<string> parts = Split(arg, "..");
435 if (parts.size() == 1) {
436 if (!FromString(parts[0], a)) a = nan;
437 b = a;
438 } else if (parts.size() == 2) {
439 if (!FromString(parts[0], a) || !FromString(parts[1], b)) {
440 a = b = nan;
441 }
442 } else {
443 a = b = nan;
444 }
445 if (IsNaN(a) || IsNaN(b)) {
446 FatalError("Invalid -label argument: " << arg);
447 }
448 ops.push_back(UniquePtr<Op>(new MaskInsideInterval(a, b)));
449 } while (HAS_ARGUMENT);
450 ops.push_back(UniquePtr<Op>(new InvertMask()));
451 } else if (OPTION("-mask-all")) {
452 ops.push_back(UniquePtr<Op>(new ResetMask(false)));
453 } else if (OPTION("-reset-mask")) {
454 ops.push_back(UniquePtr<Op>(new ResetMask(true)));
455 } else if (OPTION("-invert-mask")) {
456 ops.push_back(UniquePtr<Op>(new InvertMask()));
457 } else if (OPTION("-mask")) {
458 double c;
459 do {
460 const char *arg = ARGUMENT;
461 if (FromString(arg, c)) {
462 ops.push_back(UniquePtr<Op>(new Mask(c)));
463 } else {
464 const char *fname = arg;
465 const char *aname = nullptr;
466 if (HAS_ARGUMENT) {
467 arg = ARGUMENT;
468 if (HAS_ARGUMENT) {
469 aname = arg;
470 PARSE_ARGUMENT(c);
471 } else if (!FromString(arg, c)) {
472 aname = arg, c = 0.;
473 }
474 } else {
475 c = 0.;
476 #if MIRTK_Image_WITH_VTK
477 if (dataset && arrays->HasArray(fname)) {
478 aname = fname;
479 fname = input_name;
480 }
481 #endif
482 }
483 UniquePtr<Mask> op(new Mask(fname, c));
484 if (aname) {
485 #if MIRTK_Image_WITH_VTK
486 op->ArrayName(aname);
487 op->IsCellData(cell_data);
488 #else
489 FatalError("Cannot read point set files when build without VTK or wrong usage!");
490 #endif
491 }
492 ops.push_back(UniquePtr<Op>(op.release()));
493 break;
494 }
495 } while (HAS_ARGUMENT);
496 } else if (OPTION("-threshold-outside") || OPTION("-mask-outside")) {
497 PARSE_ARGUMENT(a);
498 PARSE_ARGUMENT(b);
499 ops.push_back(UniquePtr<Op>(new MaskOutsideOpenInterval(a, b)));
500 } else if (OPTION("-threshold-outside-percentiles") || OPTION("-threshold-outside-pcts") ||
501 OPTION("-mask-outside-percentiles") || OPTION("-mask-outside-pcts")) {
502 PARSE_ARGUMENT(p);
503 Statistic *a = new Percentile(p);
504 a->Hidden(verbose < 1);
505 ops.push_back(UniquePtr<Op>(a));
506 PARSE_ARGUMENT(p);
507 Statistic *b = new Percentile(p);
508 b->Hidden(verbose < 1);
509 ops.push_back(UniquePtr<Op>(b));
510 Op *op = new MaskOutsideOpenInterval(&a->Value(), &b->Value());
511 ops.push_back(UniquePtr<Op>(op));
512 } else if (OPTION("-threshold")) {
513 PARSE_ARGUMENT(a);
514 if (HAS_ARGUMENT) PARSE_ARGUMENT(b);
515 else b = inf;
516 ops.push_back(UniquePtr<Op>(new MaskOutsideInterval(a, b)));
517 } else if (OPTION("-percentile-threshold") || OPTION("-pct-threshold")) {
518 PARSE_ARGUMENT(p);
519 Statistic *a = new Percentile(p);
520 a->Hidden(verbose < 1);
521 ops.push_back(UniquePtr<Op>(a));
522 Op *op = new MaskOutsideInterval(&a->Value(), inf);
523 ops.push_back(UniquePtr<Op>(op));
524 } else if (OPTION("-threshold-percentiles") || OPTION("-threshold-pcts")) {
525 PARSE_ARGUMENT(p);
526 Statistic *a = new Percentile(p);
527 a->Hidden(verbose < 1);
528 ops.push_back(UniquePtr<Op>(a));
529 PARSE_ARGUMENT(p);
530 Statistic *b = new Percentile(p);
531 b->Hidden(verbose < 1);
532 ops.push_back(UniquePtr<Op>(b));
533 Op *op = new MaskOutsideInterval(&a->Value(), &b->Value());
534 ops.push_back(UniquePtr<Op>(op));
535 } else if (OPTION("-threshold-inside") || OPTION("-mask-inside")) {
536 PARSE_ARGUMENT(a);
537 PARSE_ARGUMENT(b);
538 ops.push_back(UniquePtr<Op>(new MaskInsideInterval(a, b)));
539 } else if (OPTION("-threshold-inside-percentiles") || OPTION("-threshold-inside-pcts") ||
540 OPTION("-mask-inside-percentiles") || OPTION("-mask-inside-pcts")) {
541 PARSE_ARGUMENT(p);
542 Statistic *a = new Percentile(p);
543 a->Hidden(verbose < 1);
544 ops.push_back(UniquePtr<Op>(a));
545 PARSE_ARGUMENT(p);
546 Statistic *b = new Percentile(p);
547 b->Hidden(verbose < 1);
548 ops.push_back(UniquePtr<Op>(b));
549 Op *op = new MaskInsideInterval(&a->Value(), &b->Value());
550 ops.push_back(UniquePtr<Op>(op));
551 } else if (OPTION("-threshold-lt") || OPTION("-lower-threshold") || OPTION("-mask-lt")) {
552 PARSE_ARGUMENT(a);
553 ops.push_back(UniquePtr<Op>(new MaskOutsideInterval(a, inf)));
554 } else if (OPTION("-threshold-lt-percentile") || OPTION("-threshold-lt-pct") ||
555 OPTION("-lower-percentile-threshold") || OPTION("-lower-pct-threshold") ||
556 OPTION("-mask-lt-percentile") || OPTION("-mask-lt-pct")) {
557 PARSE_ARGUMENT(p);
558 Statistic *a = new Percentile(p);
559 a->Hidden(verbose < 1);
560 ops.push_back(UniquePtr<Op>(a));
561 ops.push_back(UniquePtr<Op>(new MaskOutsideInterval(&a->Value(), inf)));
562 } else if (OPTION("-threshold-le") || OPTION("-mask-below") || OPTION("-mask-le")) {
563 PARSE_ARGUMENT(a);
564 ops.push_back(UniquePtr<Op>(new MaskOutsideOpenInterval(a, inf)));
565 } else if (OPTION("-threshold-le-percentile") || OPTION("-threshold-le-pct") ||
566 OPTION("-mask-below-percentile") || OPTION("-mask-below-pct") ||
567 OPTION("-mask-le-percentile") || OPTION("-mask-le-pct")) {
568 PARSE_ARGUMENT(p);
569 Statistic *a = new Percentile(p);
570 a->Hidden(verbose < 1);
571 ops.push_back(UniquePtr<Op>(a));
572 ops.push_back(UniquePtr<Op>(new MaskOutsideOpenInterval(&a->Value(), inf)));
573 } else if (OPTION("-threshold-ge") || OPTION("-mask-above") || OPTION("-mask-ge")) {
574 PARSE_ARGUMENT(b);
575 ops.push_back(UniquePtr<Op>(new MaskOutsideOpenInterval(-inf, b)));
576 } else if (OPTION("-threshold-ge-percentile") || OPTION("-threshold-ge-pct") ||
577 OPTION("-mask-above-percentile") || OPTION("-mask-above-pct") ||
578 OPTION("-mask-ge-percentile") || OPTION("-mask-ge-pct")) {
579 PARSE_ARGUMENT(p);
580 Statistic *b = new Percentile(p);
581 b->Hidden(verbose < 1);
582 ops.push_back(UniquePtr<Op>(b));
583 ops.push_back(UniquePtr<Op>(new MaskOutsideOpenInterval(-inf, &b->Value())));
584 } else if (OPTION("-threshold-gt") || OPTION("-upper-threshold") || OPTION("-mask-gt")) {
585 PARSE_ARGUMENT(b);
586 ops.push_back(UniquePtr<Op>(new MaskOutsideInterval(-inf, b)));
587 } else if (OPTION("-threshold-gt-percentile") || OPTION("-threshold-gt-pct") ||
588 OPTION("-upper-percentile-threshold") || OPTION("-upper-pct-threshold") ||
589 OPTION("-mask-gt-percentile") || OPTION("-mask-gt-pct")) {
590 PARSE_ARGUMENT(p);
591 Statistic *b = new Percentile(p);
592 b->Hidden(verbose < 1);
593 ops.push_back(UniquePtr<Op>(b));
594 ops.push_back(UniquePtr<Op>(new MaskOutsideInterval(-inf, &b->Value())));
595 } else if (OPTION("-even")) {
596 ops.push_back(UniquePtr<Op>(new MaskOddValues()));
597 } else if (OPTION("-odd")) {
598 ops.push_back(UniquePtr<Op>(new MaskEvenValues()));
599 // Clamping
600 } else if (OPTION("-clamp")) {
601 PARSE_ARGUMENT(a);
602 PARSE_ARGUMENT(b);
603 ops.push_back(UniquePtr<Op>(new Clamp(a, b)));
604 } else if (OPTION("-clamp-percentiles") || OPTION("-clamp-pcts")) {
605 PARSE_ARGUMENT(p);
606 Statistic *a = new Percentile(p);
607 a->Hidden(verbose < 1);
608 ops.push_back(UniquePtr<Op>(a));
609 PARSE_ARGUMENT(p);
610 Statistic *b = new Percentile(p);
611 b->Hidden(verbose < 1);
612 ops.push_back(UniquePtr<Op>(b));
613 ops.push_back(UniquePtr<Op>(new Clamp(&a->Value(), &b->Value())));
614 } else if (OPTION("-clamp-lt") || OPTION("-clamp-below")) {
615 PARSE_ARGUMENT(a);
616 ops.push_back(UniquePtr<Op>(new LowerThreshold(a)));
617 } else if (OPTION("-clamp-lt-percentile") || OPTION("-clamp-lt-pct") ||
618 OPTION("-clamp-below-percentile") || OPTION("-clamp-below-pct")) {
619 PARSE_ARGUMENT(p);
620 Statistic *a = new Percentile(p);
621 a->Hidden(verbose < 1);
622 ops.push_back(UniquePtr<Op>(a));
623 ops.push_back(UniquePtr<Op>(new LowerThreshold(&a->Value())));
624 } else if (OPTION("-clamp-gt") || OPTION("-clamp-above")) {
625 PARSE_ARGUMENT(b);
626 ops.push_back(UniquePtr<Op>(new UpperThreshold(b)));
627 } else if (OPTION("-clamp-gt-percentile") || OPTION("-clamp-gt-pct") ||
628 OPTION("-clamp-above-percentile") || OPTION("-clamp-above-pct")) {
629 PARSE_ARGUMENT(p);
630 Statistic *b = new Percentile(p);
631 b->Hidden(verbose < 1);
632 ops.push_back(UniquePtr<Op>(b));
633 ops.push_back(UniquePtr<Op>(new UpperThreshold(&b->Value())));
634 } else if (OPTION("-rescale")) {
635 double min, max;
636 if (!FromString(ARGUMENT, min)) {
637 cerr << "Invalid -rescale minimum, must be a number!" << endl;
638 exit(1);
639 }
640 if (!FromString(ARGUMENT, max)) {
641 cerr << "Invalid -rescale maximum, must be a number!" << endl;
642 exit(1);
643 }
644 ops.push_back(UniquePtr<Op>(new Rescale(min, max)));
645 } else if (OPTION("-set") || OPTION("-inside")) {
646 double inside_value;
647 if (!FromString(ARGUMENT, inside_value)) {
648 cerr << "Invalid -inside value, must be a number!" << endl;
649 exit(1);
650 }
651 ops.push_back(UniquePtr<Op>(new SetInsideValue(inside_value)));
652 } else if (OPTION("-pad") || OPTION("-outside")) {
653 double outside_value;
654 if (!FromString(ARGUMENT, outside_value)) {
655 cerr << "Invalid -outside value, must be a number!" << endl;
656 exit(1);
657 }
658 ops.push_back(UniquePtr<Op>(new SetOutsideValue(outside_value)));
659 // Data transformations
660 } else if (OPTION("-binarize")) {
661 PARSE_ARGUMENT(a);
662 if (HAS_ARGUMENT) PARSE_ARGUMENT(b);
663 else b = inf;
664 ops.push_back(UniquePtr<Op>(new Binarize(a, b)));
665 } else if (OPTION("-map")) {
666 UniquePtr<Map> map(new Map());
667 do {
668 const char * const arg1 = ARGUMENT;
669 const char * const arg2 = ARGUMENT;
670 if (!FromString(arg1, a) || !FromString(arg2, b)) {
671 FatalError("Arguments of -map option must be pairs of two numbers (i.e., number of arguments must be even)!");
672 }
673 map->Insert(a, b);
674 } while (HAS_ARGUMENT);
675 ops.push_back(UniquePtr<Op>(map.release()));
676 } else if (OPTION("-add") || OPTION("-plus") || OPTION("+")) {
677 const char *arg = ARGUMENT;
678 double c;
679 if (FromString(arg, c)) {
680 ops.push_back(UniquePtr<Op>(new Add(c)));
681 } else {
682 const char *fname = arg;
683 const char *aname = nullptr;
684 if (HAS_ARGUMENT) {
685 aname = ARGUMENT;
686 } else {
687 #if MIRTK_Image_WITH_VTK
688 if (dataset && arrays->HasArray(fname)) {
689 aname = fname;
690 fname = input_name;
691 }
692 #endif
693 }
694 UniquePtr<Add> op(new Add(fname));
695 if (aname) {
696 #if MIRTK_Image_WITH_VTK
697 op->ArrayName(aname);
698 op->IsCellData(cell_data);
699 #else
700 FatalError("Cannot read scalars from point set file when build without VTK or wrong usage!");
701 #endif
702 }
703 ops.push_back(UniquePtr<Op>(op.release()));
704 }
705 } else if (OPTION("-sub") || OPTION("-subtract") || OPTION("-minus") || OPTION("-")) {
706 const char *arg = ARGUMENT;
707 double c;
708 if (FromString(arg, c)) {
709 ops.push_back(UniquePtr<Op>(new Sub(c)));
710 } else {
711 const char *fname = arg;
712 const char *aname = nullptr;
713 if (HAS_ARGUMENT) {
714 aname = ARGUMENT;
715 } else {
716 #if MIRTK_Image_WITH_VTK
717 if (dataset && arrays->HasArray(fname)) {
718 aname = fname;
719 fname = input_name;
720 }
721 #endif
722 }
723 UniquePtr<Sub> op(new Sub(fname));
724 if (aname) {
725 #if MIRTK_Image_WITH_VTK
726 op->ArrayName(aname);
727 op->IsCellData(cell_data);
728 #else
729 FatalError("Cannot read point set files when build without VTK or wrong usage!");
730 #endif
731 }
732 ops.push_back(UniquePtr<Op>(op.release()));
733 }
734 } else if (OPTION("-mul") || OPTION("-multiply-by") || OPTION("-times") || OPTION("*")) {
735 const char *arg = ARGUMENT;
736 double c;
737 if (FromString(arg, c)) {
738 ops.push_back(UniquePtr<Op>(new Mul(c)));
739 } else {
740 const char *fname = arg;
741 const char *aname = nullptr;
742 if (HAS_ARGUMENT) {
743 aname = ARGUMENT;
744 } else {
745 #if MIRTK_Image_WITH_VTK
746 if (dataset && arrays->HasArray(fname)) {
747 aname = fname;
748 fname = input_name;
749 }
750 #endif
751 }
752 UniquePtr<Mul> op(new Mul(fname));
753 if (aname) {
754 #if MIRTK_Image_WITH_VTK
755 op->ArrayName(aname);
756 op->IsCellData(cell_data);
757 #else
758 FatalError("Cannot read point set files when build without VTK or wrong usage!");
759 #endif
760 }
761 ops.push_back(UniquePtr<Op>(op.release()));
762 }
763 } else if (OPTION("-div") || OPTION("-divide-by") || OPTION("-over") || OPTION("/")) {
764 const char *arg = ARGUMENT;
765 double c;
766 if (ToLower(arg) == "sum") {
767 Statistic *a = new Sum();
768 a->Hidden(verbose < 1);
769 ops.push_back(UniquePtr<Op>(a));
770 ops.push_back(UniquePtr<Op>(new Div(&a->Value())));
771 } else if (FromString(arg, c)) {
772 if (fequal(c, .0)) {
773 cerr << "Invalid -div argument, value must not be zero!" << endl;
774 exit(1);
775 }
776 ops.push_back(UniquePtr<Op>(new Div(c)));
777 } else {
778 const char *fname = arg;
779 const char *aname = nullptr;
780 if (HAS_ARGUMENT) {
781 aname = ARGUMENT;
782 } else {
783 #if MIRTK_Image_WITH_VTK
784 if (dataset && arrays->HasArray(fname)) {
785 aname = fname;
786 fname = input_name;
787 }
788 #endif
789 }
790 UniquePtr<Div> op(new Div(fname));
791 if (aname) {
792 #if MIRTK_Image_WITH_VTK
793 op->ArrayName(aname);
794 op->IsCellData(cell_data);
795 #else
796 FatalError("Cannot read point set files when build without VTK or wrong usage!");
797 #endif
798 }
799 ops.push_back(UniquePtr<Op>(op.release()));
800 }
801 } else if (OPTION("-div-with-zero")) {
802 const char *arg = ARGUMENT;
803 double c;
804 if (ToLower(arg) == "sum") {
805 Statistic *a = new Sum();
806 a->Hidden(verbose < 1);
807 ops.push_back(UniquePtr<Op>(a));
808 ops.push_back(UniquePtr<Op>(new DivWithZero(&a->Value())));
809 } else if (FromString(arg, c)) {
810 ops.push_back(UniquePtr<Op>(new DivWithZero(c)));
811 } else {
812 const char *fname = arg;
813 const char *aname = nullptr;
814 if (HAS_ARGUMENT) {
815 aname = ARGUMENT;
816 } else {
817 #if MIRTK_Image_WITH_VTK
818 if (dataset && arrays->HasArray(fname)) {
819 aname = fname;
820 fname = input_name;
821 }
822 #endif
823 }
824 UniquePtr<DivWithZero> op(new DivWithZero(fname));
825 if (aname) {
826 #if MIRTK_Image_WITH_VTK
827 op->ArrayName(aname);
828 op->IsCellData(cell_data);
829 #else
830 FatalError("Cannot read point set files when build without VTK or wrong usage!");
831 #endif
832 }
833 ops.push_back(UniquePtr<Op>(op.release()));
834 }
835 } else if (OPTION("-abs")) {
836 ops.push_back(UniquePtr<Op>(new Abs()));
837 } else if (OPTION("-pow") || OPTION("-power")) {
838 const char *arg = ARGUMENT;
839 double exponent;
840 if (!FromString(arg, exponent)) {
841 cerr << "Invalid -power value, must be a number!" << endl;
842 exit(1);
843 }
844 ops.push_back(UniquePtr<Op>(new Pow(exponent)));
845 } else if (OPTION("-sqrt")) {
846 ops.push_back(UniquePtr<Op>(new Pow(.5)));
847 } else if (OPTION("-square") || OPTION("-sq")) {
848 ops.push_back(UniquePtr<Op>(new Pow(2.0)));
849 } else if (OPTION("-exp")) {
850 ops.push_back(UniquePtr<Op>(new Exp()));
851 } else if (OPTION("-log") || OPTION("-log2") || OPTION("-loge") || OPTION("-log10") || OPTION("-lb") || OPTION("-ln") || OPTION("-lg")) {
852 a = numeric_limits<double>::min();
853 if (HAS_ARGUMENT) {
854 PARSE_ARGUMENT(a);
855 if (a <= .0) {
856 cerr << "Invalid -log threshold argument, must be a positive number" << endl;
857 exit(1);
858 }
859 }
860 Op *op = nullptr;
861 if (strcmp(OPTNAME, "-log") == 0) {
862 if (HAS_ARGUMENT) {
863 double base;
864 if (!FromString(ARGUMENT, base)) {
865 char c;
866 if (!FromString(ARGUMENT, c) || c != 'e') {
867 cerr << "Invalid -log base argument, must be a positive number or character e" << endl;
868 exit(1);
869 }
870 op = new Ln(a);
871 } else {
872 op = new Log(base, a);
873 }
874 } else {
875 op = new Ln(a);
876 }
877 } else if (strcmp(OPTNAME, "-log2") == 0 || strcmp(OPTNAME, "-lb") == 0) {
878 op = new Lb(a);
879 } else if (strcmp(OPTNAME, "-log10") == 0 || strcmp(OPTNAME, "-lg") == 0) {
880 op = new Lg(a);
881 } else if (strcmp(OPTNAME, "-loge") == 0 || strcmp(OPTNAME, "-ln") == 0) {
882 op = new Ln(a);
883 }
884 ops.push_back(UniquePtr<Op>(op));
885 } else if (OPTION("-mod") || OPTION("-fmod")) {
886 const char *arg = ARGUMENT;
887 double denominator;
888 if (!FromString(arg, denominator) || abs(denominator) < 1e-12) {
889 cerr << "Invalid -mod value, must be a non-zero number!" << endl;
890 exit(1);
891 }
892 ops.push_back(UniquePtr<Op>(new Mod(denominator)));
893 } else if (OPTION("-floor")) {
894 ops.push_back(UniquePtr<Op>(new Floor()));
895 } else if (OPTION("-ceil")) {
896 ops.push_back(UniquePtr<Op>(new Ceil()));
897 } else if (OPTION("-round")) {
898 ops.push_back(UniquePtr<Op>(new Round()));
899 } else if (OPTION("=")) {
900 const char *fname = ARGUMENT;
901 #if MIRTK_Image_WITH_VTK
902 ops.push_back(UniquePtr<Op>(new Write(fname, datatype, attr, dataset, scalars_name, scalars_name)));
903 #else
904 ops.push_back(UniquePtr<Op>(new Write(fname, datatype, attr)));
905 #endif
906 } else if (OPTION("-o") || OPTION("-out") || OPTION("-output")) {
907 const char *fname = ARGUMENT;
908 int dtype = datatype;
909 #if MIRTK_Image_WITH_VTK
910 const char *output_scalars_name = scalars_name;
911 #endif
912 if (HAS_ARGUMENT) {
913 const char *arg = ARGUMENT;
914 dtype = ToDataType(arg);
915 if (dtype == MIRTK_VOXEL_UNKNOWN) {
916 cerr << "Invalid -out data type " << arg << endl;
917 exit(1);
918 }
919 if (HAS_ARGUMENT) {
920 #if MIRTK_Image_WITH_VTK
921 output_scalars_name = ARGUMENT;
922 #else
923 Warning("Output scalars array name argument of -output option ignored");
924 #endif
925 }
926 }
927 #if MIRTK_Image_WITH_VTK
928 ops.push_back(UniquePtr<Op>(new Write(fname, dtype, attr, dataset, scalars_name, output_scalars_name, cell_data)));
929 #else
930 ops.push_back(UniquePtr<Op>(new Write(fname, dtype, attr)));
931 #endif
932 // Data statistics
933 } else if (OPTION("-median")) {
934 ops.push_back(UniquePtr<Op>(new Median()));
935 } else if (OPTION("-mean") || OPTION("-average") || OPTION("-avg")) {
936 ops.push_back(UniquePtr<Op>(new Mean()));
937 } else if (OPTION("-sigma") || OPTION("-stddev") || OPTION("-stdev") || OPTION("-std") || OPTION("-sd")) {
938 ops.push_back(UniquePtr<Op>(new StDev()));
939 } else if (OPTION("-normal-distribution") ||
940 OPTION("-mean+sigma") || OPTION("-mean+stddev") || OPTION("-mean+stdev") || OPTION("-mean+std") || OPTION("-mean+sd") ||
941 OPTION("-avg+sigma") || OPTION("-avg+stddev") || OPTION("-avg+stdev") || OPTION("-avg+std") || OPTION("-avg+sd")) {
942 ops.push_back(UniquePtr<Op>(new NormalDistribution()));
943 } else if (OPTION("-variance") || OPTION("-var")) {
944 ops.push_back(UniquePtr<Op>(new Var()));
945 } else if (OPTION("-mean-absolute-difference") || OPTION("-mean-absolute-deviation") || OPTION("-mad") || OPTION("-mad-mean")) {
946 ops.push_back(UniquePtr<Op>(new MeanAbsoluteDifference()));
947 } else if (OPTION("-median-absolute-difference") || OPTION("-median-absolute-deviation") || OPTION("-mad-median")) {
948 ops.push_back(UniquePtr<Op>(new MedianAbsoluteDifference()));
949 } else if (OPTION("-minimum") || OPTION("-min")) {
950 ops.push_back(UniquePtr<Op>(new Min()));
951 } else if (OPTION("-maximum") || OPTION("-max")) {
952 ops.push_back(UniquePtr<Op>(new Max()));
953 } else if (OPTION("-extrema") || OPTION("-minmax")) {
954 ops.push_back(UniquePtr<Op>(new Extrema()));
955 } else if (OPTION("-range")) {
956 ops.push_back(UniquePtr<Op>(new Range()));
957 } else if (OPTION("-percentile") || OPTION("-pct") || OPTION("-p")) {
958 do {
959 int p;
960 if (FromString(ARGUMENT, p) && 0 <= p && p <= 100) {
961 ops.push_back(UniquePtr<Op>(new Percentile(p)));
962 } else {
963 cerr << "Invalid -percentile value, must be integer in the range [0, 100]!" << endl;
964 exit(1);
965 }
966 } while (HAS_ARGUMENT);
967 } else if (OPTION("-lower-percentile-mean") || OPTION("-lpctavg")) {
968 do {
969 int p;
970 if (FromString(ARGUMENT, p) && 0 <= p && p <= 100) {
971 ops.push_back(UniquePtr<Op>(new LowerPercentileMean(p)));
972 } else {
973 cerr << "Invalid -lower-percentile-mean value, must be integer in the range [0, 100]!" << endl;
974 exit(1);
975 }
976 } while (HAS_ARGUMENT);
977 } else if (OPTION("-upper-percentile-mean") || OPTION("-upctavg")) {
978 do {
979 int p;
980 if (FromString(ARGUMENT, p) && 0 <= p && p <= 100) {
981 ops.push_back(UniquePtr<Op>(new UpperPercentileMean(p)));
982 } else {
983 cerr << "Invalid -upper-percentile-mean value, must be integer in the range [0, 100]!" << endl;
984 exit(1);
985 }
986 } while (HAS_ARGUMENT);
987 } else if (OPTION("-sum")) {
988 ops.push_back(UniquePtr<Op>(new Sum()));
989 } else if (OPTION("-count")) {
990 ops.push_back(UniquePtr<Op>(new Count()));
991 } else if (OPTION("-delimiter") || OPTION("-delim") || OPTION("-d") || OPTION("-sep")) {
992 delimiter = ARGUMENT;
993 } else if (OPTION("-precision") || OPTION("-digits")) {
994 if (!FromString(ARGUMENT, digits) || digits < 0) {
995 cerr << "Invalid -precision argument, value must be non-negative integer!" << endl;
996 exit(1);
997 }
998 } else {
999 HANDLE_COMMON_OR_UNKNOWN_OPTION();
1000 }
1001 }
1002
1003 // If delimiter explicitly set to empty string, use none
1004 if (delimiter && delimiter[0] == '\0') delimiter = NULL;
1005
1006 // Default statistics to compute
1007 if (ops.empty()) {
1008 ops.push_back(UniquePtr<Statistic>(new Mean()));
1009 ops.push_back(UniquePtr<Statistic>(new StDev()));
1010 ops.push_back(UniquePtr<Statistic>(new Extrema()));
1011 ops.push_back(UniquePtr<Statistic>(new Range()));
1012 }
1013
1014 // Initial data mask
1015 UniquePtr<bool[]> mask(new bool[n]);
1016 for (int i = 0; i < n; ++i) {
1017 if (IsNaN(data[i])) {
1018 mask[i] = false;
1019 } else {
1020 mask[i] = true;
1021 }
1022 }
1023
1024 // Process input data, either transform it or compute statistics from it
1025 for (size_t i = 0; i < ops.size(); ++i) {
1026 ops[i]->Process(n, data.get(), mask.get());
1027 }
1028 mask.reset();
1029
1030 // Open output file to append to or use STDOUT if none specified
1031 ofstream ofs;
1032 if (append_name) {
1033 if (print_header) {
1034 ifstream ifs(append_name);
1035 if (ifs.is_open()) {
1036 print_header = false;
1037 ifs.close();
1038 }
1039 }
1040 ofs.open(append_name, ios_base::app);
1041 if (!ofs.is_open()) {
1042 FatalError("Cannot append to file " << append_name);
1043 }
1044 }
1045 ostream &out = (ofs.is_open() ? ofs : cout);
1046
1047 // Print column names if requested
1048 if (delimiter && print_header) {
1049 size_t c = 0;
1050 for (size_t i = 0; i < prefix.size(); ++i, ++c) {
1051 if (c > 0) out << delimiter;
1052 if (c < header.size()) out << header[c];
1053 }
1054 for (size_t i = 0; i < ops.size(); ++i) {
1055 Statistic *stat = dynamic_cast<Statistic *>(ops[i].get());
1056 if (stat != nullptr && !stat->Hidden()) {
1057 for (size_t j = 0; j < stat->Names().size(); ++j, ++c) {
1058 if (c > 0) out << delimiter;
1059 if (c < header.size()) out << header[c];
1060 else out << stat->Names()[j];
1061 }
1062 }
1063 }
1064 out << endl;
1065 }
1066
1067 // Print image statistics
1068 if (delimiter) {
1069 for (size_t i = 0; i < prefix.size(); ++i) {
1070 if (i > 0) out << delimiter;
1071 out << prefix[i];
1072 }
1073 bool first = prefix.empty();
1074 for (size_t i = 0; i < ops.size(); ++i) {
1075 Statistic *stat = dynamic_cast<Statistic *>(ops[i].get());
1076 if (stat != nullptr && !stat->Hidden() && !stat->Names().empty()) {
1077 if (!first) out << delimiter;
1078 else first = false;
1079 stat->PrintValues(out, digits, delimiter);
1080 }
1081 }
1082 // No newline at end of row if printing results to STDOUT which in this
1083 // case is usually assigned to a string in a calling script
1084 if (print_header || ofs.is_open()) out << endl;
1085 } else {
1086 string prefix_string;
1087 for (size_t i = 0; i < prefix.size(); ++i) {
1088 if (i > 0) prefix_string += ' ';
1089 prefix_string += prefix[i];
1090 }
1091 for (size_t i = 0; i < ops.size(); ++i) {
1092 Statistic *stat = dynamic_cast<Statistic *>(ops[i].get());
1093 if (stat != nullptr && !stat->Hidden()) {
1094 stat->Print(out, digits, prefix_string.c_str());
1095 }
1096 }
1097 }
1098
1099 ofs.close();
1100 return 0;
1101 }
1102