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