OGS
VtkRaster.cpp
Go to the documentation of this file.
1
15#include "VtkRaster.h"
16
17#include <vtkBMPReader.h>
18#include <vtkImageData.h>
19#include <vtkImageImport.h>
20#include <vtkImageReader2.h>
21#include <vtkJPEGReader.h>
22#include <vtkPNGReader.h>
23#include <vtkTIFFReader.h>
24
25#include <QFileInfo>
26#include <algorithm>
27#include <cmath>
28#include <limits>
29
30#ifdef GEOTIFF_FOUND
31#include <geo_tiffp.h>
32#include <xtiffio.h>
33#endif
34
35#include <memory>
36
38#include "BaseLib/FileTools.h"
39#include "BaseLib/Logging.h"
40#include "BaseLib/StringTools.h"
41#include "GeoLib/Raster.h"
42
43vtkImageAlgorithm* VtkRaster::loadImage(const std::string& fileName)
44{
45 QFileInfo fileInfo(QString::fromStdString(fileName));
46
47 std::unique_ptr<GeoLib::Raster> raster(nullptr);
48 if (fileInfo.suffix().toLower() == "asc")
49 {
50 raster.reset(
52 }
53 else if (fileInfo.suffix().toLower() == "grd")
54 {
55 raster.reset(
57 }
58 if (raster)
59 {
60 return VtkRaster::loadImageFromArray(raster->begin(),
61 raster->getHeader());
62 }
63 if ((fileInfo.suffix().toLower() == "tif") ||
64 (fileInfo.suffix().toLower() == "tiff"))
65 {
66#ifdef GEOTIFF_FOUND
67 return loadImageFromTIFF(fileName);
68#else
69 ERR("VtkRaster::loadImage(): GeoTiff file format not supported in this "
70 "version! Trying to parse as Tiff-file.");
71 return loadImageFromFile(fileName);
72#endif
73 }
74
75 return loadImageFromFile(fileName);
76}
77vtkImageImport* VtkRaster::loadImageFromArray(double const* const data_array,
79{
80 const unsigned length = header.n_rows * header.n_cols * header.n_depth;
81 auto* data = new float[length * 2];
82 float max_val =
83 static_cast<float>(*std::max_element(data_array, data_array + length));
84 for (unsigned j = 0; j < length; ++j)
85 {
86 data[j * 2] = static_cast<float>(data_array[j]);
87 if (fabs(data[j * 2] - header.no_data) <
88 std::numeric_limits<double>::epsilon())
89 {
90 data[j * 2] = max_val;
91 data[j * 2 + 1] = 0;
92 }
93 else
94 {
95 data[j * 2 + 1] = max_val;
96 }
97 }
98
99 vtkImageImport* image = vtkImageImport::New();
100 image->SetDataSpacing(header.cell_size, header.cell_size, header.cell_size);
101 image->SetDataOrigin(header.origin[0] + (header.cell_size / 2.0),
102 header.origin[1] + (header.cell_size / 2.0),
103 0); // translate whole mesh by half a pixel in x and y
104 image->SetWholeExtent(0, header.n_cols - 1, 0, header.n_rows - 1, 0,
105 header.n_depth - 1);
106 image->SetDataExtent(0, header.n_cols - 1, 0, header.n_rows - 1, 0,
107 header.n_depth - 1);
108 image->SetDataExtentToWholeExtent();
109 image->SetDataScalarTypeToFloat();
110 image->SetNumberOfScalarComponents(2);
111 image->SetImportVoidPointer(data, 0);
112 image->Update();
113
114 return image;
115}
116
117#ifdef GEOTIFF_FOUND
118vtkImageAlgorithm* VtkRaster::loadImageFromTIFF(const std::string& fileName)
119{
120 TIFF* tiff = XTIFFOpen(fileName.c_str(), "r");
121
122 if (tiff)
123 {
124 GTIF* geoTiff = GTIFNew(tiff);
125
126 int version[3];
127 int count(0);
128 GTIFDirectoryInfo(geoTiff, version, &count);
129
130 if (count == 0)
131 WARN("VtkRaster::loadImageFromTIFF - file is not georeferenced.");
132
133 if (geoTiff)
134 {
135 double x0 = 0.0;
136 double y0 = 0.0;
137 double cellsize = 1.0;
138 int imgWidth = 0;
139 int imgHeight = 0;
140 int nImages = 0;
141 int pntCount = 0;
142 double* pnts = nullptr;
143
144 // get actual number of images in the tiff file
145 do
146 {
147 ++nImages;
148 } while (TIFFReadDirectory(tiff));
149 if (nImages > 1)
150 INFO(
151 "VtkRaster::loadImageFromTIFF() - File contains {:d} "
152 "images. This method is not tested for this case.",
153 nImages);
154
155 // get image size
156 TIFFGetField(tiff, TIFFTAG_IMAGEWIDTH, &imgWidth);
157 TIFFGetField(tiff, TIFFTAG_IMAGELENGTH, &imgHeight);
158
159 // get cellsize
160 // Note: GeoTiff allows anisotropic pixels. This is not supported
161 // here and equilateral pixels are assumed.
162 if (TIFFGetField(tiff, GTIFF_PIXELSCALE, &pntCount, &pnts))
163 {
164 if (pnts[0] != pnts[1])
165 WARN(
166 "VtkRaster::loadImageFromTIFF(): Original raster data "
167 "has anisotrop pixel size!");
168 cellsize = pnts[0];
169 }
170
171 // get upper left point / origin
172 if (TIFFGetField(tiff, GTIFF_TIEPOINTS, &pntCount, &pnts))
173 {
174 x0 = pnts[3];
175 y0 = pnts[4] -
176 (imgHeight * cellsize); // the origin should be the lower
177 // left corner of the img
178 }
179
180 // read pixel values
181 auto* pixVal = static_cast<uint32*>(
182 _TIFFmalloc(imgWidth * imgHeight * sizeof(uint32)));
183 if ((imgWidth > 0) && (imgHeight > 0))
184 {
185 if (!TIFFReadRGBAImage(tiff, imgWidth, imgHeight, pixVal, 0))
186 {
187 ERR("VtkRaster::loadImageFromTIFF(): reading GeoTIFF "
188 "file.");
189 _TIFFfree(pixVal);
190 GTIFFree(geoTiff);
191 XTIFFClose(tiff);
192 return nullptr;
193 }
194 }
195
196 // check for colormap
197 uint16 photometric;
198 TIFFGetField(tiff, TIFFTAG_PHOTOMETRIC, &photometric);
199 // read colormap
200 uint16 *cmap_red = nullptr, *cmap_green = nullptr,
201 *cmap_blue = nullptr;
202 int colormap_used = TIFFGetField(tiff, TIFFTAG_COLORMAP, &cmap_red,
203 &cmap_green, &cmap_blue);
204
205 auto* data = new float[imgWidth * imgHeight * 4];
206 auto* pxl(new int[4]);
207 for (int j = 0; j < imgHeight; ++j)
208 {
209 int lineindex = j * imgWidth;
210 for (int i = 0; i < imgWidth; ++i)
211 { // scale intensities and set nodata values to white (i.e. the
212 // background colour)
213 unsigned pxl_idx(lineindex + i);
214 unsigned pos = 4 * (pxl_idx);
215 if (photometric == 1 && colormap_used == 1)
216 {
217 int idx = TIFFGetR(pixVal[pxl_idx]);
218 data[pos] = static_cast<float>(cmap_red[idx] >> 8);
219 data[pos + 1] =
220 static_cast<float>(cmap_green[idx] >> 8);
221 data[pos + 2] = static_cast<float>(cmap_blue[idx] >> 8);
222 data[pos + 3] = 1;
223 }
224 else
225 {
226 data[pos] =
227 static_cast<float>(TIFFGetR(pixVal[pxl_idx]));
228 data[pos + 1] =
229 static_cast<float>(TIFFGetG(pixVal[pxl_idx]));
230 data[pos + 2] =
231 static_cast<float>(TIFFGetB(pixVal[pxl_idx]));
232 data[pos + 3] =
233 static_cast<float>(TIFFGetA(pixVal[pxl_idx]));
234 }
235 }
236 }
237 delete[] pxl;
238
239 // set transparency values according to maximum pixel value
240 if (photometric == 1)
241 {
242 float max_val(0);
243 unsigned nPixels = 4 * imgWidth * imgHeight;
244 for (unsigned j = 0; j < nPixels; ++j)
245 {
246 if (data[j] > max_val)
247 {
248 max_val = data[j];
249 }
250 }
251
252 for (unsigned j = 0; j < nPixels; j += 4)
253 {
254 data[j + 3] = max_val;
255 }
256 }
257
258 vtkImageImport* image = vtkImageImport::New();
259 image->SetDataOrigin(x0, y0, 0);
260 image->SetDataSpacing(cellsize, cellsize, cellsize);
261 image->SetWholeExtent(0, imgWidth - 1, 0, imgHeight - 1, 0, 0);
262 image->SetDataExtent(0, imgWidth - 1, 0, imgHeight - 1, 0, 0);
263 image->SetDataExtentToWholeExtent();
264 image->SetDataScalarTypeToFloat();
265 image->SetNumberOfScalarComponents(4);
266 image->SetImportVoidPointer(data, 0);
267 image->Update();
268
269 _TIFFfree(pixVal);
270 GTIFFree(geoTiff);
271 XTIFFClose(tiff);
272 return image;
273 }
274
275 XTIFFClose(tiff);
276 ERR("VtkRaster::loadImageFromTIFF() - File not recognised as "
277 "GeoTIFF-Image.");
278 return nullptr;
279 }
280
281 ERR("VtkRaster::loadImageFromTIFF() - File not recognised as TIFF-Image.");
282 return nullptr;
283}
284#endif
285
286vtkImageReader2* VtkRaster::loadImageFromFile(const std::string& fileName)
287{
288 QString file_name(QString::fromStdString(fileName));
289 QFileInfo fi(file_name);
290 vtkImageReader2* image(nullptr);
291
292 if (fi.suffix().toLower() == "png")
293 {
294 image = vtkPNGReader::New();
295 }
296 else if ((fi.suffix().toLower() == "tif") ||
297 (fi.suffix().toLower() == "tiff"))
298 {
299 image = vtkTIFFReader::New();
300 }
301 else if ((fi.suffix().toLower() == "jpg") ||
302 (fi.suffix().toLower() == "jpeg"))
303 {
304 image = vtkJPEGReader::New();
305 }
306 else if (fi.suffix().toLower() == "bmp")
307 {
308 image = vtkBMPReader::New();
309 }
310 else
311 {
312 ERR("VtkRaster::readImageFromFile(): File format not supported, please "
313 "convert to BMP, JPG, PNG or TIFF.");
314 return nullptr;
315 }
316
317 image->SetFileName(fileName.c_str());
318 image->GetOutput()->AllocateScalars(VTK_FLOAT, 1);
319 image->Update();
320 readWorldFile(fileName, image);
321 return image;
322}
323
324std::string VtkRaster::findWorldFile(std::string const& filename)
325{
326 std::string const no_ext = BaseLib::dropFileExtension(filename);
327
328 constexpr std::array supported_extensions = {
329 ".pgw", ".pngw", ".pgwx", ".jgw", ".jpgw", ".jgwx", ".tfw",
330 ".tifw", ".tfwx", ".bpw", ".bmpw", ".bpwx", ".wld"};
331
332 auto const res =
333 std::find_if(supported_extensions.begin(), supported_extensions.end(),
334 [&no_ext](auto const& ext) -> bool
335 { return BaseLib::IsFileExisting(no_ext + ext); });
336 if (res != supported_extensions.end())
337 {
338 return no_ext + *res;
339 }
340
341 // no world file found
342 return {};
343}
344
345bool VtkRaster::readWorldFile(std::string const& filename,
346 vtkImageReader2* image)
347{
348 std::string const world_file = findWorldFile(filename);
349 if (world_file.empty())
350 {
351 WARN("No world file found. Image is not georeferenced.");
352 return false;
353 }
354
355 std::ifstream in(world_file.c_str());
356 if (!in.is_open())
357 {
358 ERR("VtkRaster::readWorldFile(): Could not open file {:s}.", filename);
359 return false;
360 }
361
362 std::string line;
363 // x-scaling
364 if (!std::getline(in, line))
365 {
366 return false;
367 }
368 double const delta_x = BaseLib::str2number<double>(line);
369 // 2x rotation (ignored)
370 if (!(std::getline(in, line) && std::getline(in, line)))
371 {
372 return false;
373 }
374 // negative y-scaling
375 if (!std::getline(in, line))
376 {
377 return false;
378 }
379 double const delta_y = BaseLib::str2number<double>(line);
380 if (delta_x != -delta_y)
381 WARN(
382 "Anisotropic pixel size detected ({:f} vs {:f}). An isotropic "
383 "spacing of {:f} is assumed, be aware results may be wrong.",
384 delta_x, delta_y, delta_x);
385 // x-translation
386 if (!std::getline(in, line))
387 {
388 return false;
389 }
390 double const x0 = BaseLib::str2number<double>(line);
391 // y-translation
392 if (!std::getline(in, line))
393 {
394 return false;
395 }
396 double const y0 = BaseLib::str2number<double>(line);
397
398 int extent[3];
399 image->GetOutput()->GetDimensions(extent);
400 image->SetDataSpacing(delta_x, delta_x, delta_x);
401 // for GIS the origin is on the lower left, for VTK it is on the upper left
402 double const vtk_y0 = y0 + (extent[1] * delta_y);
403 image->SetDataOrigin(x0, vtk_y0, 0);
404 return true;
405}
Definition of the AsciiRasterInterface class.
Filename manipulation routines.
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:34
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:44
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:39
Definition of the GeoLib::Raster class.
Definition of string helper functions.
Definition of the VtkRaster class.
static GeoLib::Raster * getRasterFromSurferFile(std::string const &fname)
Reads a Surfer GRD raster file.
static GeoLib::Raster * getRasterFromASCFile(std::string const &fname)
Reads an ArcGis ASC raster file.
static std::string findWorldFile(const std::string &filename)
Definition: VtkRaster.cpp:324
static vtkImageAlgorithm * loadImageFromTIFF(const std::string &fileName)
Definition: VtkRaster.cpp:118
static vtkImageReader2 * loadImageFromFile(const std::string &fileName)
Definition: VtkRaster.cpp:286
static bool readWorldFile(std::string const &filename, vtkImageReader2 *image)
Definition: VtkRaster.cpp:345
static vtkImageImport * loadImageFromArray(double const *const data_array, GeoLib::RasterHeader header)
Returns a VtkImageAlgorithm from an array of pixel values and some image meta data.
Definition: VtkRaster.cpp:77
static vtkImageAlgorithm * loadImage(const std::string &fileName)
Loads an image- or raster-file into an vtkImageAlgorithm-Object.
Definition: VtkRaster.cpp:43
std::string dropFileExtension(std::string const &filename)
Definition: FileTools.cpp:169
Contains the relevant information when storing a geoscientific raster data.
Definition: Raster.h:27
std::size_t n_depth
Definition: Raster.h:30
double cell_size
Definition: Raster.h:32
MathLib::Point3d origin
Definition: Raster.h:31
std::size_t n_cols
Definition: Raster.h:28
std::size_t n_rows
Definition: Raster.h:29