Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // MOOSE includes
11 : #include "ImageSampler.h"
12 : #include "MooseApp.h"
13 : #include "ImageMesh.h"
14 :
15 : #include "libmesh/mesh_tools.h"
16 :
17 : using namespace libMesh;
18 :
19 : InputParameters
20 14932 : ImageSampler::validParams()
21 : {
22 : // Define the general parameters
23 14932 : InputParameters params = emptyInputParameters();
24 14932 : params += FileRangeBuilder::validParams();
25 :
26 14932 : params.addParam<Point>("origin", "Origin of the image (defaults to mesh origin)");
27 14932 : params.addParam<Point>("dimensions",
28 : "x,y,z dimensions of the image (defaults to mesh dimensions)");
29 14932 : params.addParam<unsigned int>(
30 : "component",
31 : "The image RGB-component to return, leaving this blank will result in a greyscale value "
32 : "for the image to be created. The component number is zero based, i.e. 0 returns the first "
33 : "(RED) component of the image.");
34 :
35 : // Shift and Scale (application of these occurs prior to threshold)
36 14932 : params.addParam<double>("shift", 0, "Value to add to all pixels; occurs prior to scaling");
37 44796 : params.addParam<double>(
38 29864 : "scale", 1, "Multiplier to apply to all pixel values; occurs after shifting");
39 14932 : params.addParamNamesToGroup("shift scale", "Rescale");
40 :
41 : // Threshold parameters
42 14932 : params.addParam<double>("threshold", "The threshold value");
43 44796 : params.addParam<double>(
44 29864 : "upper_value", 1, "The value to set for data greater than the threshold value");
45 44796 : params.addParam<double>(
46 29864 : "lower_value", 0, "The value to set for data less than the threshold value");
47 14932 : params.addParamNamesToGroup("threshold upper_value lower_value", "Threshold");
48 :
49 : // Flip image
50 14932 : params.addParam<bool>("flip_x", false, "Flip the image along the x-axis");
51 14932 : params.addParam<bool>("flip_y", false, "Flip the image along the y-axis");
52 14932 : params.addParam<bool>("flip_z", false, "Flip the image along the z-axis");
53 14932 : params.addParamNamesToGroup("flip_x flip_y flip_z", "Flip");
54 :
55 14932 : return params;
56 0 : }
57 :
58 349 : ImageSampler::ImageSampler(const InputParameters & parameters)
59 : : FileRangeBuilder(parameters),
60 : #ifdef LIBMESH_HAVE_VTK
61 349 : _data(nullptr),
62 349 : _algorithm(nullptr),
63 : #endif
64 349 : _is_pars(parameters),
65 349 : _is_console(
66 : (parameters.getCheckedPointerParam<MooseApp *>("_moose_app"))->getOutputWarehouse()),
67 349 : _flip({{_is_pars.get<bool>("flip_x"),
68 349 : _is_pars.get<bool>("flip_y"),
69 698 : _is_pars.get<bool>("flip_z")}})
70 :
71 : {
72 : #ifndef LIBMESH_HAVE_VTK
73 : // This should be impossible to reach, the registration of ImageSampler is also guarded with
74 : // LIBMESH_HAVE_VTK
75 : mooseError("libMesh must be configured with VTK enabled to utilize ImageSampler");
76 : #endif
77 349 : }
78 :
79 : void
80 346 : ImageSampler::setupImageSampler(MooseMesh & mesh)
81 : {
82 : // Don't warn that mesh or _is_pars are unused when VTK is not enabled.
83 346 : libmesh_ignore(mesh);
84 346 : libmesh_ignore(_is_pars);
85 :
86 : #ifdef LIBMESH_HAVE_VTK
87 : // Get access to the Mesh object
88 346 : BoundingBox bbox = MeshTools::create_bounding_box(mesh.getMesh());
89 :
90 : // Set the dimensions from the Mesh if not set by the User
91 346 : if (_is_pars.isParamValid("dimensions"))
92 104 : _physical_dims = _is_pars.get<Point>("dimensions");
93 :
94 : else
95 : {
96 242 : _physical_dims(0) = bbox.max()(0) - bbox.min()(0);
97 : #if LIBMESH_DIM > 1
98 242 : _physical_dims(1) = bbox.max()(1) - bbox.min()(1);
99 : #endif
100 : #if LIBMESH_DIM > 2
101 242 : _physical_dims(2) = bbox.max()(2) - bbox.min()(2);
102 : #endif
103 : }
104 :
105 : // Set the origin from the Mesh if not set in the input file
106 346 : if (_is_pars.isParamValid("origin"))
107 104 : _origin = _is_pars.get<Point>("origin");
108 : else
109 : {
110 242 : _origin(0) = bbox.min()(0);
111 : #if LIBMESH_DIM > 1
112 242 : _origin(1) = bbox.min()(1);
113 : #endif
114 : #if LIBMESH_DIM > 2
115 242 : _origin(2) = bbox.min()(2);
116 : #endif
117 : }
118 :
119 : // An array of filenames, to be filled in
120 346 : std::vector<std::string> filenames;
121 :
122 : // The file suffix, to be determined
123 346 : std::string file_suffix;
124 :
125 : // Try to parse our own file range parameters. If that fails, then
126 : // see if the associated Mesh is an ImageMesh and use its. If that
127 : // also fails, then we have to throw an error...
128 : //
129 : // The parseFileRange method sets parameters, thus a writable reference to the InputParameters
130 : // object must be obtained from the warehouse. Generally, this should be avoided, but
131 : // this is a special case.
132 346 : if (_status != 0)
133 : {
134 : // We don't have parameters, so see if we can get them from ImageMesh
135 26 : ImageMesh * image_mesh = dynamic_cast<ImageMesh *>(&mesh);
136 26 : if (!image_mesh)
137 0 : mooseError("No file range parameters were provided and the Mesh is not an ImageMesh.");
138 :
139 : // Get the ImageMesh's parameters. This should work, otherwise
140 : // errors would already have been thrown...
141 26 : filenames = image_mesh->filenames();
142 26 : file_suffix = image_mesh->fileSuffix();
143 : }
144 : else
145 : {
146 : // Use our own parameters (using 'this' b/c of conflicts with filenames the local variable)
147 320 : filenames = this->filenames();
148 320 : file_suffix = fileSuffix();
149 : }
150 :
151 : // Storage for the file names
152 346 : _files = vtkSmartPointer<vtkStringArray>::New();
153 :
154 1243 : for (const auto & filename : filenames)
155 897 : _files->InsertNextValue(filename);
156 :
157 : // Error if no files where located
158 346 : if (_files->GetNumberOfValues() == 0)
159 8 : mooseError("No image file(s) located");
160 :
161 : // Read the image stack. Hurray for VTK not using polymorphism in a
162 : // smart way... we actually have to explicitly create the type of
163 : // reader based on the file extension, using an if-statement...
164 338 : if (file_suffix == "png")
165 338 : _image = vtkSmartPointer<vtkPNGReader>::New();
166 0 : else if (file_suffix == "tiff" || file_suffix == "tif")
167 0 : _image = vtkSmartPointer<vtkTIFFReader>::New();
168 : else
169 0 : mooseError("Un-supported file type '", file_suffix, "'");
170 :
171 : // Now that _image is set up, actually read the images
172 : // Indicate that data read has started
173 338 : _is_console << "Reading image(s)..." << std::endl;
174 :
175 : // Extract the data
176 338 : _image->SetFileNames(_files);
177 338 : _image->Update();
178 338 : _data = _image->GetOutput();
179 338 : _algorithm = _image->GetOutputPort();
180 :
181 : // Set the image dimensions and voxel size member variable
182 338 : int * dims = _data->GetDimensions();
183 1352 : for (unsigned int i = 0; i < 3; ++i)
184 : {
185 1014 : _dims.push_back(dims[i]);
186 1014 : _voxel.push_back(_physical_dims(i) / _dims[i]);
187 : }
188 :
189 : // Set the dimensions of the image and bounding box
190 338 : _data->SetSpacing(_voxel[0], _voxel[1], _voxel[2]);
191 338 : _data->SetOrigin(_origin(0), _origin(1), _origin(2));
192 338 : _bounding_box.min() = _origin;
193 338 : _bounding_box.max() = _origin + _physical_dims;
194 :
195 : // Indicate data read is completed
196 338 : _is_console << " ...image read finished" << std::endl;
197 :
198 : // Set the component parameter
199 : // If the parameter is not set then vtkMagnitude() will applied
200 338 : if (_is_pars.isParamValid("component"))
201 : {
202 17 : unsigned int n = _data->GetNumberOfScalarComponents();
203 17 : _component = _is_pars.get<unsigned int>("component");
204 17 : if (_component >= n)
205 4 : mooseError("'component' parameter must be empty or have a value of 0 to ", n - 1);
206 : }
207 : else
208 321 : _component = 0;
209 :
210 : // Apply filters, the toggling on and off of each filter is handled internally
211 334 : vtkMagnitude();
212 334 : vtkShiftAndScale();
213 334 : vtkThreshold();
214 : #endif
215 334 : }
216 :
217 : Real
218 2022640 : ImageSampler::sample(const Point & p) const
219 : {
220 : #ifdef LIBMESH_HAVE_VTK
221 :
222 : // Do nothing if the point is outside of the image domain
223 2022640 : if (!_bounding_box.contains_point(p))
224 39648 : return 0.0;
225 :
226 : // Determine pixel coordinates
227 1982992 : std::vector<int> x(3, 0);
228 7931968 : for (int i = 0; i < LIBMESH_DIM; ++i)
229 : {
230 : // Compute position, only if voxel size is greater than zero
231 5948976 : if (_voxel[i] != 0)
232 : {
233 5143584 : x[i] = std::floor((p(i) - _origin(i)) / _voxel[i]);
234 :
235 : // If the point falls on the mesh extents the index needs to be decreased by one
236 5143584 : if (x[i] == _dims[i])
237 110640 : x[i]--;
238 :
239 : // flip
240 5143584 : if (_flip[i])
241 77440 : x[i] = _dims[i] - x[i] - 1;
242 : }
243 : }
244 :
245 : // Return the image data at the given point
246 1982992 : return _data->GetScalarComponentAsDouble(x[0], x[1], x[2], _component);
247 :
248 : #else
249 : libmesh_ignore(p); // avoid un-used parameter warnings
250 : return 0.0;
251 : #endif
252 1982992 : }
253 :
254 : void
255 334 : ImageSampler::vtkMagnitude()
256 : {
257 : #ifdef LIBMESH_HAVE_VTK
258 : // Do nothing if 'component' is set
259 334 : if (_is_pars.isParamValid("component"))
260 13 : return;
261 :
262 : // Apply the greyscale filtering
263 321 : _magnitude_filter = vtkSmartPointer<vtkImageMagnitude>::New();
264 321 : _magnitude_filter->SetInputConnection(_algorithm);
265 321 : _magnitude_filter->Update();
266 :
267 : // Update the pointers
268 321 : _data = _magnitude_filter->GetOutput();
269 321 : _algorithm = _magnitude_filter->GetOutputPort();
270 : #endif
271 : }
272 :
273 : void
274 334 : ImageSampler::vtkShiftAndScale()
275 : {
276 : #ifdef LIBMESH_HAVE_VTK
277 : // Capture the parameters
278 334 : double shift = _is_pars.get<double>("shift");
279 334 : double scale = _is_pars.get<double>("scale");
280 :
281 : // Do nothing if shift and scale are not set
282 334 : if (shift == 0 && scale == 1)
283 295 : return;
284 :
285 : // Perform the scaling and offset actions
286 39 : _shift_scale_filter = vtkSmartPointer<vtkImageShiftScale>::New();
287 39 : _shift_scale_filter->SetOutputScalarTypeToDouble();
288 :
289 39 : _shift_scale_filter->SetInputConnection(_algorithm);
290 39 : _shift_scale_filter->SetShift(shift);
291 39 : _shift_scale_filter->SetScale(scale);
292 39 : _shift_scale_filter->Update();
293 :
294 : // Update the pointers
295 39 : _data = _shift_scale_filter->GetOutput();
296 39 : _algorithm = _shift_scale_filter->GetOutputPort();
297 : #endif
298 : }
299 :
300 : void
301 334 : ImageSampler::vtkThreshold()
302 : {
303 : #ifdef LIBMESH_HAVE_VTK
304 : // Do nothing if threshold not set
305 334 : if (!_is_pars.isParamValid("threshold"))
306 260 : return;
307 :
308 : // Error if both upper and lower are not set
309 74 : if (!_is_pars.isParamValid("upper_value") || !_is_pars.isParamValid("lower_value"))
310 0 : mooseError("When thresholding is applied, both the upper_value and lower_value parameters must "
311 : "be set");
312 :
313 : // Create the thresholding object
314 74 : _image_threshold = vtkSmartPointer<vtkImageThreshold>::New();
315 :
316 : // Set the data source
317 74 : _image_threshold->SetInputConnection(_algorithm);
318 :
319 : // Setup the thresholding options
320 74 : _image_threshold->ThresholdByUpper(_is_pars.get<Real>("threshold"));
321 74 : _image_threshold->ReplaceInOn();
322 74 : _image_threshold->SetInValue(_is_pars.get<Real>("upper_value"));
323 74 : _image_threshold->ReplaceOutOn();
324 74 : _image_threshold->SetOutValue(_is_pars.get<Real>("lower_value"));
325 74 : _image_threshold->SetOutputScalarTypeToDouble();
326 :
327 : // Perform the thresholding
328 74 : _image_threshold->Update();
329 :
330 : // Update the pointers
331 74 : _data = _image_threshold->GetOutput();
332 74 : _algorithm = _image_threshold->GetOutputPort();
333 : #endif
334 : }
|
