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