1 /*
2  * Medical Image Registration ToolKit (MIRTK)
3  *
4  * Copyright 2013-2015 Imperial College London
5  * Copyright 2013-2015 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 
21 /// @file  mirtk/DataOp.h
22 /// @brief Defines base class and I/O functions for arbitrary 1D data sequences
23 ///
24 /// Functions to manipulate the data are defined in mirtkDataFunctions.h.
25 /// Statistics of the data sequence such as mean and variance or percentiles
26 /// can be computed using the operators found in mirtkDataStatistics.h.
27 /// The data operators are used in particular by the calculate.cc tool for
28 /// which they were originally developed. They were added to the linear
29 /// algebra library because they are useful to compute common statistics or
30 /// perform basic mathematical operations on a data sequence such as an image
31 /// or attributes of a VTK point set.
32 ///
33 /// @sa mirtkDataFunctions.h
34 /// @sa mirtkDataStatistics.h
35 #ifndef MIRTK_DataOp_H
36 #define MIRTK_DataOp_H
37 
38 #include "mirtk/ImageConfig.h"
39 
40 #include "mirtk/Math.h"
41 #include "mirtk/Stream.h"
42 #include "mirtk/Voxel.h"
43 #include "mirtk/ImageAttributes.h"
44 
45 #if MIRTK_Image_WITH_VTK
46 #  include "vtkSmartPointer.h"
47 #  include "vtkDataSet.h"
48 #  include "vtkDataArray.h"
49 #endif
50 
51 
52 namespace mirtk { namespace data {
53 
54 
55 // =============================================================================
56 // Base class of data operations
57 // =============================================================================
58 
59 // -----------------------------------------------------------------------------
60 /// Base class of all data operations
61 class Op
62 {
63 public:
64 
65   /// Destructor
~Op()66   virtual ~Op() {}
67 
68   /// Process given data
69   virtual void Process(int, double *, bool * = NULL) = 0;
70 
71 #if MIRTK_Image_WITH_VTK
72   /// Process given vtkDataArray
73   virtual void Process(vtkDataArray *data, bool *mask = nullptr)
74   {
75     const int n = static_cast<int>(data->GetNumberOfTuples() * data->GetNumberOfComponents());
76     if (data->GetDataType() == VTK_DOUBLE) {
77       this->Process(n, reinterpret_cast<double *>(data->GetVoidPointer(0)), mask);
78     } else {
79       UniquePtr<double[]> _data(new double[n]);
80       double *tuple = _data.get();
81       for (vtkIdType i = 0; i < data->GetNumberOfTuples(); ++i) {
82         data->GetTuple(i, tuple);
83         tuple += data->GetNumberOfComponents();
84       }
85       this->Process(n, _data.get(), mask);
86       tuple = _data.get();
87       for (vtkIdType i = 0; i < data->GetNumberOfTuples(); ++i) {
88         data->SetTuple(i, tuple);
89         tuple += data->GetNumberOfComponents();
90       }
91     }
92   }
93 #endif
94 };
95 
96 // =============================================================================
97 // I/O functions
98 // =============================================================================
99 
100 // -----------------------------------------------------------------------------
101 /// Enumeration of supported input data file types
102 enum DataFileType
103 {
104   UnknownDataFile,
105   IMAGE,
106   LEGACY_VTK,
107   XML_VTK
108 };
109 
110 /// Get (or guess) type of input file
111 DataFileType FileType(const char *name);
112 
113 // -----------------------------------------------------------------------------
114 /// Read data sequence from any supported input file type
115 int Read(const char *name,
116          UniquePtr<double[]> &data,
117          int *dtype = nullptr,
118          ImageAttributes *attr = nullptr,
119          #if MIRTK_Image_WITH_VTK
120            vtkSmartPointer<vtkDataSet> *dataset = nullptr,
121          #else
122            void * = nullptr,
123          #endif
124          const char *scalars_name = nullptr,
125          bool cell_data = false);
126 
127 // -----------------------------------------------------------------------------
128 /// Write data sequence
129 class Write : public Op
130 {
131   /// Name of output file
132   mirtkPublicAttributeMacro(string, FileName);
133 
134 #if MIRTK_Image_WITH_VTK
135 
136   /// VTK input dataset whose scalar data was modified
137   mirtkPublicAttributeMacro(vtkSmartPointer<vtkDataSet>, DataSet);
138 
139   /// Name of input data array
140   mirtkPublicAttributeMacro(string, ArrayName);
141 
142   /// Name of output data array
143   mirtkPublicAttributeMacro(string, OutputName);
144 
145   /// Whether to add data as cell data instead of point data
146   mirtkPublicAttributeMacro(bool, AsCellData);
147 
148 #endif // MIRTK_Image_WITH_VTK
149 
150   /// Attributes of input image whose data was modified
151   mirtkPublicAttributeMacro(ImageAttributes, Attributes);
152 
153   /// Output data type
154   mirtkPublicAttributeMacro(int, DataType);
155 
156 public:
157 
158   /// Constructor
159 #if MIRTK_Image_WITH_VTK
160 
161   Write(const char *fname, int dtype = MIRTK_VOXEL_DOUBLE,
162         ImageAttributes attr = ImageAttributes(),
163         vtkDataSet *dataset     = nullptr,
164         const char *array_name  = nullptr,
165         const char *output_name = nullptr,
166         bool        cell_data   = false)
167   :
_FileName(fname)168     _FileName(fname),
169     _DataSet(dataset),
170     _AsCellData(cell_data),
171     _Attributes(attr),
172     _DataType(dtype)
173   {
174     if (array_name)  _ArrayName  = array_name;
175     if (output_name) _OutputName = output_name;
176   }
177 
178 #else // MIRTK_Image_WITH_VTK
179 
180   Write(const char *fname, int dtype = MIRTK_VOXEL_DOUBLE,
181         ImageAttributes attr = ImageAttributes())
182   :
183     _FileName(fname),
184     _Attributes(attr),
185     _DataType(dtype)
186   {}
187 
188 #endif // MIRTK_Image_WITH_VTK
189 
190   /// Process given data
191   virtual void Process(int n, double *data, bool * = nullptr);
192 
193 };
194 
195 // =============================================================================
196 // Auxiliary macros for subclass implementation
197 // =============================================================================
198 
199 // -----------------------------------------------------------------------------
200 // Add Calculate function that takes a vtkDataArray as argument
201 // and computes a single return value
202 //
203 // Subclass must implement:
204 //   template <class T> double Calculate(int n, const T *data, const bool *mask)
205 #if MIRTK_Image_WITH_VTK
206   #define mirtkCalculateVtkDataArray1() \
207     static double Calculate(vtkDataArray *data, const bool *mask = NULL) \
208     { \
209       const void *p = data->GetVoidPointer(0); \
210       const int   n = static_cast<int>(data->GetNumberOfTuples() * \
211                                        data->GetNumberOfComponents()); \
212       switch (data->GetDataType()) { \
213         case VTK_SHORT:  return Calculate(n, reinterpret_cast<const short  *>(p), mask); \
214         case VTK_INT:    return Calculate(n, reinterpret_cast<const int    *>(p), mask); \
215         case VTK_FLOAT:  return Calculate(n, reinterpret_cast<const float  *>(p), mask); \
216         case VTK_DOUBLE: return Calculate(n, reinterpret_cast<const double *>(p), mask); \
217         default: \
218           cerr << "Unsupported vtkDataArray type: " << data->GetDataType() << endl; \
219           exit(1); \
220       } \
221       return numeric_limits<double>::quiet_NaN(); \
222     }
223 #else
224   #define mirtkCalculateVtkDataArray1()
225 #endif
226 
227 // -----------------------------------------------------------------------------
228 // Add Calculate function that takes a vtkDataArray as argument
229 // and computes two return values
230 //
231 // Subclass must implement:
232 //   template <class T>
233 //   void Calculate(double &, double &, int n, const T *data, const bool *mask)
234 #if MIRTK_Image_WITH_VTK
235   #define mirtkCalculateVtkDataArray2() \
236     static void Calculate(double &v1, double &v2, vtkDataArray *data, const bool *mask = NULL) \
237     { \
238       const void *p = data->GetVoidPointer(0); \
239       const int   n = static_cast<int>(data->GetNumberOfTuples() * \
240                                        data->GetNumberOfComponents()); \
241       switch (data->GetDataType()) { \
242         case VTK_SHORT:  Calculate(v1, v2, n, reinterpret_cast<const short  *>(p), mask); break; \
243         case VTK_INT:    Calculate(v1, v2, n, reinterpret_cast<const int    *>(p), mask); break; \
244         case VTK_FLOAT:  Calculate(v1, v2, n, reinterpret_cast<const float  *>(p), mask); break; \
245         case VTK_DOUBLE: Calculate(v1, v2, n, reinterpret_cast<const double *>(p), mask); break; \
246         default: \
247           cerr << "Unsupported vtkDataArray type: " << data->GetDataType() << endl; \
248           exit(1); \
249       } \
250     }
251 #else
252   #define mirtkCalculateVtkDataArray2()
253 #endif
254 
255 
256 } } // namespace mirtk::data
257 
258 #endif
259