Cogs.Core
GeoTiff.cpp
1#include "GeoTiff.h"
2
3#include "Foundation/Logging/Logger.h"
4
5#include <tiffio.hxx>
6
7#include <cmath>
8
9namespace {
10 using namespace Cogs::Core::TerrainProvider;
11 Cogs::Logging::Log logger = Cogs::Logging::getLogger("GeoTiff");
12
13 // LIBTIFF IO callbacks
14
15
16 static tmsize_t readProc(thandle_t clientData, void* ptr, tmsize_t bytes)
17 {
18 auto* self = (GeoTiffView*)clientData;
19 auto byteCount = std::min(size_t(bytes), self->size - self->offset);
20 std::memcpy(ptr, (uint8_t*)(self->ptr) + self->offset, byteCount);
21 self->offset += byteCount;
22 return tmsize_t(byteCount);
23 }
24
25 static tmsize_t writeProc(thandle_t, void*, tmsize_t)
26 {
27 assert(false && "Write should not be invoked");
28 return 0;
29 }
30
31 static toff_t seekProc(thandle_t clientData, toff_t off, int whence) {
32 auto* self = (GeoTiffView*)clientData;
33 switch (whence) {
34 case SEEK_SET:
35 self->offset = std::min<toff_t>(self->size, off);
36 break;
37 case SEEK_CUR:
38 self->offset = std::min<toff_t>(self->size, self->offset + off);
39 break;
40 case SEEK_END:
41 self->offset = self->size;
42 break;
43 default:
44 assert(false && "Invalid whence");
45 break;
46 }
47 return 0;
48 }
49
50 static toff_t sizeProc(thandle_t clientData) {
51 auto* self = (GeoTiffView*)clientData;
52 return self->size;
53 }
54
55 static int closeProc(thandle_t) { return 0; }
56
57 static int mapFileProc(thandle_t clientData, void** base, toff_t* size)
58 {
59 auto* self = (GeoTiffView*)clientData;
60 *base = (void*)self->ptr;
61 *size = self->size;
62 return 1;
63 }
64
65 static void unmapFileProc(thandle_t /*clientData*/, void* /*base*/, toff_t /*size*/) { }
66
67 struct IOWriteWrapper
68 {
69 Cogs::Memory::MemoryBuffer& out;
70 size_t offset = 0;
71 size_t size = 0;
72
73 void accommodate(size_t high)
74 {
75 while (out.size() < high) {
76 out.resize(out.size() + 1024 * 1024);
77 }
78 }
79
80 static tmsize_t writeProc(thandle_t clientData, void* ptr, tmsize_t bytes)
81 {
82 auto* self = reinterpret_cast<IOWriteWrapper*>(clientData);
83 self->accommodate(self->offset + bytes);
84 std::memcpy(static_cast<uint8_t*>(self->out.data()) + self->offset, ptr, bytes);
85 self->offset += bytes;
86 if (self->size < self->offset) {
87 self->size = self->offset;
88 }
89 return bytes;
90 }
91
92 static toff_t sizeProc(thandle_t clientData) {
93 auto* self = reinterpret_cast<IOWriteWrapper*>(clientData);
94 return self->size;
95 }
96
97 static int closeProc(thandle_t clientData) {
98 auto* self = reinterpret_cast<IOWriteWrapper*>(clientData);
99 self->out.resize(self->size);
100 return 0;
101 }
102
103 static tmsize_t readProc(thandle_t /*lientData*/, void* /*ptr*/, tmsize_t /*bytes*/)
104 {
105 assert(false);
106 return 0;
107 }
108
109 static toff_t seekProc(thandle_t clientData, toff_t off, int whence) {
110 auto* self = reinterpret_cast<IOWriteWrapper*>(clientData);
111
112 switch (whence) {
113 case SEEK_SET:
114 //self->offset = std::min(self->size, off);
115 self->offset = off;
116 break;
117 case SEEK_CUR:
118 //self->offset = std::min(self->size, self->offset + off);
119 self->offset = self->offset + off;
120 break;
121 case SEEK_END:
122 self->offset = self->size;
123 break;
124 default:
125 assert(false && "Invalid whence");
126 break;
127 }
128 if (self->size < self->offset) {
129 self->size = self->offset;
130 self->accommodate(self->size);
131 }
132
133 return self->offset;
134 }
135 };
136
137 // Registering of custom tags (=GEOTIFF tags)
138
139#define TIFFTAG_GDAL_METADATA 42112 // XML fragment
140#define TIFFTAG_GDAL_NODATA 42113 // text string with number
141#define TIFFTAG_MODELPIXELSCALE 33550 // 3 doubles
142#define TIFFTAG_MODELTIEPOINT 33922 // 6 double's per tie point
143#define TIFFTAG_MODELTRANSFORMATION 34264 // 16 double's row-major order matrix
144#define TIFFTAG_GEOKEYDIRECTORY 34735
145
146#define TIFFTAG_INTERGRAPH_MATRIX 33920
147#define TIFFTAG_JPL_CARTO_IFD 34263
148#define TIFFTAG_GEODOUBLEPARAMS 34736
149#define TIFFTAG_GEOASCIIPARAMS 34737
150
151 static const TIFFFieldInfo extraFieldInfo[] = {
152 { TIFFTAG_GDAL_METADATA, -1, -1, TIFF_ASCII, FIELD_CUSTOM, 1, 0, const_cast<char*>("GDALMetaData") },
153 { TIFFTAG_GDAL_NODATA, -1, -1, TIFF_ASCII, FIELD_CUSTOM, 1, 0, const_cast<char*>("GDALNoData") },
154 { TIFFTAG_MODELTIEPOINT, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, 1, 1, const_cast<char*>("ModelTiepointTag") },
155 { TIFFTAG_MODELPIXELSCALE, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, 1, 1, const_cast<char*>("ModelPixelScaleTag") },
156 { TIFFTAG_MODELTRANSFORMATION, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, 1, 1, const_cast<char*>("ModelTransformation") },
157 { TIFFTAG_GEOKEYDIRECTORY, -1, -1, TIFF_SHORT, FIELD_CUSTOM, 1, 1, const_cast<char*>("GeoKeyDirectory") },
158 { TIFFTAG_INTERGRAPH_MATRIX, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, 1, 1, const_cast<char*>("Intergraph TransformationMatrix") },
159 { TIFFTAG_GEODOUBLEPARAMS, -1, -1, TIFF_DOUBLE, FIELD_CUSTOM, 1, 1, const_cast<char*>("GeoDoubleParams") },
160 { TIFFTAG_GEOASCIIPARAMS, -1, -1, TIFF_ASCII, FIELD_CUSTOM, 1, 0, const_cast<char*>("GeoASCIIParams") },
161 { TIFFTAG_JPL_CARTO_IFD, 1, 1, TIFF_LONG, FIELD_CUSTOM, 1, 1, const_cast<char*>("JPL Carto IFD offset") },
162 };
163
164 static TIFFExtendProc previousExtender;
165
166 static void defaultDirectory(TIFF* tif)
167 {
168 TIFFMergeFieldInfo(tif, extraFieldInfo, sizeof(extraFieldInfo) / sizeof(extraFieldInfo[0]));
169 if (previousExtender) {
170 previousExtender(tif);
171 }
172 }
173
174 static struct CustomTIFFFIelds
175 {
176 CustomTIFFFIelds()
177 {
178 LOG_DEBUG(logger, "Registering custom TIFF fields");
179 previousExtender = TIFFSetTagExtender(defaultDirectory);
180 }
181 } customTIFFFields;
182
183}
184
185Cogs::Core::TerrainProvider::GeoTiffView::GeoTiffView(const void* ptr, size_t size)
186 : ptr(ptr), size(size), offset(0)
187{
188 tif = TIFFClientOpen("contents", "r", this,
189 readProc, writeProc,
190 seekProc, closeProc, sizeProc,
191 mapFileProc, unmapFileProc);
192 if (tif == nullptr) {
193 LOG_ERROR(logger, "Failed to open TIFF");
194 return;
195 }
196
197 // Inspect tags
198 uint16_t* tmp_u16p = nullptr;
199 uint32_t tmp_u32 = 0;
200 uint16_t tmp_u16 = 0;
201 double tmp_f64 = 0.0;
202 double* tmp_f64p = nullptr;
203 const char* tmp_str = nullptr;
204
205 // Image size
206 if (!TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &tmp_u32)) { LOG_ERROR(logger, "missing TIFFTAG_IMAGEWIDTH"); return; } else { width = tmp_u32; }
207 if (!TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &tmp_u32)) { LOG_ERROR(logger, "missing TIFFTAG_IMAGELENGTH"); return; } else { height = tmp_u32; }
208
209 // Pixel format
210 uint16_t bps = 0;
211 uint16_t ssp = 0;
212 uint16_t fmt = 0; // SAMPLEFORMAT_UINT, SAMPLEFORMAT_INT, SAMPLEFORMAT_IEEEFP, SAMPLEFORMAT_VOID, SAMPLEFORMAT_COMPLEXINT, SAMPLEFORMAT_COMPLEXIEEEFP
213 if (!TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &bps)) { LOG_ERROR(logger, "missing TIFFTAG_BITSPERSAMPLE"); return; }
214 if (!TIFFGetField(tif, TIFFTAG_SAMPLESPERPIXEL, &ssp)) { LOG_ERROR(logger, "missing TIFFTAG_SAMPLESPERPIXEL"); return; }
215 if (!TIFFGetField(tif, TIFFTAG_SAMPLEFORMAT, &fmt)) { LOG_ERROR(logger, "missing TIFFTAG_SAMPLEFORMAT"); return; }
216 if ((bps == 32) && (ssp == 1) && (fmt == SAMPLEFORMAT_IEEEFP)) {
217 format = Cogs::TextureFormat::R32_FLOAT;
218 }
219 else {
220 LOG_ERROR(logger, "Unsupported pixel format bits_per_sample=%u, samples_per_pixel=%u, sample_format=%u", unsigned(bps), unsigned(ssp), unsigned(fmt));
221 return;
222 }
223
224 // min/max sample value
225 if (TIFFGetField(tif, TIFFTAG_SMINSAMPLEVALUE, &tmp_f64)) { minSample = float(tmp_f64); }
226 else if (TIFFGetField(tif, TIFFTAG_MINSAMPLEVALUE, &tmp_u16)) { minSample = tmp_u16; }
227 if (TIFFGetField(tif, TIFFTAG_SMAXSAMPLEVALUE, &tmp_f64)) { maxSample = float(tmp_f64); }
228 else if (TIFFGetField(tif, TIFFTAG_MAXSAMPLEVALUE, &tmp_u16)) { maxSample = tmp_u16; }
229
230 // No-data is encoded as an ascii string
231 if (TIFFGetField(tif, TIFFTAG_GDAL_NODATA, &tmp_str)) {
232 if (strcmp("nan", tmp_str) == 0) {
233 noData = std::numeric_limits<float>::quiet_NaN();
234 }
235 else {
236 char* end = nullptr;
237 noData = std::strtof(tmp_str, &end);
238 if (!end || *end != '\0') {
239 LOG_WARNING(logger, "Failed to parse nodata value '%s'", tmp_str);
240 }
241 }
242 }
243
244 if (TIFFGetField(tif, TIFFTAG_MODELPIXELSCALE, &tmp_u16, &tmp_f64p)) {
245 if (tmp_u16 == 3) {
246 scale[0] = tmp_f64p[0];
247 scale[1] = tmp_f64p[1];
248 scale[2] = tmp_f64p[2];
249 }
250 else {
251 LOG_WARNING(logger, "TIFFTAG_MODELPIXELSACE has unexpected count %u, expected 3.", unsigned(tmp_u16));
252 }
253 }
254
255 if (TIFFGetField(tif, TIFFTAG_MODELTIEPOINT, &tmp_u16, &tmp_f64p)) {
256 if (tmp_u16 == 6) {
257 for (size_t i = 0; i < 6; i++) { tiePoint[i] = tmp_f64p[i]; }
258 origin[0] = tiePoint[3] - (scale[0] == 0.0 ? 1.0 : scale[0]) * tiePoint[0];
259 origin[1] = tiePoint[4] - (scale[1] == 0.0 ? 1.0 : scale[1]) * tiePoint[1];
260 origin[2] = tiePoint[5] - (scale[2] == 0.0 ? 1.0 : scale[2]) * tiePoint[2];
261 }
262 else {
263 LOG_WARNING(logger, "TIFFTAG_MODELTIEPOINT has unexpected count %u, expected 6", unsigned(tmp_u16));
264 }
265 }
266
267 if (TIFFGetField(tif, TIFFTAG_MODELTRANSFORMATION, &tmp_u16, &tmp_f64p)) {
268 if (tmp_u16 == 16) {
269 scale[0] = tmp_f64p[0];
270 scale[1] = tmp_f64p[5];
271 scale[2] = tmp_f64p[10];
272 origin[0] = tmp_f64p[3];
273 origin[1] = tmp_f64p[7];
274 origin[2] = tmp_f64p[11];
275 }
276 else {
277 LOG_WARNING(logger, "TIFFTAG_MODELTIEPOINT has unexpected count %u, expected 16", unsigned(tmp_u16));
278 }
279 }
280
281 if (TIFFGetField(tif, TIFFTAG_GDAL_METADATA, &tmp_str)) {
282 LOG_DEBUG(logger, "Unprocessed GDAL metadata: %s", tmp_str);
283 }
284
285 if (TIFFGetField(tif, TIFFTAG_GEOKEYDIRECTORY, &tmp_u16, &tmp_u16p)) {
286 LOG_DEBUG(logger, "Unprocessed %u geokeys: %u, %u, %u, %u,...", unsigned(tmp_u16),
287 unsigned(tmp_u16p[0]), unsigned(tmp_u16p[1]),
288 unsigned(tmp_u16p[2]), unsigned(tmp_u16p[3]));
289 }
290
291#if 0
292 LOG_DEBUG(logger, "width=%u, height=%u, noData=%e, min=%f, max=%f, scale=[%f, %f, %f], origin=[%f, %f, %f]",
293 width, height,
294 noData.value_or(std::numeric_limits<float>::infinity()),
295 minSample.value_or(std::numeric_limits<float>::infinity()),
296 maxSample.value_or(std::numeric_limits<float>::infinity()),
297 scale[0], scale[1], scale[2],
298 origin[0], origin[1], origin[2]);
299#endif
300 valid = true;
301}
302
303Cogs::Core::TerrainProvider::GeoTiffView::~GeoTiffView()
304{
305 if (tif) { TIFFClose(tif); }
306}
307
308bool Cogs::Core::TerrainProvider::GeoTiffView::getFullImage(Memory::MemoryBuffer& output, float suggested_nodata)
309{
310 if (!valid) return false;
311
312 switch (format) {
313
314 case Cogs::TextureFormat::R32_FLOAT: {
315 if (TIFFIsTiled(tif)) {
316 uint32_t tileWidth;
317 uint32_t tileHeight;
318
319 TIFFGetField(tif, TIFFTAG_TILEWIDTH, &tileWidth);
320 TIFFGetField(tif, TIFFTAG_TILELENGTH, &tileHeight);
321
322 auto buf = _TIFFmalloc(TIFFTileSize(tif));
323 if (buf == nullptr) return false;
324
325 size_t widthMultipleOfTile = static_cast<size_t>(std::ceil((float)width / tileWidth) * tileWidth);
326 size_t heightMultipleOfTile = static_cast<size_t>(std::ceil((float)height / tileHeight) * tileHeight);
327 output.resize(sizeof(float) * widthMultipleOfTile * heightMultipleOfTile);
328
329 float* rowStart = reinterpret_cast<float*>(output.data());
330
331 for (uint32_t y = 0; y < height; y += tileHeight, rowStart += (tileHeight * width)) {
332 float* tileStart = rowStart;
333
334 for (uint32_t x = 0; x < width; x += tileWidth, tileStart += tileWidth) {
335 if (TIFFReadTile(tif, buf, x, y, 0, 0) > 0) {
336 float* write = tileStart;
337 float* read = reinterpret_cast<float*>(buf);
338
339 for (uint32_t dy = tileHeight; dy--; write += width, read += tileWidth) {
340 memcpy(write, read, tileWidth * sizeof(float));
341 }
342 }
343 }
344 }
345 _TIFFfree(buf);
346 }
347 else {
348 auto* buf = _TIFFmalloc(TIFFStripSize(tif));
349
350 if (buf == nullptr) return false;
351
352 size_t o = 0;
353 output.resize(sizeof(float) * width * height);
354 for (tstrip_t strip = 0, noOfStrips = TIFFNumberOfStrips(tif); strip < noOfStrips; strip++) {
355 auto byteCount = (size_t)TIFFReadEncodedStrip(tif, strip, buf, (tsize_t)-1);
356 //size_t n = std::min(output.size() - o, byteCount);
357 std::memcpy((uint8_t*)output.data() + o, buf, byteCount);
358 o += byteCount;
359 assert(o <= output.size());
360 }
361 _TIFFfree(buf);
362 }
363
364 auto* oo = (float*)output.data();
365 auto* ee = oo + size_t(width) * size_t(height);
366 float match_nodata = noData.value_or(suggested_nodata);
367 float replace_nodata = std::numeric_limits<float>::quiet_NaN();
368 float min = std::numeric_limits<float>::max();
369 float max = -std::numeric_limits<float>::max();
370 if (std::isfinite(match_nodata)) { // regular value
371 for (auto* p = oo; p < ee; p++) {
372 float v = *p;
373 if (v == noData) {
374 *p = replace_nodata;
375 }
376 else {
377 if (min > v) min = v;
378 if (max < v) max = v;
379 }
380 }
381 }
382 else {
383 for (auto * p = oo; p < ee; p++) {
384 float v = *p;
385 if (!std::isfinite(v)) {
386 *p = replace_nodata;
387 }
388 else {
389 if (min > v) min = v;
390 if (max < v) max = v;
391 }
392 }
393 }
394 if(min < max) {
395 minSample = min;
396 maxSample = max;
397 }
398 break;
399 }
400 default:
401 return false;
402 }
403 return true;
404}
405
406
407bool Cogs::Core::TerrainProvider::createTiff(Memory::MemoryBuffer& out, const Memory::MemoryBuffer& in, unsigned w, unsigned h, Cogs::TextureFormat format)
408{
409 IOWriteWrapper ww{ out };
410
411 uint16_t spp, bpp, fmt;
412 uint16_t compression = COMPRESSION_DEFLATE;
413 uint16_t predictor;
414 uint16_t photometric;
415
416 switch (format) {
417 case Cogs::Format::R32_FLOAT:
418 spp = 1;
419 bpp = 32;
420 fmt = SAMPLEFORMAT_IEEEFP;
421 predictor = PREDICTOR_FLOATINGPOINT;
422 photometric = PHOTOMETRIC_MINISBLACK;
423 break;
424 case Cogs::Format::R8G8B8A8_UINT:
425 spp = 4;
426 bpp = 8;
427 fmt = SAMPLEFORMAT_UINT;
428 predictor = PREDICTOR_HORIZONTAL;
429 photometric = PHOTOMETRIC_RGB;
430 break;
431
432 default:
433 LOG_ERROR(logger, "Unsupported format %u", (unsigned)format);
434 return false;
435 }
436
437 auto tif = TIFFClientOpen("contents", "w", &ww,
438 IOWriteWrapper::readProc, IOWriteWrapper::writeProc,
439 IOWriteWrapper::seekProc, IOWriteWrapper::closeProc, IOWriteWrapper::sizeProc,
440 nullptr, nullptr);
441 if (tif == nullptr) {
442 LOG_ERROR(logger, "Failed to open TIFF");
443 return false;
444 }
445
446 TIFFSetField(tif, TIFFTAG_IMAGEWIDTH, (uint32_t)w);
447 TIFFSetField(tif, TIFFTAG_IMAGELENGTH, (uint32_t)h);
448 TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL, spp);
449 TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, bpp);
450 TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT, fmt);
451 TIFFSetField(tif, TIFFTAG_COMPRESSION, compression);
452 TIFFSetField(tif, TIFFTAG_PREDICTOR, predictor);
453 TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, photometric);
454 TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
455 TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_TOPLEFT);
456 TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, TIFFDefaultStripSize(tif, w * spp));
457
458 auto* ptr = static_cast<const uint8_t*>(in.data());
459 for (unsigned j = 0; j < h; j++) {
460 if (TIFFWriteScanline(tif, (void*)ptr, j, 0) != 1) {
461 LOG_ERROR(logger, "Error writing scanline");
462 return false;
463 }
464 ptr += w * ((bpp * spp) / 8);
465 }
466 TIFFClose(tif);
467 return true;
468}
Cogs::Logging::Log
Log implementation class.
Definition: LogManager.h:139
Cogs::Logging::getLogger
constexpr Log getLogger(const char(&name)[LEN]) noexcept
Definition: LogManager.h:180
Generated by doxygen 1.9.6