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