Line data Source code
1 : /******************************************************************************
2 : * terragendataset.cpp,v 1.2
3 : *
4 : * Project: Terragen(tm) TER Driver
5 : * Purpose: Reader for Terragen TER documents
6 : * Author: Ray Gardener, Daylon Graphics Ltd.
7 : *
8 : * Portions of this module derived from GDAL drivers by
9 : * Frank Warmerdam, see http://www.gdal.org
10 :
11 : rcg apr 19/06 Fixed bug with hf size being misread by one
12 : if xpts/ypts tags not included in file.
13 : Added Create() support.
14 : Treat pixels as points.
15 :
16 : rcg jun 6/07 Better heightscale/baseheight determination
17 : when writing.
18 :
19 : *
20 : ******************************************************************************
21 : * Copyright (c) 2006-2007 Daylon Graphics Ltd.
22 : *
23 : * SPDX-License-Identifier: MIT
24 : ******************************************************************************
25 : */
26 :
27 : /*
28 : Terragen format notes:
29 :
30 : Based on official Planetside specs.
31 :
32 : All distances along all three axes are in
33 : terrain units, which are 30m by default.
34 : If a SCAL chunk is present, however, it
35 : can indicate something other than 30.
36 : Note that uniform scaling should be used.
37 :
38 : The offset (base height) in the ALTW chunk
39 : is in terrain units, and the scale (height scale)
40 : is a normalized value using unsigned 16-bit notation.
41 : The physical terrain value for a read pixel is
42 : hv' = hv * scale / 65536 + offset.
43 : It still needs to be scaled by SCAL to
44 : get to meters.
45 :
46 : For writing:
47 :
48 : SCAL = gridpost distance in meters
49 : hv_px = hv_m / SCAL
50 : span_px = span_m / SCAL
51 : offset = see TerragenDataset::write_header()
52 : scale = see TerragenDataset::write_header()
53 : physical hv =
54 : (hv_px - offset) * 65536.0/scale
55 :
56 : We tell callers that:
57 :
58 : Elevations are Int16 when reading,
59 : and Float32 when writing. We need logical
60 : elevations when writing so that we can
61 : encode them with as much precision as possible
62 : when going down to physical 16-bit ints.
63 : Implementing band::SetScale/SetOffset won't work because
64 : it requires callers to know format write details.
65 : So we've added two Create() options that let the
66 : caller tell us the span's logical extent, and with
67 : those two values we can convert to physical pixels.
68 :
69 : band::GetUnitType() returns meters.
70 : band::GetScale() returns SCAL * (scale/65536)
71 : band::GetOffset() returns SCAL * offset
72 : ds::GetSpatialRef() returns a local CS
73 : using meters.
74 : ds::GetGeoTransform() returns a scale matrix
75 : having SCAL sx,sy members.
76 :
77 : ds::SetGeoTransform() lets us establish the
78 : size of ground pixels.
79 : ds::_SetProjection() lets us establish what
80 : units ground measures are in (also needed
81 : to calc the size of ground pixels).
82 : band::SetUnitType() tells us what units
83 : the given Float32 elevations are in.
84 : band::SetScale() is unused.
85 : band::SetOffset() is unused.
86 : */
87 :
88 : #include "cpl_string.h"
89 : #include "gdal_frmts.h"
90 : #include "gdal_pam.h"
91 : #include "gdal_driver.h"
92 : #include "gdal_drivermanager.h"
93 : #include "gdal_openinfo.h"
94 : #include "gdal_cpp_functions.h"
95 : #include "ogr_spatialref.h"
96 :
97 : #include <cmath>
98 :
99 : #include <algorithm>
100 :
101 : const double kdEarthCircumPolar = 40007849;
102 : const double kdEarthCircumEquat = 40075004;
103 :
104 1 : static double average(double a, double b)
105 : {
106 1 : return 0.5 * (a + b);
107 : }
108 :
109 0 : static double degrees_to_radians(double d)
110 : {
111 0 : return d * (M_PI / 180);
112 : }
113 :
114 2 : static bool approx_equal(double a, double b)
115 : {
116 2 : const double epsilon = 1e-5;
117 2 : return std::abs(a - b) <= epsilon;
118 : }
119 :
120 : /************************************************************************/
121 : /* ==================================================================== */
122 : /* TerragenDataset */
123 : /* ==================================================================== */
124 : /************************************************************************/
125 :
126 : class TerragenRasterBand;
127 :
128 : class TerragenDataset final : public GDALPamDataset
129 : {
130 : friend class TerragenRasterBand;
131 :
132 : double m_dScale, m_dOffset,
133 : m_dSCAL, // 30.0 normally, from SCAL chunk
134 : m_dGroundScale, m_dMetersPerGroundUnit, m_dMetersPerElevUnit,
135 : m_dLogSpan[2], m_span_m[2], m_span_px[2];
136 :
137 : GDALGeoTransform m_gt{};
138 :
139 : VSILFILE *m_fp;
140 : vsi_l_offset m_nDataOffset;
141 :
142 : GInt16 m_nHeightScale;
143 : GInt16 m_nBaseHeight;
144 :
145 : char *m_pszFilename;
146 : OGRSpatialReference m_oSRS{};
147 : char m_szUnits[32];
148 :
149 : bool m_bIsGeo;
150 :
151 : int LoadFromFile();
152 :
153 : public:
154 : TerragenDataset();
155 : ~TerragenDataset() override;
156 :
157 : static GDALDataset *Open(GDALOpenInfo *);
158 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
159 : int nBandsIn, GDALDataType eType,
160 : char **papszOptions);
161 :
162 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
163 : CPLErr SetGeoTransform(const GDALGeoTransform >) override;
164 : const OGRSpatialReference *GetSpatialRef() const override;
165 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
166 :
167 : protected:
168 : bool get(GInt16 &);
169 : bool get(GUInt16 &);
170 : bool get(float &);
171 : bool put(GInt16);
172 : bool put(float);
173 :
174 9 : bool skip(size_t n)
175 : {
176 9 : return 0 == VSIFSeekL(m_fp, n, SEEK_CUR);
177 : }
178 :
179 1 : bool pad(size_t n)
180 : {
181 1 : return skip(n);
182 : }
183 :
184 : bool read_next_tag(char *);
185 : bool write_next_tag(const char *);
186 : static bool tag_is(const char *szTag, const char *);
187 :
188 : bool write_header(void);
189 : };
190 :
191 : /************************************************************************/
192 : /* ==================================================================== */
193 : /* TerragenRasterBand */
194 : /* ==================================================================== */
195 : /************************************************************************/
196 :
197 : class TerragenRasterBand final : public GDALPamRasterBand
198 : {
199 : friend class TerragenDataset;
200 :
201 : void *m_pvLine;
202 : bool m_bFirstTime;
203 :
204 : public:
205 : explicit TerragenRasterBand(TerragenDataset *);
206 :
207 10 : ~TerragenRasterBand() override
208 5 : {
209 5 : if (m_pvLine != nullptr)
210 5 : CPLFree(m_pvLine);
211 10 : }
212 :
213 : // Geomeasure support.
214 : CPLErr IReadBlock(int, int, void *) override;
215 : const char *GetUnitType() override;
216 : double GetOffset(int *pbSuccess = nullptr) override;
217 : double GetScale(int *pbSuccess = nullptr) override;
218 :
219 : CPLErr IWriteBlock(int, int, void *) override;
220 : CPLErr SetUnitType(const char *) override;
221 : };
222 :
223 : /************************************************************************/
224 : /* TerragenRasterBand() */
225 : /************************************************************************/
226 :
227 5 : TerragenRasterBand::TerragenRasterBand(TerragenDataset *poDSIn)
228 5 : : m_pvLine(CPLMalloc(sizeof(GInt16) * poDSIn->GetRasterXSize())),
229 5 : m_bFirstTime(true)
230 : {
231 5 : poDS = poDSIn;
232 5 : nBand = 1;
233 :
234 5 : eDataType = poDSIn->GetAccess() == GA_ReadOnly ? GDT_Int16 : GDT_Float32;
235 :
236 5 : nBlockXSize = poDSIn->GetRasterXSize();
237 5 : nBlockYSize = 1;
238 5 : }
239 :
240 : /************************************************************************/
241 : /* IReadBlock() */
242 : /************************************************************************/
243 :
244 40 : CPLErr TerragenRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
245 : void *pImage)
246 : {
247 : // CPLAssert( sizeof(float) == sizeof(GInt32) );
248 40 : CPLAssert(nBlockXOff == 0);
249 40 : CPLAssert(pImage != nullptr);
250 :
251 40 : TerragenDataset &ds = *cpl::down_cast<TerragenDataset *>(poDS);
252 :
253 : /* -------------------------------------------------------------------- */
254 : /* Seek to scanline.
255 : Terragen is a bottom-top format, so we have to
256 : invert the row location.
257 : -------------------------------------------------------------------- */
258 40 : const size_t rowbytes = nBlockXSize * sizeof(GInt16);
259 :
260 40 : if (0 != VSIFSeekL(ds.m_fp,
261 40 : ds.m_nDataOffset +
262 40 : (ds.GetRasterYSize() - 1 - nBlockYOff) * rowbytes,
263 : SEEK_SET))
264 : {
265 0 : CPLError(CE_Failure, CPLE_FileIO, "Terragen Seek failed:%s",
266 0 : VSIStrerror(errno));
267 0 : return CE_Failure;
268 : }
269 :
270 : /* -------------------------------------------------------------------- */
271 : /* Read the scanline into the line buffer. */
272 : /* -------------------------------------------------------------------- */
273 :
274 40 : if (VSIFReadL(pImage, rowbytes, 1, ds.m_fp) != 1)
275 : {
276 0 : CPLError(CE_Failure, CPLE_FileIO, "Terragen read failed:%s",
277 0 : VSIStrerror(errno));
278 0 : return CE_Failure;
279 : }
280 :
281 : /* -------------------------------------------------------------------- */
282 : /* Swap on MSB platforms. */
283 : /* -------------------------------------------------------------------- */
284 : #ifdef CPL_MSB
285 : GDALSwapWords(pImage, sizeof(GInt16), nRasterXSize, sizeof(GInt16));
286 : #endif
287 :
288 40 : return CE_None;
289 : }
290 :
291 : /************************************************************************/
292 : /* GetUnitType() */
293 : /************************************************************************/
294 0 : const char *TerragenRasterBand::GetUnitType()
295 : {
296 : // todo: Return elevation units.
297 : // For Terragen documents, it is the same as the ground units.
298 0 : TerragenDataset *poGDS = cpl::down_cast<TerragenDataset *>(poDS);
299 :
300 0 : return poGDS->m_szUnits;
301 : }
302 :
303 : /************************************************************************/
304 : /* GetScale() */
305 : /************************************************************************/
306 :
307 1 : double TerragenRasterBand::GetScale(int *pbSuccess)
308 : {
309 1 : const TerragenDataset &ds = *cpl::down_cast<TerragenDataset *>(poDS);
310 1 : if (pbSuccess != nullptr)
311 0 : *pbSuccess = TRUE;
312 :
313 1 : return ds.m_dScale;
314 : }
315 :
316 : /************************************************************************/
317 : /* GetOffset() */
318 : /************************************************************************/
319 :
320 1 : double TerragenRasterBand::GetOffset(int *pbSuccess)
321 : {
322 1 : const TerragenDataset &ds = *cpl::down_cast<TerragenDataset *>(poDS);
323 1 : if (pbSuccess != nullptr)
324 0 : *pbSuccess = TRUE;
325 :
326 1 : return ds.m_dOffset;
327 : }
328 :
329 : /************************************************************************/
330 : /* IWriteBlock() */
331 : /************************************************************************/
332 :
333 20 : CPLErr TerragenRasterBand::IWriteBlock(CPL_UNUSED int nBlockXOff,
334 : int nBlockYOff, void *pImage)
335 : {
336 20 : CPLAssert(nBlockXOff == 0);
337 20 : CPLAssert(pImage != nullptr);
338 20 : CPLAssert(m_pvLine != nullptr);
339 :
340 20 : const size_t pixelsize = sizeof(GInt16);
341 :
342 20 : TerragenDataset &ds = *cpl::down_cast<TerragenDataset *>(poDS);
343 20 : if (m_bFirstTime)
344 : {
345 1 : m_bFirstTime = false;
346 1 : ds.write_header();
347 1 : ds.m_nDataOffset = VSIFTellL(ds.m_fp);
348 : }
349 20 : const size_t rowbytes = nBlockXSize * pixelsize;
350 :
351 20 : GInt16 *pLine = reinterpret_cast<GInt16 *>(m_pvLine);
352 :
353 20 : if (0 == VSIFSeekL(ds.m_fp,
354 20 : ds.m_nDataOffset +
355 : // Terragen is Y inverted.
356 20 : (ds.GetRasterYSize() - 1 - nBlockYOff) * rowbytes,
357 : SEEK_SET))
358 : {
359 : // Convert each float32 to int16.
360 20 : float *pfImage = reinterpret_cast<float *>(pImage);
361 420 : for (size_t x = 0; x < static_cast<size_t>(nBlockXSize); x++)
362 : {
363 400 : const double f = pfImage[x] * ds.m_dMetersPerElevUnit / ds.m_dSCAL;
364 400 : const GInt16 hv = static_cast<GInt16>(
365 400 : (f - ds.m_nBaseHeight) * 65536.0 / ds.m_nHeightScale
366 : /*+ ds.m_nShift*/);
367 400 : pLine[x] = hv;
368 : }
369 :
370 : #ifdef CPL_MSB
371 : GDALSwapWords(m_pvLine, pixelsize, nBlockXSize, pixelsize);
372 : #endif
373 20 : if (1 == VSIFWriteL(m_pvLine, rowbytes, 1, ds.m_fp))
374 20 : return CE_None;
375 : }
376 :
377 0 : return CE_Failure;
378 : }
379 :
380 0 : CPLErr TerragenRasterBand::SetUnitType(const char *psz)
381 : {
382 0 : TerragenDataset &ds = *cpl::down_cast<TerragenDataset *>(poDS);
383 :
384 0 : if (EQUAL(psz, "m"))
385 0 : ds.m_dMetersPerElevUnit = 1.0;
386 0 : else if (EQUAL(psz, "ft"))
387 0 : ds.m_dMetersPerElevUnit = 0.3048;
388 0 : else if (EQUAL(psz, "sft"))
389 0 : ds.m_dMetersPerElevUnit = 1200.0 / 3937.0;
390 : else
391 0 : return CE_Failure;
392 :
393 0 : return CE_None;
394 : }
395 :
396 : /************************************************************************/
397 : /* ==================================================================== */
398 : /* TerragenDataset */
399 : /* ==================================================================== */
400 : /************************************************************************/
401 :
402 : /************************************************************************/
403 : /* TerragenDataset() */
404 : /************************************************************************/
405 :
406 54 : TerragenDataset::TerragenDataset()
407 : : m_dScale(0.0), m_dOffset(0.0), m_dSCAL(30.0), m_dGroundScale(0.0),
408 : m_dMetersPerGroundUnit(1.0), m_dMetersPerElevUnit(1.0), m_fp(nullptr),
409 : m_nDataOffset(0), m_nHeightScale(0), m_nBaseHeight(0),
410 54 : m_pszFilename(nullptr), m_bIsGeo(false)
411 : {
412 54 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
413 :
414 54 : m_dLogSpan[0] = 0.0;
415 54 : m_dLogSpan[1] = 0.0;
416 :
417 54 : m_gt[0] = 0.0;
418 54 : m_gt[1] = m_dSCAL;
419 54 : m_gt[2] = 0.0;
420 54 : m_gt[3] = 0.0;
421 54 : m_gt[4] = 0.0;
422 54 : m_gt[5] = m_dSCAL;
423 54 : m_span_m[0] = 0.0;
424 54 : m_span_m[1] = 0.0;
425 54 : m_span_px[0] = 0.0;
426 54 : m_span_px[1] = 0.0;
427 54 : memset(m_szUnits, 0, sizeof(m_szUnits));
428 54 : }
429 :
430 : /************************************************************************/
431 : /* ~TerragenDataset() */
432 : /************************************************************************/
433 :
434 108 : TerragenDataset::~TerragenDataset()
435 :
436 : {
437 54 : FlushCache(true);
438 :
439 54 : CPLFree(m_pszFilename);
440 :
441 54 : if (m_fp != nullptr)
442 5 : VSIFCloseL(m_fp);
443 108 : }
444 :
445 1 : bool TerragenDataset::write_header()
446 : {
447 : char szHeader[16];
448 : // cppcheck-suppress bufferNotZeroTerminated
449 1 : memcpy(szHeader, "TERRAGENTERRAIN ", sizeof(szHeader));
450 :
451 1 : if (1 != VSIFWriteL(reinterpret_cast<void *>(szHeader), sizeof(szHeader), 1,
452 : m_fp))
453 : {
454 0 : CPLError(CE_Failure, CPLE_FileIO,
455 : "Couldn't write to Terragen file %s.\n"
456 : "Is file system full?",
457 : m_pszFilename);
458 :
459 0 : return false;
460 : }
461 :
462 : // --------------------------------------------------------------------
463 : // Write out the heightfield dimensions, etc.
464 : // --------------------------------------------------------------------
465 :
466 1 : const int nXSize = GetRasterXSize();
467 1 : const int nYSize = GetRasterYSize();
468 :
469 1 : write_next_tag("SIZE");
470 1 : put(static_cast<GInt16>(std::min(nXSize, nYSize) - 1));
471 1 : pad(sizeof(GInt16));
472 :
473 1 : if (nXSize != nYSize)
474 : {
475 0 : write_next_tag("XPTS");
476 0 : put(static_cast<GInt16>(nXSize));
477 0 : pad(sizeof(GInt16));
478 0 : write_next_tag("YPTS");
479 0 : put(static_cast<GInt16>(nYSize));
480 0 : pad(sizeof(GInt16));
481 : }
482 :
483 1 : if (m_bIsGeo)
484 : {
485 : /*
486 : With a geographic projection (degrees),
487 : m_dGroundScale will be in degrees and
488 : m_dMetersPerGroundUnit is undefined.
489 : So we're going to estimate a m_dMetersPerGroundUnit
490 : value here (i.e., meters per degree).
491 :
492 : We figure out the degree size of one
493 : pixel, and then the latitude degrees
494 : of the heightfield's center. The circumference of
495 : the latitude's great circle lets us know how
496 : wide the pixel is in meters, and we
497 : average that with the pixel's meter breadth,
498 : which is based on the polar circumference.
499 : */
500 :
501 : /* const double m_dDegLongPerPixel =
502 : fabs(m_gt[1]); */
503 :
504 0 : const double m_dDegLatPerPixel = std::abs(m_gt[5]);
505 :
506 : /* const double m_dCenterLongitude =
507 : m_gt[0] +
508 : (0.5 * m_dDegLongPerPixel * (nXSize-1)); */
509 :
510 : const double m_dCenterLatitude =
511 0 : m_gt[3] + (0.5 * m_dDegLatPerPixel * (nYSize - 1));
512 :
513 : const double dLatCircum =
514 : kdEarthCircumEquat *
515 0 : std::sin(degrees_to_radians(90.0 - m_dCenterLatitude));
516 :
517 0 : const double dMetersPerDegLongitude = dLatCircum / 360;
518 : /* const double dMetersPerPixelX =
519 : (m_dDegLongPerPixel / 360) * dLatCircum; */
520 :
521 0 : const double dMetersPerDegLatitude = kdEarthCircumPolar / 360;
522 : /* const double dMetersPerPixelY =
523 : (m_dDegLatPerPixel / 360) * kdEarthCircumPolar; */
524 :
525 0 : m_dMetersPerGroundUnit =
526 0 : average(dMetersPerDegLongitude, dMetersPerDegLatitude);
527 : }
528 :
529 1 : m_dSCAL = m_dGroundScale * m_dMetersPerGroundUnit;
530 :
531 1 : if (m_dSCAL != 30.0)
532 : {
533 1 : const float sc = static_cast<float>(m_dSCAL);
534 1 : write_next_tag("SCAL");
535 1 : put(sc);
536 1 : put(sc);
537 1 : put(sc);
538 : }
539 :
540 1 : if (!write_next_tag("ALTW"))
541 : {
542 0 : CPLError(CE_Failure, CPLE_FileIO,
543 : "Couldn't write to Terragen file %s.\n"
544 : "Is file system full?",
545 : m_pszFilename);
546 :
547 0 : return false;
548 : }
549 :
550 : // Compute physical scales and offsets.
551 1 : m_span_m[0] = m_dLogSpan[0] * m_dMetersPerElevUnit;
552 1 : m_span_m[1] = m_dLogSpan[1] * m_dMetersPerElevUnit;
553 :
554 1 : m_span_px[0] = m_span_m[0] / m_dSCAL;
555 1 : m_span_px[1] = m_span_m[1] / m_dSCAL;
556 :
557 1 : const double span_px = m_span_px[1] - m_span_px[0];
558 1 : m_nHeightScale = static_cast<GInt16>(span_px);
559 1 : if (m_nHeightScale == 0)
560 0 : m_nHeightScale++;
561 :
562 : // TODO(schwehr): Make static functions.
563 : #define P2L_PX(n, hs, bh) (static_cast<double>(n) / 65536.0 * (hs) + (bh))
564 :
565 : #define L2P_PX(n, hs, bh) (static_cast<int>(((n) - (bh)) * 65536.0 / (hs)))
566 :
567 : // Increase the heightscale until the physical span
568 : // fits within a 16-bit range. The smaller the logical span,
569 : // the more necessary this becomes.
570 1 : int hs = m_nHeightScale;
571 1 : int bh = 0;
572 4 : for (; hs <= 32767; hs++)
573 : {
574 4 : double prevdelta = 1.0e30;
575 229383 : for (bh = -32768; bh <= 32767; bh++)
576 : {
577 229380 : const int nValley = L2P_PX(m_span_px[0], hs, bh);
578 229380 : if (nValley < -32768)
579 98293 : continue;
580 131087 : const int nPeak = L2P_PX(m_span_px[1], hs, bh);
581 131087 : if (nPeak > 32767)
582 131082 : continue;
583 :
584 : // now see how closely the baseheight gets
585 : // to the pixel span.
586 5 : const double d = P2L_PX(nValley, hs, bh);
587 5 : const double delta = std::abs(d - m_span_px[0]);
588 5 : if (delta < prevdelta) // Converging?
589 4 : prevdelta = delta;
590 : else
591 : {
592 : // We're diverging, so use the previous bh
593 : // and stop looking.
594 1 : bh--;
595 1 : break;
596 : }
597 : }
598 4 : if (bh != 32768)
599 1 : break;
600 : }
601 1 : if (hs == 32768)
602 : {
603 0 : CPLError(CE_Failure, CPLE_FileIO,
604 : "Couldn't write to Terragen file %s.\n"
605 : "Cannot find adequate heightscale/baseheight combination.",
606 : m_pszFilename);
607 :
608 0 : return false;
609 : }
610 :
611 1 : m_nHeightScale = static_cast<GInt16>(hs);
612 1 : m_nBaseHeight = static_cast<GInt16>(bh);
613 :
614 : // m_nHeightScale is the one that gives us the
615 : // widest use of the 16-bit space. However, there
616 : // might be larger heightscales that, even though
617 : // the reduce the space usage, give us a better fit
618 : // for preserving the span extents.
619 :
620 1 : return put(m_nHeightScale) && put(m_nBaseHeight);
621 : }
622 :
623 : /************************************************************************/
624 : /* get() */
625 : /************************************************************************/
626 :
627 8 : bool TerragenDataset::get(GInt16 &value)
628 : {
629 8 : if (1 == VSIFReadL(&value, sizeof(value), 1, m_fp))
630 : {
631 8 : CPL_LSBPTR16(&value);
632 8 : return true;
633 : }
634 0 : return false;
635 : }
636 :
637 4 : bool TerragenDataset::get(GUInt16 &value)
638 : {
639 4 : if (1 == VSIFReadL(&value, sizeof(value), 1, m_fp))
640 : {
641 4 : CPL_LSBPTR16(&value);
642 4 : return true;
643 : }
644 0 : return false;
645 : }
646 :
647 12 : bool TerragenDataset::get(float &value)
648 : {
649 12 : if (1 == VSIFReadL(&value, sizeof(value), 1, m_fp))
650 : {
651 12 : CPL_LSBPTR32(&value);
652 12 : return true;
653 : }
654 0 : return false;
655 : }
656 :
657 : /************************************************************************/
658 : /* put() */
659 : /************************************************************************/
660 :
661 3 : bool TerragenDataset::put(GInt16 n)
662 : {
663 3 : CPL_LSBPTR16(&n);
664 3 : return 1 == VSIFWriteL(&n, sizeof(n), 1, m_fp);
665 : }
666 :
667 3 : bool TerragenDataset::put(float f)
668 : {
669 3 : CPL_LSBPTR32(&f);
670 3 : return 1 == VSIFWriteL(&f, sizeof(f), 1, m_fp);
671 : }
672 :
673 : /************************************************************************/
674 : /* tag stuff */
675 : /************************************************************************/
676 :
677 16 : bool TerragenDataset::read_next_tag(char *szTag)
678 : {
679 16 : return 1 == VSIFReadL(szTag, 4, 1, m_fp);
680 : }
681 :
682 3 : bool TerragenDataset::write_next_tag(const char *szTag)
683 : {
684 3 : return 1 == VSIFWriteL(reinterpret_cast<void *>(const_cast<char *>(szTag)),
685 3 : 4, 1, m_fp);
686 : }
687 :
688 40 : bool TerragenDataset::tag_is(const char *szTag, const char *sz)
689 : {
690 40 : return 0 == memcmp(szTag, sz, 4);
691 : }
692 :
693 : /************************************************************************/
694 : /* LoadFromFile() */
695 : /************************************************************************/
696 :
697 4 : int TerragenDataset::LoadFromFile()
698 : {
699 4 : m_dSCAL = 30.0;
700 4 : m_nDataOffset = 0;
701 :
702 4 : if (0 != VSIFSeekL(m_fp, 16, SEEK_SET))
703 0 : return FALSE;
704 :
705 : char szTag[4];
706 4 : if (!read_next_tag(szTag) || !tag_is(szTag, "SIZE"))
707 0 : return FALSE;
708 :
709 : GUInt16 nSize;
710 4 : if (!get(nSize) || !skip(2))
711 0 : return FALSE;
712 :
713 : // Set dimensions to SIZE chunk. If we don't
714 : // encounter XPTS/YPTS chunks, we can assume
715 : // the terrain to be square.
716 4 : GUInt16 xpts = nSize + 1;
717 4 : GUInt16 ypts = nSize + 1;
718 :
719 12 : while (read_next_tag(szTag))
720 : {
721 8 : if (tag_is(szTag, "XPTS"))
722 : {
723 0 : get(xpts);
724 0 : if (xpts < nSize || !skip(2))
725 0 : return FALSE;
726 0 : continue;
727 : }
728 :
729 8 : if (tag_is(szTag, "YPTS"))
730 : {
731 0 : get(ypts);
732 0 : if (ypts < nSize || !skip(2))
733 0 : return FALSE;
734 0 : continue;
735 : }
736 :
737 8 : if (tag_is(szTag, "SCAL"))
738 : {
739 4 : float sc[3] = {0.0f};
740 4 : get(sc[0]);
741 4 : get(sc[1]);
742 4 : get(sc[2]);
743 4 : m_dSCAL = sc[1];
744 4 : continue;
745 : }
746 :
747 4 : if (tag_is(szTag, "CRAD"))
748 : {
749 0 : if (!skip(sizeof(float)))
750 0 : return FALSE;
751 0 : continue;
752 : }
753 4 : if (tag_is(szTag, "CRVM"))
754 : {
755 0 : if (!skip(sizeof(GUInt32)))
756 0 : return FALSE;
757 0 : continue;
758 : }
759 4 : if (tag_is(szTag, "ALTW"))
760 : {
761 4 : get(m_nHeightScale);
762 4 : get(m_nBaseHeight);
763 4 : m_nDataOffset = VSIFTellL(m_fp);
764 4 : if (!skip(static_cast<size_t>(xpts) * static_cast<size_t>(ypts) *
765 : sizeof(GInt16)))
766 0 : return FALSE;
767 4 : continue;
768 : }
769 0 : if (tag_is(szTag, "EOF "))
770 : {
771 0 : break;
772 : }
773 : }
774 :
775 4 : if (xpts == 0 || ypts == 0 || m_nDataOffset == 0)
776 0 : return FALSE;
777 :
778 4 : nRasterXSize = xpts;
779 4 : nRasterYSize = ypts;
780 :
781 : // todo: sanity check: do we have enough pixels?
782 :
783 : // Cache realworld scaling and offset.
784 4 : m_dScale = m_dSCAL / 65536 * m_nHeightScale;
785 4 : m_dOffset = m_dSCAL * m_nBaseHeight;
786 4 : strcpy(m_szUnits, "m");
787 :
788 : // Make our projection to have origin at the
789 : // NW corner, and groundscale to match elev scale
790 : // (i.e., uniform voxels).
791 4 : m_gt[0] = 0.0;
792 4 : m_gt[1] = m_dSCAL;
793 4 : m_gt[2] = 0.0;
794 4 : m_gt[3] = 0.0;
795 4 : m_gt[4] = 0.0;
796 4 : m_gt[5] = m_dSCAL;
797 :
798 : /* -------------------------------------------------------------------- */
799 : /* Set projection. */
800 : /* -------------------------------------------------------------------- */
801 : // Terragen files as of Apr 2006 are partially georeferenced,
802 : // we can declare a local coordsys that uses meters.
803 4 : m_oSRS.SetLocalCS("Terragen world space");
804 4 : m_oSRS.SetLinearUnits("m", 1.0);
805 :
806 4 : return TRUE;
807 : }
808 :
809 : /************************************************************************/
810 : /* SetSpatialRef() */
811 : /************************************************************************/
812 :
813 1 : CPLErr TerragenDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
814 : {
815 : // Terragen files aren't really georeferenced, but
816 : // we should get the projection's linear units so
817 : // that we can scale elevations correctly.
818 :
819 : // m_dSCAL = 30.0; // default
820 :
821 1 : m_oSRS.Clear();
822 1 : if (poSRS)
823 1 : m_oSRS = *poSRS;
824 :
825 : /* -------------------------------------------------------------------- */
826 : /* Linear units. */
827 : /* -------------------------------------------------------------------- */
828 1 : m_bIsGeo = poSRS != nullptr && m_oSRS.IsGeographic() != FALSE;
829 1 : if (m_bIsGeo)
830 : {
831 : // The caller is using degrees. We need to convert
832 : // to meters, otherwise we can't derive a SCAL
833 : // value to scale elevations with.
834 0 : m_bIsGeo = true;
835 : }
836 : else
837 : {
838 1 : const double dfLinear = m_oSRS.GetLinearUnits();
839 :
840 1 : if (approx_equal(dfLinear, 0.3048))
841 0 : m_dMetersPerGroundUnit = 0.3048;
842 1 : else if (approx_equal(dfLinear, CPLAtof(SRS_UL_US_FOOT_CONV)))
843 0 : m_dMetersPerGroundUnit = CPLAtof(SRS_UL_US_FOOT_CONV);
844 : else
845 1 : m_dMetersPerGroundUnit = 1.0;
846 : }
847 :
848 1 : return CE_None;
849 : }
850 :
851 : /************************************************************************/
852 : /* GetSpatialRef() */
853 : /************************************************************************/
854 :
855 1 : const OGRSpatialReference *TerragenDataset::GetSpatialRef() const
856 : {
857 1 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
858 : }
859 :
860 : /************************************************************************/
861 : /* SetGeoTransform() */
862 : /************************************************************************/
863 :
864 1 : CPLErr TerragenDataset::SetGeoTransform(const GDALGeoTransform >)
865 : {
866 1 : m_gt = gt;
867 :
868 : // Average the projection scales.
869 1 : m_dGroundScale = average(fabs(m_gt[1]), fabs(m_gt[5]));
870 1 : return CE_None;
871 : }
872 :
873 : /************************************************************************/
874 : /* GetGeoTransform() */
875 : /************************************************************************/
876 :
877 1 : CPLErr TerragenDataset::GetGeoTransform(GDALGeoTransform >) const
878 : {
879 1 : gt = m_gt;
880 1 : return CE_None;
881 : }
882 :
883 : /************************************************************************/
884 : /* Create() */
885 : /************************************************************************/
886 50 : GDALDataset *TerragenDataset::Create(const char *pszFilename, int nXSize,
887 : int nYSize, int nBandsIn,
888 : GDALDataType eType, char **papszOptions)
889 : {
890 50 : TerragenDataset *poDS = new TerragenDataset();
891 :
892 50 : poDS->eAccess = GA_Update;
893 :
894 50 : poDS->m_pszFilename = CPLStrdup(pszFilename);
895 :
896 : // --------------------------------------------------------------------
897 : // Verify input options.
898 : // --------------------------------------------------------------------
899 50 : const char *pszValue = CSLFetchNameValue(papszOptions, "MINUSERPIXELVALUE");
900 50 : if (pszValue != nullptr)
901 1 : poDS->m_dLogSpan[0] = CPLAtof(pszValue);
902 :
903 50 : pszValue = CSLFetchNameValue(papszOptions, "MAXUSERPIXELVALUE");
904 50 : if (pszValue != nullptr)
905 1 : poDS->m_dLogSpan[1] = CPLAtof(pszValue);
906 :
907 50 : if (poDS->m_dLogSpan[1] <= poDS->m_dLogSpan[0])
908 : {
909 49 : CPLError(CE_Failure, CPLE_AppDefined,
910 : "Inverted, flat, or unspecified span for Terragen file.");
911 :
912 49 : delete poDS;
913 49 : return nullptr;
914 : }
915 :
916 1 : if (eType != GDT_Float32)
917 : {
918 0 : CPLError(CE_Failure, CPLE_AppDefined,
919 : "Attempt to create Terragen dataset with a non-float32\n"
920 : "data type (%s).\n",
921 : GDALGetDataTypeName(eType));
922 :
923 0 : delete poDS;
924 0 : return nullptr;
925 : }
926 :
927 1 : if (nBandsIn != 1)
928 : {
929 0 : CPLError(CE_Failure, CPLE_NotSupported,
930 : "Terragen driver doesn't support %d bands. Must be 1.\n",
931 : nBandsIn);
932 :
933 0 : delete poDS;
934 0 : return nullptr;
935 : }
936 :
937 : // --------------------------------------------------------------------
938 : // Try to create the file.
939 : // --------------------------------------------------------------------
940 :
941 1 : poDS->m_fp = VSIFOpenL(pszFilename, "wb+");
942 :
943 1 : if (poDS->m_fp == nullptr)
944 : {
945 0 : CPLError(CE_Failure, CPLE_OpenFailed,
946 : "Attempt to create file `%s' failed.\n", pszFilename);
947 0 : delete poDS;
948 0 : return nullptr;
949 : }
950 :
951 1 : poDS->nRasterXSize = nXSize;
952 1 : poDS->nRasterYSize = nYSize;
953 :
954 : // Don't bother writing the header here; the first
955 : // call to IWriteBlock will do that instead, since
956 : // the elevation data's location depends on the
957 : // header size.
958 :
959 : // --------------------------------------------------------------------
960 : // Instance a band.
961 : // --------------------------------------------------------------------
962 1 : poDS->SetBand(1, new TerragenRasterBand(poDS));
963 :
964 : // VSIFClose( poDS->m_fp );
965 :
966 : // return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
967 1 : return GDALDataset::FromHandle(poDS);
968 : }
969 :
970 : /************************************************************************/
971 : /* Open() */
972 : /************************************************************************/
973 :
974 37812 : GDALDataset *TerragenDataset::Open(GDALOpenInfo *poOpenInfo)
975 :
976 : {
977 : // The file should have at least 32 header bytes
978 37812 : if (poOpenInfo->nHeaderBytes < 32 || poOpenInfo->fpL == nullptr)
979 32245 : return nullptr;
980 :
981 5567 : if (!EQUALN(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
982 : "TERRAGENTERRAIN ", 16))
983 5563 : return nullptr;
984 :
985 : /* -------------------------------------------------------------------- */
986 : /* Create a corresponding GDALDataset. */
987 : /* -------------------------------------------------------------------- */
988 4 : TerragenDataset *poDS = new TerragenDataset();
989 4 : poDS->eAccess = poOpenInfo->eAccess;
990 4 : poDS->m_fp = poOpenInfo->fpL;
991 4 : poOpenInfo->fpL = nullptr;
992 :
993 : /* -------------------------------------------------------------------- */
994 : /* Read the file. */
995 : /* -------------------------------------------------------------------- */
996 4 : if (!poDS->LoadFromFile())
997 : {
998 0 : delete poDS;
999 0 : return nullptr;
1000 : }
1001 :
1002 : /* -------------------------------------------------------------------- */
1003 : /* Create band information objects. */
1004 : /* -------------------------------------------------------------------- */
1005 4 : poDS->SetBand(1, new TerragenRasterBand(poDS));
1006 :
1007 4 : poDS->SetMetadataItem(GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT);
1008 :
1009 : /* -------------------------------------------------------------------- */
1010 : /* Initialize any PAM information. */
1011 : /* -------------------------------------------------------------------- */
1012 4 : poDS->SetDescription(poOpenInfo->pszFilename);
1013 4 : poDS->TryLoadXML();
1014 :
1015 : /* -------------------------------------------------------------------- */
1016 : /* Support overviews. */
1017 : /* -------------------------------------------------------------------- */
1018 4 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
1019 :
1020 4 : return poDS;
1021 : }
1022 :
1023 : /************************************************************************/
1024 : /* GDALRegister_Terragen() */
1025 : /************************************************************************/
1026 :
1027 2038 : void GDALRegister_Terragen()
1028 :
1029 : {
1030 2038 : if (GDALGetDriverByName("Terragen") != nullptr)
1031 283 : return;
1032 :
1033 1755 : GDALDriver *poDriver = new GDALDriver();
1034 :
1035 1755 : poDriver->SetDescription("Terragen");
1036 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1037 1755 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ter");
1038 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Terragen heightfield");
1039 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1040 1755 : "drivers/raster/terragen.html");
1041 :
1042 1755 : poDriver->SetMetadataItem(
1043 : GDAL_DMD_CREATIONOPTIONLIST,
1044 : "<CreationOptionList>"
1045 : " <Option name='MINUSERPIXELVALUE' type='float' description='Lowest "
1046 : "logical elevation'/>"
1047 : " <Option name='MAXUSERPIXELVALUE' type='float' description='Highest "
1048 : "logical elevation'/>"
1049 1755 : "</CreationOptionList>");
1050 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1051 :
1052 1755 : poDriver->pfnOpen = TerragenDataset::Open;
1053 1755 : poDriver->pfnCreate = TerragenDataset::Create;
1054 :
1055 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
1056 : }
|
