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