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