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 14988 : ImageSampler::validParams()
21 : {
22 : // Define the general parameters
23 14988 : InputParameters params = emptyInputParameters();
24 14988 : params += FileRangeBuilder::validParams();
25 :
26 14988 : params.addParam<Point>("origin", "Origin of the image (defaults to mesh origin)");
27 14988 : params.addParam<Point>("dimensions",
28 : "x,y,z dimensions of the image (defaults to mesh dimensions)");
29 14988 : 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 14988 : params.addParam<double>("shift", 0, "Value to add to all pixels; occurs prior to scaling");
37 44964 : params.addParam<double>(
38 29976 : "scale", 1, "Multiplier to apply to all pixel values; occurs after shifting");
39 14988 : params.addParamNamesToGroup("shift scale", "Rescale");
40 :
41 : // Threshold parameters
42 14988 : params.addParam<double>("threshold", "The threshold value");
43 44964 : params.addParam<double>(
44 29976 : "upper_value", 1, "The value to set for data greater than the threshold value");
45 44964 : params.addParam<double>(
46 29976 : "lower_value", 0, "The value to set for data less than the threshold value");
47 14988 : params.addParamNamesToGroup("threshold upper_value lower_value", "Threshold");
48 :
49 : // Flip image
50 14988 : params.addParam<bool>("flip_x", false, "Flip the image along the x-axis");
51 14988 : params.addParam<bool>("flip_y", false, "Flip the image along the y-axis");
52 14988 : params.addParam<bool>("flip_z", false, "Flip the image along the z-axis");
53 14988 : params.addParamNamesToGroup("flip_x flip_y flip_z", "Flip");
54 :
55 14988 : return params;
56 0 : }
57 :
58 377 : ImageSampler::ImageSampler(const InputParameters & parameters)
59 : : FileRangeBuilder(parameters),
60 : #ifdef LIBMESH_HAVE_VTK
61 377 : _data(nullptr),
62 377 : _algorithm(nullptr),
63 : #endif
64 377 : _is_pars(parameters),
65 377 : _is_console(
66 : (parameters.getCheckedPointerParam<MooseApp *>("_moose_app"))->getOutputWarehouse()),
67 377 : _flip({{_is_pars.get<bool>("flip_x"),
68 377 : _is_pars.get<bool>("flip_y"),
69 754 : _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 377 : }
78 :
79 : void
80 374 : ImageSampler::setupImageSampler(MooseMesh & mesh)
81 : {
82 : // Don't warn that mesh or _is_pars are unused when VTK is not enabled.
83 374 : libmesh_ignore(mesh);
84 374 : libmesh_ignore(_is_pars);
85 :
86 : #ifdef LIBMESH_HAVE_VTK
87 : // Get access to the Mesh object
88 374 : BoundingBox bbox = MeshTools::create_bounding_box(mesh.getMesh());
89 :
90 : // Set the dimensions from the Mesh if not set by the User
91 374 : if (_is_pars.isParamValid("dimensions"))
92 112 : _physical_dims = _is_pars.get<Point>("dimensions");
93 :
94 : else
95 : {
96 262 : _physical_dims(0) = bbox.max()(0) - bbox.min()(0);
97 : #if LIBMESH_DIM > 1
98 262 : _physical_dims(1) = bbox.max()(1) - bbox.min()(1);
99 : #endif
100 : #if LIBMESH_DIM > 2
101 262 : _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 374 : if (_is_pars.isParamValid("origin"))
107 112 : _origin = _is_pars.get<Point>("origin");
108 : else
109 : {
110 262 : _origin(0) = bbox.min()(0);
111 : #if LIBMESH_DIM > 1
112 262 : _origin(1) = bbox.min()(1);
113 : #endif
114 : #if LIBMESH_DIM > 2
115 262 : _origin(2) = bbox.min()(2);
116 : #endif
117 : }
118 :
119 : // An array of filenames, to be filled in
120 374 : std::vector<std::string> filenames;
121 :
122 : // The file suffix, to be determined
123 374 : 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 374 : if (_status != 0)
133 : {
134 : // We don't have parameters, so see if we can get them from ImageMesh
135 28 : ImageMesh * image_mesh = dynamic_cast<ImageMesh *>(&mesh);
136 28 : 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 28 : filenames = image_mesh->filenames();
142 28 : 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 346 : filenames = this->filenames();
148 346 : file_suffix = fileSuffix();
149 : }
150 :
151 : // Storage for the file names
152 374 : _files = vtkSmartPointer<vtkStringArray>::New();
153 :
154 1342 : for (const auto & filename : filenames)
155 968 : _files->InsertNextValue(filename);
156 :
157 : // Error if no files where located
158 374 : 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 366 : if (file_suffix == "png")
165 366 : _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 366 : _is_console << "Reading image(s)..." << std::endl;
174 :
175 : // Extract the data
176 366 : _image->SetFileNames(_files);
177 366 : _image->Update();
178 366 : _data = _image->GetOutput();
179 366 : _algorithm = _image->GetOutputPort();
180 :
181 : // Set the image dimensions and voxel size member variable
182 366 : int * dims = _data->GetDimensions();
183 1464 : for (unsigned int i = 0; i < 3; ++i)
184 : {
185 1098 : _dims.push_back(dims[i]);
186 1098 : _voxel.push_back(_physical_dims(i) / _dims[i]);
187 : }
188 :
189 : // Set the dimensions of the image and bounding box
190 366 : _data->SetSpacing(_voxel[0], _voxel[1], _voxel[2]);
191 366 : _data->SetOrigin(_origin(0), _origin(1), _origin(2));
192 366 : _bounding_box.min() = _origin;
193 366 : _bounding_box.max() = _origin + _physical_dims;
194 :
195 : // Indicate data read is completed
196 366 : _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 366 : if (_is_pars.isParamValid("component"))
201 : {
202 18 : unsigned int n = _data->GetNumberOfScalarComponents();
203 18 : _component = _is_pars.get<unsigned int>("component");
204 18 : if (_component >= n)
205 4 : mooseError("'component' parameter must be empty or have a value of 0 to ", n - 1);
206 : }
207 : else
208 348 : _component = 0;
209 :
210 : // Apply filters, the toggling on and off of each filter is handled internally
211 362 : vtkMagnitude();
212 362 : vtkShiftAndScale();
213 362 : vtkThreshold();
214 : #endif
215 362 : }
216 :
217 : Real
218 2277060 : 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 2277060 : if (!_bounding_box.contains_point(p))
224 44604 : return 0.0;
225 :
226 : // Determine pixel coordinates
227 2232456 : std::vector<int> x(3, 0);
228 8929824 : for (int i = 0; i < LIBMESH_DIM; ++i)
229 : {
230 : // Compute position, only if voxel size is greater than zero
231 6697368 : if (_voxel[i] != 0)
232 : {
233 5789712 : 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 5789712 : if (x[i] == _dims[i])
237 124512 : x[i]--;
238 :
239 : // flip
240 5789712 : if (_flip[i])
241 87120 : x[i] = _dims[i] - x[i] - 1;
242 : }
243 : }
244 :
245 : // Return the image data at the given point
246 2232456 : 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 2232456 : }
253 :
254 : void
255 362 : ImageSampler::vtkMagnitude()
256 : {
257 : #ifdef LIBMESH_HAVE_VTK
258 : // Do nothing if 'component' is set
259 362 : if (_is_pars.isParamValid("component"))
260 14 : return;
261 :
262 : // Apply the greyscale filtering
263 348 : _magnitude_filter = vtkSmartPointer<vtkImageMagnitude>::New();
264 348 : _magnitude_filter->SetInputConnection(_algorithm);
265 348 : _magnitude_filter->Update();
266 :
267 : // Update the pointers
268 348 : _data = _magnitude_filter->GetOutput();
269 348 : _algorithm = _magnitude_filter->GetOutputPort();
270 : #endif
271 : }
272 :
273 : void
274 362 : ImageSampler::vtkShiftAndScale()
275 : {
276 : #ifdef LIBMESH_HAVE_VTK
277 : // Capture the parameters
278 362 : double shift = _is_pars.get<double>("shift");
279 362 : double scale = _is_pars.get<double>("scale");
280 :
281 : // Do nothing if shift and scale are not set
282 362 : if (shift == 0 && scale == 1)
283 320 : return;
284 :
285 : // Perform the scaling and offset actions
286 42 : _shift_scale_filter = vtkSmartPointer<vtkImageShiftScale>::New();
287 42 : _shift_scale_filter->SetOutputScalarTypeToDouble();
288 :
289 42 : _shift_scale_filter->SetInputConnection(_algorithm);
290 42 : _shift_scale_filter->SetShift(shift);
291 42 : _shift_scale_filter->SetScale(scale);
292 42 : _shift_scale_filter->Update();
293 :
294 : // Update the pointers
295 42 : _data = _shift_scale_filter->GetOutput();
296 42 : _algorithm = _shift_scale_filter->GetOutputPort();
297 : #endif
298 : }
299 :
300 : void
301 362 : ImageSampler::vtkThreshold()
302 : {
303 : #ifdef LIBMESH_HAVE_VTK
304 : // Do nothing if threshold not set
305 362 : if (!_is_pars.isParamValid("threshold"))
306 280 : return;
307 :
308 : // Error if both upper and lower are not set
309 82 : 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 82 : _image_threshold = vtkSmartPointer<vtkImageThreshold>::New();
315 :
316 : // Set the data source
317 82 : _image_threshold->SetInputConnection(_algorithm);
318 :
319 : // Setup the thresholding options
320 82 : _image_threshold->ThresholdByUpper(_is_pars.get<Real>("threshold"));
321 82 : _image_threshold->ReplaceInOn();
322 82 : _image_threshold->SetInValue(_is_pars.get<Real>("upper_value"));
323 82 : _image_threshold->ReplaceOutOn();
324 82 : _image_threshold->SetOutValue(_is_pars.get<Real>("lower_value"));
325 82 : _image_threshold->SetOutputScalarTypeToDouble();
326 :
327 : // Perform the thresholding
328 82 : _image_threshold->Update();
329 :
330 : // Update the pointers
331 82 : _data = _image_threshold->GetOutput();
332 82 : _algorithm = _image_threshold->GetOutputPort();
333 : #endif
334 : }
|
