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