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