1 : /******************************************************************************
2 : * $Id: hfaopen.cpp 18348 2009-12-19 14:32:39Z rouault $
3 : *
4 : * Project: Erdas Imagine (.img) Translator
5 : * Purpose: Supporting functions for HFA (.img) ... main (C callable) API
6 : * that is not dependent on GDAL (just CPL).
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 1999, Intergraph Corporation
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ******************************************************************************
30 : *
31 : * hfaopen.cpp
32 : *
33 : * Supporting routines for reading Erdas Imagine (.imf) Heirarchical
34 : * File Architecture files. This is intended to be a library independent
35 : * of the GDAL core, but dependent on the Common Portability Library.
36 : *
37 : */
38 :
39 : #include "hfa_p.h"
40 : #include "cpl_conv.h"
41 : #include <limits.h>
42 :
43 : CPL_CVSID("$Id: hfaopen.cpp 18348 2009-12-19 14:32:39Z rouault $");
44 :
45 :
46 : static const char *apszAuxMetadataItems[] = {
47 :
48 : // node/entry field_name metadata_key type
49 :
50 : "Statistics", "dminimum", "STATISTICS_MINIMUM", "Esta_Statistics",
51 : "Statistics", "dmaximum", "STATISTICS_MAXIMUM", "Esta_Statistics",
52 : "Statistics", "dmean", "STATISTICS_MEAN", "Esta_Statistics",
53 : "Statistics", "dmedian", "STATISTICS_MEDIAN", "Esta_Statistics",
54 : "Statistics", "dmode", "STATISTICS_MODE", "Esta_Statistics",
55 : "Statistics", "dstddev", "STATISTICS_STDDEV", "Esta_Statistics",
56 : "HistogramParameters", "lBinFunction.numBins", "STATISTICS_HISTONUMBINS","Eimg_StatisticsParameters830",
57 : "HistogramParameters", "dBinFunction.minLimit", "STATISTICS_HISTOMIN", "Eimg_StatisticsParameters830",
58 : "HistogramParameters", "dBinFunction.maxLimit", "STATISTICS_HISTOMAX", "Eimg_StatisticsParameters830",
59 : "StatisticsParameters", "lSkipFactorX", "STATISTICS_SKIPFACTORX", "",
60 : "StatisticsParameters", "lSkipFactorY", "STATISTICS_SKIPFACTORY", "",
61 : "StatisticsParameters", "dExcludedValues", "STATISTICS_EXCLUDEDVALUES","",
62 : "", "elayerType", "LAYER_TYPE", "",
63 : NULL
64 : };
65 :
66 :
67 449 : const char ** GetHFAAuxMetaDataList()
68 : {
69 449 : return apszAuxMetadataItems;
70 : }
71 :
72 :
73 : /************************************************************************/
74 : /* HFAGetDictionary() */
75 : /************************************************************************/
76 :
77 280 : static char * HFAGetDictionary( HFAHandle hHFA )
78 :
79 : {
80 280 : int nDictMax = 100;
81 280 : char *pszDictionary = (char *) CPLMalloc(nDictMax);
82 280 : int nDictSize = 0;
83 :
84 280 : VSIFSeekL( hHFA->fp, hHFA->nDictionaryPos, SEEK_SET );
85 :
86 1031179 : while( TRUE )
87 : {
88 1031459 : if( nDictSize >= nDictMax-1 )
89 : {
90 1390 : nDictMax = nDictSize * 2 + 100;
91 1390 : pszDictionary = (char *) CPLRealloc(pszDictionary, nDictMax );
92 : }
93 :
94 3173761 : if( VSIFReadL( pszDictionary + nDictSize, 1, 1, hHFA->fp ) < 1
95 1031438 : || pszDictionary[nDictSize] == '\0'
96 1030340 : || (nDictSize > 2 && pszDictionary[nDictSize-2] == ','
97 80524 : && pszDictionary[nDictSize-1] == '.') )
98 : break;
99 :
100 1031179 : nDictSize++;
101 : }
102 :
103 280 : pszDictionary[nDictSize] = '\0';
104 :
105 :
106 280 : return( pszDictionary );
107 : }
108 :
109 : /************************************************************************/
110 : /* HFAOpen() */
111 : /************************************************************************/
112 :
113 280 : HFAHandle HFAOpen( const char * pszFilename, const char * pszAccess )
114 :
115 : {
116 : FILE *fp;
117 : char szHeader[16];
118 : HFAInfo_t *psInfo;
119 : GUInt32 nHeaderPos;
120 :
121 : /* -------------------------------------------------------------------- */
122 : /* Open the file. */
123 : /* -------------------------------------------------------------------- */
124 456 : if( EQUAL(pszAccess,"r") || EQUAL(pszAccess,"rb" ) )
125 176 : fp = VSIFOpenL( pszFilename, "rb" );
126 : else
127 104 : fp = VSIFOpenL( pszFilename, "r+b" );
128 :
129 : /* should this be changed to use some sort of CPLFOpen() which will
130 : set the error? */
131 280 : if( fp == NULL )
132 : {
133 : CPLError( CE_Failure, CPLE_OpenFailed,
134 : "File open of %s failed.",
135 0 : pszFilename );
136 :
137 0 : return NULL;
138 : }
139 :
140 : /* -------------------------------------------------------------------- */
141 : /* Read and verify the header. */
142 : /* -------------------------------------------------------------------- */
143 280 : if( VSIFReadL( szHeader, 16, 1, fp ) < 1 )
144 : {
145 : CPLError( CE_Failure, CPLE_AppDefined,
146 : "Attempt to read 16 byte header failed for\n%s.",
147 0 : pszFilename );
148 :
149 0 : return NULL;
150 : }
151 :
152 280 : if( !EQUALN(szHeader,"EHFA_HEADER_TAG",15) )
153 : {
154 : CPLError( CE_Failure, CPLE_AppDefined,
155 : "File %s is not an Imagine HFA file ... header wrong.",
156 0 : pszFilename );
157 :
158 0 : return NULL;
159 : }
160 :
161 : /* -------------------------------------------------------------------- */
162 : /* Create the HFAInfo_t */
163 : /* -------------------------------------------------------------------- */
164 280 : psInfo = (HFAInfo_t *) CPLCalloc(sizeof(HFAInfo_t),1);
165 :
166 280 : psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
167 280 : psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
168 280 : psInfo->fp = fp;
169 456 : if( EQUAL(pszAccess,"r") || EQUAL(pszAccess,"rb" ) )
170 176 : psInfo->eAccess = HFA_ReadOnly;
171 : else
172 104 : psInfo->eAccess = HFA_Update;
173 280 : psInfo->bTreeDirty = FALSE;
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Where is the header? */
177 : /* -------------------------------------------------------------------- */
178 280 : VSIFReadL( &nHeaderPos, sizeof(GInt32), 1, fp );
179 : HFAStandard( 4, &nHeaderPos );
180 :
181 : /* -------------------------------------------------------------------- */
182 : /* Read the header. */
183 : /* -------------------------------------------------------------------- */
184 280 : VSIFSeekL( fp, nHeaderPos, SEEK_SET );
185 :
186 280 : VSIFReadL( &(psInfo->nVersion), sizeof(GInt32), 1, fp );
187 : HFAStandard( 4, &(psInfo->nVersion) );
188 :
189 280 : VSIFReadL( szHeader, 4, 1, fp ); /* skip freeList */
190 :
191 280 : VSIFReadL( &(psInfo->nRootPos), sizeof(GInt32), 1, fp );
192 : HFAStandard( 4, &(psInfo->nRootPos) );
193 :
194 280 : VSIFReadL( &(psInfo->nEntryHeaderLength), sizeof(GInt16), 1, fp );
195 : HFAStandard( 2, &(psInfo->nEntryHeaderLength) );
196 :
197 280 : VSIFReadL( &(psInfo->nDictionaryPos), sizeof(GInt32), 1, fp );
198 : HFAStandard( 4, &(psInfo->nDictionaryPos) );
199 :
200 : /* -------------------------------------------------------------------- */
201 : /* Collect file size. */
202 : /* -------------------------------------------------------------------- */
203 280 : VSIFSeekL( fp, 0, SEEK_END );
204 280 : psInfo->nEndOfFile = (GUInt32) VSIFTellL( fp );
205 :
206 : /* -------------------------------------------------------------------- */
207 : /* Instantiate the root entry. */
208 : /* -------------------------------------------------------------------- */
209 280 : psInfo->poRoot = new HFAEntry( psInfo, psInfo->nRootPos, NULL, NULL );
210 :
211 : /* -------------------------------------------------------------------- */
212 : /* Read the dictionary */
213 : /* -------------------------------------------------------------------- */
214 280 : psInfo->pszDictionary = HFAGetDictionary( psInfo );
215 560 : psInfo->poDictionary = new HFADictionary( psInfo->pszDictionary );
216 :
217 : /* -------------------------------------------------------------------- */
218 : /* Collect band definitions. */
219 : /* -------------------------------------------------------------------- */
220 280 : HFAParseBandInfo( psInfo );
221 :
222 280 : return psInfo;
223 : }
224 :
225 : /************************************************************************/
226 : /* HFACreateDependent() */
227 : /* */
228 : /* Create a .rrd file for the named file if it does not exist, */
229 : /* or return the existing dependent if it already exists. */
230 : /************************************************************************/
231 :
232 1 : HFAInfo_t *HFACreateDependent( HFAInfo_t *psBase )
233 :
234 : {
235 1 : if( psBase->psDependent != NULL )
236 0 : return psBase->psDependent;
237 :
238 : /* -------------------------------------------------------------------- */
239 : /* Create desired RRD filename. */
240 : /* -------------------------------------------------------------------- */
241 1 : CPLString oBasename = CPLGetBasename( psBase->pszFilename );
242 : CPLString oRRDFilename =
243 1 : CPLFormFilename( psBase->pszPath, oBasename, "rrd" );
244 :
245 : /* -------------------------------------------------------------------- */
246 : /* Does this file already exist? If so, re-use it. */
247 : /* -------------------------------------------------------------------- */
248 1 : FILE *fp = VSIFOpenL( oRRDFilename, "rb" );
249 1 : if( fp != NULL )
250 : {
251 0 : VSIFCloseL( fp );
252 0 : psBase->psDependent = HFAOpen( oRRDFilename, "rb" );
253 : }
254 :
255 : /* -------------------------------------------------------------------- */
256 : /* Otherwise create it now. */
257 : /* -------------------------------------------------------------------- */
258 : HFAInfo_t *psDep;
259 1 : psDep = psBase->psDependent = HFACreateLL( oRRDFilename );
260 :
261 : /* -------------------------------------------------------------------- */
262 : /* Add the DependentFile node with the pointer back to the */
263 : /* parent. When working from an .aux file we really want the */
264 : /* .rrd to point back to the original file, not the .aux file. */
265 : /* -------------------------------------------------------------------- */
266 1 : HFAEntry *poEntry = psBase->poRoot->GetNamedChild("DependentFile");
267 1 : const char *pszDependentFile = NULL;
268 1 : if( poEntry != NULL )
269 0 : pszDependentFile = poEntry->GetStringField( "dependent.string" );
270 1 : if( pszDependentFile == NULL )
271 1 : pszDependentFile = psBase->pszFilename;
272 :
273 : HFAEntry *poDF = new HFAEntry( psDep, "DependentFile",
274 1 : "Eimg_DependentFile", psDep->poRoot );
275 :
276 1 : poDF->MakeData( strlen(pszDependentFile) + 50 );
277 1 : poDF->SetPosition();
278 1 : poDF->SetStringField( "dependent.string", pszDependentFile );
279 :
280 1 : return psDep;
281 : }
282 :
283 : /************************************************************************/
284 : /* HFAGetDependent() */
285 : /************************************************************************/
286 :
287 13 : HFAInfo_t *HFAGetDependent( HFAInfo_t *psBase, const char *pszFilename )
288 :
289 : {
290 13 : if( EQUAL(pszFilename,psBase->pszFilename) )
291 7 : return psBase;
292 :
293 6 : if( psBase->psDependent != NULL )
294 : {
295 1 : if( EQUAL(pszFilename,psBase->psDependent->pszFilename) )
296 1 : return psBase->psDependent;
297 : else
298 0 : return NULL;
299 : }
300 :
301 : /* -------------------------------------------------------------------- */
302 : /* Try to open the dependent file. */
303 : /* -------------------------------------------------------------------- */
304 : char *pszDependent;
305 : FILE *fp;
306 5 : const char* pszMode = psBase->eAccess == HFA_Update ? "r+b" : "rb";
307 :
308 : pszDependent = CPLStrdup(
309 5 : CPLFormFilename( psBase->pszPath, pszFilename, NULL ) );
310 :
311 5 : fp = VSIFOpenL( pszDependent, pszMode );
312 5 : if( fp != NULL )
313 : {
314 3 : VSIFCloseL( fp );
315 3 : psBase->psDependent = HFAOpen( pszDependent, pszMode );
316 : }
317 :
318 5 : CPLFree( pszDependent );
319 :
320 5 : return psBase->psDependent;
321 : }
322 :
323 :
324 : /************************************************************************/
325 : /* HFAParseBandInfo() */
326 : /* */
327 : /* This is used by HFAOpen() and HFACreate() to initialize the */
328 : /* band structures. */
329 : /************************************************************************/
330 :
331 375 : CPLErr HFAParseBandInfo( HFAInfo_t *psInfo )
332 :
333 : {
334 : HFAEntry *poNode;
335 :
336 : /* -------------------------------------------------------------------- */
337 : /* Find the first band node. */
338 : /* -------------------------------------------------------------------- */
339 375 : psInfo->nBands = 0;
340 375 : poNode = psInfo->poRoot->GetChild();
341 1755 : while( poNode != NULL )
342 : {
343 1005 : if( EQUAL(poNode->GetType(),"Eimg_Layer")
344 : && poNode->GetIntField("width") > 0
345 : && poNode->GetIntField("height") > 0 )
346 : {
347 548 : if( psInfo->nBands == 0 )
348 : {
349 368 : psInfo->nXSize = poNode->GetIntField("width");
350 368 : psInfo->nYSize = poNode->GetIntField("height");
351 : }
352 180 : else if( poNode->GetIntField("width") != psInfo->nXSize
353 : || poNode->GetIntField("height") != psInfo->nYSize )
354 : {
355 0 : return CE_Failure;
356 : }
357 :
358 : psInfo->papoBand = (HFABand **)
359 : CPLRealloc(psInfo->papoBand,
360 548 : sizeof(HFABand *) * (psInfo->nBands+1));
361 548 : psInfo->papoBand[psInfo->nBands] = new HFABand( psInfo, poNode );
362 548 : if (psInfo->papoBand[psInfo->nBands]->nWidth == 0)
363 : {
364 0 : delete psInfo->papoBand[psInfo->nBands];
365 0 : return CE_Failure;
366 : }
367 548 : psInfo->nBands++;
368 : }
369 :
370 1005 : poNode = poNode->GetNext();
371 : }
372 :
373 375 : return CE_None;
374 : }
375 :
376 : /************************************************************************/
377 : /* HFAClose() */
378 : /************************************************************************/
379 :
380 376 : void HFAClose( HFAHandle hHFA )
381 :
382 : {
383 : int i;
384 :
385 376 : if( hHFA->bTreeDirty || hHFA->poDictionary->bDictionaryTextDirty )
386 179 : HFAFlush( hHFA );
387 :
388 376 : if( hHFA->psDependent != NULL )
389 3 : HFAClose( hHFA->psDependent );
390 :
391 376 : delete hHFA->poRoot;
392 :
393 376 : VSIFCloseL( hHFA->fp );
394 :
395 376 : if( hHFA->poDictionary != NULL )
396 376 : delete hHFA->poDictionary;
397 :
398 376 : CPLFree( hHFA->pszDictionary );
399 376 : CPLFree( hHFA->pszFilename );
400 376 : CPLFree( hHFA->pszIGEFilename );
401 376 : CPLFree( hHFA->pszPath );
402 :
403 924 : for( i = 0; i < hHFA->nBands; i++ )
404 : {
405 548 : delete hHFA->papoBand[i];
406 : }
407 :
408 376 : CPLFree( hHFA->papoBand );
409 :
410 376 : if( hHFA->pProParameters != NULL )
411 : {
412 : Eprj_ProParameters *psProParms = (Eprj_ProParameters *)
413 32 : hHFA->pProParameters;
414 :
415 32 : CPLFree( psProParms->proExeName );
416 32 : CPLFree( psProParms->proName );
417 32 : CPLFree( psProParms->proSpheroid.sphereName );
418 :
419 32 : CPLFree( psProParms );
420 : }
421 :
422 376 : if( hHFA->pDatum != NULL )
423 : {
424 32 : CPLFree( ((Eprj_Datum *) hHFA->pDatum)->datumname );
425 32 : CPLFree( ((Eprj_Datum *) hHFA->pDatum)->gridname );
426 32 : CPLFree( hHFA->pDatum );
427 : }
428 :
429 376 : if( hHFA->pMapInfo != NULL )
430 : {
431 96 : CPLFree( ((Eprj_MapInfo *) hHFA->pMapInfo)->proName );
432 96 : CPLFree( ((Eprj_MapInfo *) hHFA->pMapInfo)->units );
433 96 : CPLFree( hHFA->pMapInfo );
434 : }
435 :
436 376 : CPLFree( hHFA );
437 376 : }
438 :
439 : /************************************************************************/
440 : /* HFARemove() */
441 : /* Used from HFADelete() function. */
442 : /************************************************************************/
443 :
444 0 : CPLErr HFARemove( const char *pszFilename )
445 :
446 : {
447 : VSIStatBufL sStat;
448 :
449 0 : if( VSIStatL( pszFilename, &sStat ) == 0 && VSI_ISREG( sStat.st_mode ) )
450 : {
451 0 : if( VSIUnlink( pszFilename ) == 0 )
452 0 : return CE_None;
453 : else
454 : {
455 : CPLError( CE_Failure, CPLE_AppDefined,
456 0 : "Attempt to unlink %s failed.\n", pszFilename );
457 0 : return CE_Failure;
458 : }
459 : }
460 : else
461 : {
462 : CPLError( CE_Failure, CPLE_AppDefined,
463 0 : "Unable to delete %s, not a file.\n", pszFilename );
464 0 : return CE_Failure;
465 : }
466 : }
467 :
468 : /************************************************************************/
469 : /* HFADelete() */
470 : /************************************************************************/
471 :
472 0 : CPLErr HFADelete( const char *pszFilename )
473 :
474 : {
475 0 : HFAInfo_t *psInfo = HFAOpen( pszFilename, "rb" );
476 0 : HFAEntry *poDMS = NULL;
477 0 : HFAEntry *poLayer = NULL;
478 0 : HFAEntry *poNode = NULL;
479 :
480 0 : if( psInfo != NULL )
481 : {
482 0 : poNode = psInfo->poRoot->GetChild();
483 0 : while( ( poNode != NULL ) && ( poLayer == NULL ) )
484 : {
485 0 : if( EQUAL(poNode->GetType(),"Eimg_Layer") )
486 : {
487 0 : poLayer = poNode;
488 : }
489 0 : poNode = poNode->GetNext();
490 : }
491 :
492 0 : if( poLayer != NULL )
493 0 : poDMS = poLayer->GetNamedChild( "ExternalRasterDMS" );
494 :
495 0 : if ( poDMS )
496 : {
497 : const char *pszRawFilename =
498 0 : poDMS->GetStringField( "fileName.string" );
499 :
500 0 : if( pszRawFilename != NULL )
501 : HFARemove( CPLFormFilename( psInfo->pszPath,
502 0 : pszRawFilename, NULL ) );
503 : }
504 :
505 0 : HFAClose( psInfo );
506 : }
507 0 : return HFARemove( pszFilename );
508 : }
509 :
510 : /************************************************************************/
511 : /* HFAGetRasterInfo() */
512 : /************************************************************************/
513 :
514 277 : CPLErr HFAGetRasterInfo( HFAHandle hHFA, int * pnXSize, int * pnYSize,
515 : int * pnBands )
516 :
517 : {
518 277 : if( pnXSize != NULL )
519 277 : *pnXSize = hHFA->nXSize;
520 277 : if( pnYSize != NULL )
521 277 : *pnYSize = hHFA->nYSize;
522 277 : if( pnBands != NULL )
523 277 : *pnBands = hHFA->nBands;
524 277 : return CE_None;
525 : }
526 :
527 : /************************************************************************/
528 : /* HFAGetBandInfo() */
529 : /************************************************************************/
530 :
531 419 : CPLErr HFAGetBandInfo( HFAHandle hHFA, int nBand, int * pnDataType,
532 : int * pnBlockXSize, int * pnBlockYSize,
533 : int *pnCompressionType )
534 :
535 : {
536 419 : if( nBand < 0 || nBand > hHFA->nBands )
537 : {
538 : CPLAssert( FALSE );
539 0 : return CE_Failure;
540 : }
541 :
542 419 : HFABand *poBand = hHFA->papoBand[nBand-1];
543 :
544 419 : if( pnDataType != NULL )
545 419 : *pnDataType = poBand->nDataType;
546 :
547 419 : if( pnBlockXSize != NULL )
548 419 : *pnBlockXSize = poBand->nBlockXSize;
549 :
550 419 : if( pnBlockYSize != NULL )
551 419 : *pnBlockYSize = poBand->nBlockYSize;
552 :
553 : /* -------------------------------------------------------------------- */
554 : /* Get compression code from RasterDMS. */
555 : /* -------------------------------------------------------------------- */
556 419 : if( pnCompressionType != NULL )
557 : {
558 : HFAEntry *poDMS;
559 :
560 419 : *pnCompressionType = 0;
561 :
562 419 : poDMS = poBand->poNode->GetNamedChild( "RasterDMS" );
563 :
564 419 : if( poDMS != NULL )
565 374 : *pnCompressionType = poDMS->GetIntField( "compressionType" );
566 : }
567 :
568 419 : return( CE_None );
569 : }
570 :
571 : /************************************************************************/
572 : /* HFAGetBandNoData() */
573 : /* */
574 : /* returns TRUE if value is set, otherwise FALSE. */
575 : /************************************************************************/
576 :
577 142 : int HFAGetBandNoData( HFAHandle hHFA, int nBand, double *pdfNoData )
578 :
579 : {
580 142 : if( nBand < 0 || nBand > hHFA->nBands )
581 : {
582 : CPLAssert( FALSE );
583 0 : return CE_Failure;
584 : }
585 :
586 142 : HFABand *poBand = hHFA->papoBand[nBand-1];
587 :
588 142 : *pdfNoData = poBand->dfNoData;
589 142 : return poBand->bNoDataSet;
590 : }
591 :
592 : /************************************************************************/
593 : /* HFASetBandNoData() */
594 : /* */
595 : /* attempts to set a no-data value on the given band */
596 : /************************************************************************/
597 :
598 4 : CPLErr HFASetBandNoData( HFAHandle hHFA, int nBand, double dfValue )
599 :
600 : {
601 4 : if ( nBand < 0 || nBand > hHFA->nBands )
602 : {
603 : CPLAssert( FALSE );
604 0 : return CE_Failure;
605 : }
606 :
607 4 : HFABand *poBand = hHFA->papoBand[nBand - 1];
608 :
609 4 : return poBand->SetNoDataValue( dfValue );
610 : }
611 :
612 : /************************************************************************/
613 : /* HFAGetOverviewCount() */
614 : /************************************************************************/
615 :
616 13 : int HFAGetOverviewCount( HFAHandle hHFA, int nBand )
617 :
618 : {
619 : HFABand *poBand;
620 :
621 13 : if( nBand < 0 || nBand > hHFA->nBands )
622 : {
623 : CPLAssert( FALSE );
624 0 : return CE_Failure;
625 : }
626 :
627 13 : poBand = hHFA->papoBand[nBand-1];
628 13 : poBand->LoadOverviews();
629 :
630 13 : return poBand->nOverviews;
631 : }
632 :
633 : /************************************************************************/
634 : /* HFAGetOverviewInfo() */
635 : /************************************************************************/
636 :
637 16 : CPLErr HFAGetOverviewInfo( HFAHandle hHFA, int nBand, int iOverview,
638 : int * pnXSize, int * pnYSize,
639 : int * pnBlockXSize, int * pnBlockYSize,
640 : int * pnHFADataType )
641 :
642 : {
643 : HFABand *poBand;
644 :
645 16 : if( nBand < 0 || nBand > hHFA->nBands )
646 : {
647 : CPLAssert( FALSE );
648 0 : return CE_Failure;
649 : }
650 :
651 16 : poBand = hHFA->papoBand[nBand-1];
652 16 : poBand->LoadOverviews();
653 :
654 16 : if( iOverview < 0 || iOverview >= poBand->nOverviews )
655 : {
656 : CPLAssert( FALSE );
657 0 : return CE_Failure;
658 : }
659 16 : poBand = poBand->papoOverviews[iOverview];
660 :
661 16 : if( pnXSize != NULL )
662 16 : *pnXSize = poBand->nWidth;
663 :
664 16 : if( pnYSize != NULL )
665 16 : *pnYSize = poBand->nHeight;
666 :
667 16 : if( pnBlockXSize != NULL )
668 16 : *pnBlockXSize = poBand->nBlockXSize;
669 :
670 16 : if( pnBlockYSize != NULL )
671 16 : *pnBlockYSize = poBand->nBlockYSize;
672 :
673 16 : if( pnHFADataType != NULL )
674 16 : *pnHFADataType = poBand->nDataType;
675 :
676 16 : return( CE_None );
677 : }
678 :
679 : /************************************************************************/
680 : /* HFAGetRasterBlock() */
681 : /************************************************************************/
682 :
683 0 : CPLErr HFAGetRasterBlock( HFAHandle hHFA, int nBand,
684 : int nXBlock, int nYBlock, void * pData )
685 :
686 : {
687 0 : return HFAGetRasterBlockEx(hHFA, nBand, nXBlock, nYBlock, pData, -1);
688 : }
689 :
690 : /************************************************************************/
691 : /* HFAGetRasterBlockEx() */
692 : /************************************************************************/
693 :
694 323 : CPLErr HFAGetRasterBlockEx( HFAHandle hHFA, int nBand,
695 : int nXBlock, int nYBlock, void * pData, int nDataSize )
696 :
697 : {
698 323 : if( nBand < 1 || nBand > hHFA->nBands )
699 0 : return CE_Failure;
700 :
701 323 : return( hHFA->papoBand[nBand-1]->GetRasterBlock(nXBlock,nYBlock,pData,nDataSize) );
702 : }
703 :
704 : /************************************************************************/
705 : /* HFAGetOverviewRasterBlock() */
706 : /************************************************************************/
707 :
708 0 : CPLErr HFAGetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
709 : int nXBlock, int nYBlock, void * pData )
710 :
711 : {
712 0 : return HFAGetOverviewRasterBlockEx(hHFA, nBand, iOverview, nXBlock, nYBlock, pData, -1);
713 : }
714 :
715 : /************************************************************************/
716 : /* HFAGetOverviewRasterBlockEx() */
717 : /************************************************************************/
718 :
719 12 : CPLErr HFAGetOverviewRasterBlockEx( HFAHandle hHFA, int nBand, int iOverview,
720 : int nXBlock, int nYBlock, void * pData, int nDataSize )
721 :
722 : {
723 12 : if( nBand < 1 || nBand > hHFA->nBands )
724 0 : return CE_Failure;
725 :
726 12 : if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand-1]->nOverviews )
727 0 : return CE_Failure;
728 :
729 24 : return( hHFA->papoBand[nBand-1]->papoOverviews[iOverview]->
730 24 : GetRasterBlock(nXBlock,nYBlock,pData, nDataSize) );
731 : }
732 :
733 : /************************************************************************/
734 : /* HFASetRasterBlock() */
735 : /************************************************************************/
736 :
737 111 : CPLErr HFASetRasterBlock( HFAHandle hHFA, int nBand,
738 : int nXBlock, int nYBlock, void * pData )
739 :
740 : {
741 111 : if( nBand < 1 || nBand > hHFA->nBands )
742 0 : return CE_Failure;
743 :
744 111 : return( hHFA->papoBand[nBand-1]->SetRasterBlock(nXBlock,nYBlock,pData) );
745 : }
746 :
747 : /************************************************************************/
748 : /* HFASetRasterBlock() */
749 : /************************************************************************/
750 :
751 8 : CPLErr HFASetOverviewRasterBlock( HFAHandle hHFA, int nBand, int iOverview,
752 : int nXBlock, int nYBlock, void * pData )
753 :
754 : {
755 8 : if( nBand < 1 || nBand > hHFA->nBands )
756 0 : return CE_Failure;
757 :
758 8 : if( iOverview < 0 || iOverview >= hHFA->papoBand[nBand-1]->nOverviews )
759 0 : return CE_Failure;
760 :
761 16 : return( hHFA->papoBand[nBand-1]->papoOverviews[iOverview]->
762 16 : SetRasterBlock(nXBlock,nYBlock,pData) );
763 : }
764 :
765 : /************************************************************************/
766 : /* HFAGetBandName() */
767 : /************************************************************************/
768 :
769 83 : const char * HFAGetBandName( HFAHandle hHFA, int nBand )
770 : {
771 83 : if( nBand < 1 || nBand > hHFA->nBands )
772 0 : return "";
773 :
774 83 : return( hHFA->papoBand[nBand-1]->GetBandName() );
775 : }
776 :
777 : /************************************************************************/
778 : /* HFASetBandName() */
779 : /************************************************************************/
780 :
781 111 : void HFASetBandName( HFAHandle hHFA, int nBand, const char *pszName )
782 : {
783 111 : if( nBand < 1 || nBand > hHFA->nBands )
784 0 : return;
785 :
786 111 : hHFA->papoBand[nBand-1]->SetBandName( pszName );
787 : }
788 :
789 : /************************************************************************/
790 : /* HFAGetDataTypeBits() */
791 : /************************************************************************/
792 :
793 1829 : int HFAGetDataTypeBits( int nDataType )
794 :
795 : {
796 1829 : switch( nDataType )
797 : {
798 : case EPT_u1:
799 29 : return 1;
800 :
801 : case EPT_u2:
802 0 : return 2;
803 :
804 : case EPT_u4:
805 0 : return 4;
806 :
807 : case EPT_u8:
808 : case EPT_s8:
809 794 : return 8;
810 :
811 : case EPT_u16:
812 : case EPT_s16:
813 230 : return 16;
814 :
815 : case EPT_u32:
816 : case EPT_s32:
817 : case EPT_f32:
818 225 : return 32;
819 :
820 : case EPT_f64:
821 : case EPT_c64:
822 484 : return 64;
823 :
824 : case EPT_c128:
825 67 : return 128;
826 : }
827 :
828 0 : return 0;
829 : }
830 :
831 : /************************************************************************/
832 : /* HFAGetDataTypeName() */
833 : /************************************************************************/
834 :
835 0 : const char *HFAGetDataTypeName( int nDataType )
836 :
837 : {
838 0 : switch( nDataType )
839 : {
840 : case EPT_u1:
841 0 : return "u1";
842 :
843 : case EPT_u2:
844 0 : return "u2";
845 :
846 : case EPT_u4:
847 0 : return "u4";
848 :
849 : case EPT_u8:
850 0 : return "u8";
851 :
852 : case EPT_s8:
853 0 : return "s8";
854 :
855 : case EPT_u16:
856 0 : return "u16";
857 :
858 : case EPT_s16:
859 0 : return "s16";
860 :
861 : case EPT_u32:
862 0 : return "u32";
863 :
864 : case EPT_s32:
865 0 : return "s32";
866 :
867 : case EPT_f32:
868 0 : return "f32";
869 :
870 : case EPT_f64:
871 0 : return "f64";
872 :
873 : case EPT_c64:
874 0 : return "c64";
875 :
876 : case EPT_c128:
877 0 : return "c128";
878 :
879 : default:
880 0 : return "unknown";
881 : }
882 : }
883 :
884 : /************************************************************************/
885 : /* HFAGetMapInfo() */
886 : /************************************************************************/
887 :
888 504 : const Eprj_MapInfo *HFAGetMapInfo( HFAHandle hHFA )
889 :
890 : {
891 : HFAEntry *poMIEntry;
892 : Eprj_MapInfo *psMapInfo;
893 :
894 504 : if( hHFA->nBands < 1 )
895 0 : return NULL;
896 :
897 : /* -------------------------------------------------------------------- */
898 : /* Do we already have it? */
899 : /* -------------------------------------------------------------------- */
900 504 : if( hHFA->pMapInfo != NULL )
901 52 : return( (Eprj_MapInfo *) hHFA->pMapInfo );
902 :
903 : /* -------------------------------------------------------------------- */
904 : /* Get the HFA node. */
905 : /* -------------------------------------------------------------------- */
906 452 : poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild( "Map_Info" );
907 452 : if( poMIEntry == NULL )
908 356 : return NULL;
909 :
910 : /* -------------------------------------------------------------------- */
911 : /* Allocate the structure. */
912 : /* -------------------------------------------------------------------- */
913 96 : psMapInfo = (Eprj_MapInfo *) CPLCalloc(sizeof(Eprj_MapInfo),1);
914 :
915 : /* -------------------------------------------------------------------- */
916 : /* Fetch the fields. */
917 : /* -------------------------------------------------------------------- */
918 96 : psMapInfo->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
919 :
920 : psMapInfo->upperLeftCenter.x =
921 96 : poMIEntry->GetDoubleField("upperLeftCenter.x");
922 : psMapInfo->upperLeftCenter.y =
923 96 : poMIEntry->GetDoubleField("upperLeftCenter.y");
924 :
925 : psMapInfo->lowerRightCenter.x =
926 96 : poMIEntry->GetDoubleField("lowerRightCenter.x");
927 : psMapInfo->lowerRightCenter.y =
928 96 : poMIEntry->GetDoubleField("lowerRightCenter.y");
929 :
930 : psMapInfo->pixelSize.width =
931 96 : poMIEntry->GetDoubleField("pixelSize.width");
932 : psMapInfo->pixelSize.height =
933 96 : poMIEntry->GetDoubleField("pixelSize.height");
934 :
935 96 : psMapInfo->units = CPLStrdup(poMIEntry->GetStringField("units"));
936 :
937 96 : hHFA->pMapInfo = (void *) psMapInfo;
938 :
939 96 : return psMapInfo;
940 : }
941 :
942 : /************************************************************************/
943 : /* HFAInvGeoTransform() */
944 : /************************************************************************/
945 :
946 6 : static int HFAInvGeoTransform( double *gt_in, double *gt_out )
947 :
948 : {
949 : double det, inv_det;
950 :
951 : /* we assume a 3rd row that is [1 0 0] */
952 :
953 : /* Compute determinate */
954 :
955 6 : det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
956 :
957 6 : if( fabs(det) < 0.000000000000001 )
958 0 : return 0;
959 :
960 6 : inv_det = 1.0 / det;
961 :
962 : /* compute adjoint, and devide by determinate */
963 :
964 6 : gt_out[1] = gt_in[5] * inv_det;
965 6 : gt_out[4] = -gt_in[4] * inv_det;
966 :
967 6 : gt_out[2] = -gt_in[2] * inv_det;
968 6 : gt_out[5] = gt_in[1] * inv_det;
969 :
970 6 : gt_out[0] = ( gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
971 6 : gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
972 :
973 6 : return 1;
974 : }
975 :
976 : /************************************************************************/
977 : /* HFAGetGeoTransform() */
978 : /************************************************************************/
979 :
980 275 : int HFAGetGeoTransform( HFAHandle hHFA, double *padfGeoTransform )
981 :
982 : {
983 275 : const Eprj_MapInfo *psMapInfo = HFAGetMapInfo( hHFA );
984 :
985 275 : padfGeoTransform[0] = 0.0;
986 275 : padfGeoTransform[1] = 1.0;
987 275 : padfGeoTransform[2] = 0.0;
988 275 : padfGeoTransform[3] = 0.0;
989 275 : padfGeoTransform[4] = 0.0;
990 275 : padfGeoTransform[5] = 1.0;
991 :
992 : /* -------------------------------------------------------------------- */
993 : /* Simple (north up) MapInfo approach. */
994 : /* -------------------------------------------------------------------- */
995 275 : if( psMapInfo != NULL )
996 : {
997 : padfGeoTransform[0] = psMapInfo->upperLeftCenter.x
998 96 : - psMapInfo->pixelSize.width*0.5;
999 96 : padfGeoTransform[1] = psMapInfo->pixelSize.width;
1000 96 : padfGeoTransform[2] = 0.0;
1001 96 : if( psMapInfo->upperLeftCenter.y >= psMapInfo->lowerRightCenter.y )
1002 96 : padfGeoTransform[5] = - psMapInfo->pixelSize.height;
1003 : else
1004 0 : padfGeoTransform[5] = psMapInfo->pixelSize.height;
1005 :
1006 96 : padfGeoTransform[3] = psMapInfo->upperLeftCenter.y
1007 96 : - padfGeoTransform[5]*0.5;
1008 96 : padfGeoTransform[4] = 0.0;
1009 :
1010 : // special logic to fixup odd angular units.
1011 96 : if( EQUAL(psMapInfo->units,"ds") )
1012 : {
1013 0 : padfGeoTransform[0] /= 3600.0;
1014 0 : padfGeoTransform[1] /= 3600.0;
1015 0 : padfGeoTransform[2] /= 3600.0;
1016 0 : padfGeoTransform[3] /= 3600.0;
1017 0 : padfGeoTransform[4] /= 3600.0;
1018 0 : padfGeoTransform[5] /= 3600.0;
1019 : }
1020 :
1021 96 : return TRUE;
1022 : }
1023 :
1024 : /* -------------------------------------------------------------------- */
1025 : /* Try for a MapToPixelXForm affine polynomial supporting */
1026 : /* rotated and sheared affine transformations. */
1027 : /* -------------------------------------------------------------------- */
1028 179 : if( hHFA->nBands == 0 )
1029 0 : return FALSE;
1030 :
1031 : HFAEntry *poXForm0 =
1032 179 : hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm.XForm0" );
1033 :
1034 179 : if( poXForm0 == NULL )
1035 174 : return FALSE;
1036 :
1037 5 : if( poXForm0->GetIntField( "order" ) != 1
1038 : || poXForm0->GetIntField( "numdimtransform" ) != 2
1039 : || poXForm0->GetIntField( "numdimpolynomial" ) != 2
1040 : || poXForm0->GetIntField( "termcount" ) != 3 )
1041 0 : return FALSE;
1042 :
1043 : // Verify that there aren't any further xform steps.
1044 5 : if( hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm.XForm1" )
1045 : != NULL )
1046 1 : return FALSE;
1047 :
1048 : // we should check that the exponent list is 0 0 1 0 0 1 but
1049 : // we don't because we are lazy
1050 :
1051 : // fetch geotransform values.
1052 : double adfXForm[6];
1053 :
1054 4 : adfXForm[0] = poXForm0->GetDoubleField( "polycoefvector[0]" );
1055 4 : adfXForm[1] = poXForm0->GetDoubleField( "polycoefmtx[0]" );
1056 4 : adfXForm[4] = poXForm0->GetDoubleField( "polycoefmtx[1]" );
1057 4 : adfXForm[3] = poXForm0->GetDoubleField( "polycoefvector[1]" );
1058 4 : adfXForm[2] = poXForm0->GetDoubleField( "polycoefmtx[2]" );
1059 4 : adfXForm[5] = poXForm0->GetDoubleField( "polycoefmtx[3]" );
1060 :
1061 : // invert
1062 :
1063 4 : HFAInvGeoTransform( adfXForm, padfGeoTransform );
1064 :
1065 : // Adjust origin from center of top left pixel to top left corner
1066 : // of top left pixel.
1067 :
1068 4 : padfGeoTransform[0] -= padfGeoTransform[1] * 0.5;
1069 4 : padfGeoTransform[0] -= padfGeoTransform[2] * 0.5;
1070 4 : padfGeoTransform[3] -= padfGeoTransform[4] * 0.5;
1071 4 : padfGeoTransform[3] -= padfGeoTransform[5] * 0.5;
1072 :
1073 4 : return TRUE;
1074 : }
1075 :
1076 : /************************************************************************/
1077 : /* HFASetMapInfo() */
1078 : /************************************************************************/
1079 :
1080 46 : CPLErr HFASetMapInfo( HFAHandle hHFA, const Eprj_MapInfo *poMapInfo )
1081 :
1082 : {
1083 : /* -------------------------------------------------------------------- */
1084 : /* Loop over bands, setting information on each one. */
1085 : /* -------------------------------------------------------------------- */
1086 102 : for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1087 : {
1088 : HFAEntry *poMIEntry;
1089 :
1090 : /* -------------------------------------------------------------------- */
1091 : /* Create a new Map_Info if there isn't one present already. */
1092 : /* -------------------------------------------------------------------- */
1093 56 : poMIEntry = hHFA->papoBand[iBand]->poNode->GetNamedChild( "Map_Info" );
1094 56 : if( poMIEntry == NULL )
1095 : {
1096 : poMIEntry = new HFAEntry( hHFA, "Map_Info", "Eprj_MapInfo",
1097 56 : hHFA->papoBand[iBand]->poNode );
1098 : }
1099 :
1100 56 : poMIEntry->MarkDirty();
1101 :
1102 : /* -------------------------------------------------------------------- */
1103 : /* Ensure we have enough space for all the data. */
1104 : /* -------------------------------------------------------------------- */
1105 : int nSize;
1106 : GByte *pabyData;
1107 :
1108 : nSize = 48 + 40
1109 : + strlen(poMapInfo->proName) + 1
1110 56 : + strlen(poMapInfo->units) + 1;
1111 :
1112 56 : pabyData = poMIEntry->MakeData( nSize );
1113 56 : memset( pabyData, 0, nSize );
1114 :
1115 56 : poMIEntry->SetPosition();
1116 :
1117 : /* -------------------------------------------------------------------- */
1118 : /* Write the various fields. */
1119 : /* -------------------------------------------------------------------- */
1120 56 : poMIEntry->SetStringField( "proName", poMapInfo->proName );
1121 :
1122 : poMIEntry->SetDoubleField( "upperLeftCenter.x",
1123 56 : poMapInfo->upperLeftCenter.x );
1124 : poMIEntry->SetDoubleField( "upperLeftCenter.y",
1125 56 : poMapInfo->upperLeftCenter.y );
1126 :
1127 : poMIEntry->SetDoubleField( "lowerRightCenter.x",
1128 56 : poMapInfo->lowerRightCenter.x );
1129 : poMIEntry->SetDoubleField( "lowerRightCenter.y",
1130 56 : poMapInfo->lowerRightCenter.y );
1131 :
1132 : poMIEntry->SetDoubleField( "pixelSize.width",
1133 56 : poMapInfo->pixelSize.width );
1134 : poMIEntry->SetDoubleField( "pixelSize.height",
1135 56 : poMapInfo->pixelSize.height );
1136 :
1137 56 : poMIEntry->SetStringField( "units", poMapInfo->units );
1138 : }
1139 :
1140 46 : return CE_None;
1141 : }
1142 :
1143 : /************************************************************************/
1144 : /* HFAGetPEString() */
1145 : /* */
1146 : /* Some files have a ProjectionX node contining the ESRI style */
1147 : /* PE_STRING. This function allows fetching from it. */
1148 : /************************************************************************/
1149 :
1150 275 : char *HFAGetPEString( HFAHandle hHFA )
1151 :
1152 : {
1153 275 : if( hHFA->nBands == 0 )
1154 0 : return NULL;
1155 :
1156 : /* -------------------------------------------------------------------- */
1157 : /* Get the HFA node. */
1158 : /* -------------------------------------------------------------------- */
1159 : HFAEntry *poProX;
1160 :
1161 275 : poProX = hHFA->papoBand[0]->poNode->GetNamedChild( "ProjectionX" );
1162 275 : if( poProX == NULL )
1163 229 : return NULL;
1164 :
1165 46 : const char *pszType = poProX->GetStringField( "projection.type.string" );
1166 46 : if( pszType == NULL || !EQUAL(pszType,"PE_COORDSYS") )
1167 0 : return NULL;
1168 :
1169 : /* -------------------------------------------------------------------- */
1170 : /* Use a gross hack to scan ahead to the actual projection */
1171 : /* string. We do it this way because we don't have general */
1172 : /* handling for MIFObjects. */
1173 : /* -------------------------------------------------------------------- */
1174 46 : GByte *pabyData = poProX->GetData();
1175 46 : int nDataSize = poProX->GetDataSize();
1176 :
1177 4324 : while( nDataSize > 10
1178 : && !EQUALN((const char *) pabyData,"PE_COORDSYS,.",13) ) {
1179 4232 : pabyData++;
1180 4232 : nDataSize--;
1181 : }
1182 :
1183 46 : if( nDataSize < 31 )
1184 0 : return NULL;
1185 :
1186 : /* -------------------------------------------------------------------- */
1187 : /* Skip ahead to the actual string. */
1188 : /* -------------------------------------------------------------------- */
1189 46 : pabyData += 30;
1190 46 : nDataSize -= 30;
1191 :
1192 46 : return CPLStrdup( (const char *) pabyData );
1193 : }
1194 :
1195 : /************************************************************************/
1196 : /* HFASetPEString() */
1197 : /************************************************************************/
1198 :
1199 36 : CPLErr HFASetPEString( HFAHandle hHFA, const char *pszPEString )
1200 :
1201 : {
1202 : /* -------------------------------------------------------------------- */
1203 : /* Loop over bands, setting information on each one. */
1204 : /* -------------------------------------------------------------------- */
1205 : int iBand;
1206 :
1207 82 : for( iBand = 0; iBand < hHFA->nBands; iBand++ )
1208 : {
1209 : HFAEntry *poProX;
1210 :
1211 : /* -------------------------------------------------------------------- */
1212 : /* Verify we don't already have the node, since update-in-place */
1213 : /* is likely to be more complicated. */
1214 : /* -------------------------------------------------------------------- */
1215 46 : poProX = hHFA->papoBand[iBand]->poNode->GetNamedChild( "ProjectionX" );
1216 :
1217 : /* -------------------------------------------------------------------- */
1218 : /* Create the node. */
1219 : /* -------------------------------------------------------------------- */
1220 46 : if( poProX == NULL )
1221 : {
1222 : poProX = new HFAEntry( hHFA, "ProjectionX","Eprj_MapProjection842",
1223 46 : hHFA->papoBand[iBand]->poNode );
1224 92 : if( poProX == NULL || poProX->GetTypeObject() == NULL )
1225 0 : return CE_Failure;
1226 : }
1227 :
1228 : /* -------------------------------------------------------------------- */
1229 : /* Prepare the data area with some extra space just in case. */
1230 : /* -------------------------------------------------------------------- */
1231 46 : GByte *pabyData = poProX->MakeData( 700 + strlen(pszPEString) );
1232 46 : if( !pabyData )
1233 0 : return CE_Failure;
1234 :
1235 46 : memset( pabyData, 0, 250+strlen(pszPEString) );
1236 :
1237 46 : poProX->SetPosition();
1238 :
1239 46 : poProX->SetStringField( "projection.type.string", "PE_COORDSYS" );
1240 : poProX->SetStringField( "projection.MIFDictionary.string",
1241 46 : "{0:pcstring,}Emif_String,{1:x{0:pcstring,}Emif_String,coordSys,}PE_COORDSYS,." );
1242 :
1243 : /* -------------------------------------------------------------------- */
1244 : /* Use a gross hack to scan ahead to the actual projection */
1245 : /* string. We do it this way because we don't have general */
1246 : /* handling for MIFObjects. */
1247 : /* -------------------------------------------------------------------- */
1248 46 : pabyData = poProX->GetData();
1249 46 : int nDataSize = poProX->GetDataSize();
1250 46 : GUInt32 iOffset = poProX->GetDataPos();
1251 : GUInt32 nSize;
1252 :
1253 4324 : while( nDataSize > 10
1254 : && !EQUALN((const char *) pabyData,"PE_COORDSYS,.",13) ) {
1255 4232 : pabyData++;
1256 4232 : nDataSize--;
1257 4232 : iOffset++;
1258 : }
1259 :
1260 : CPLAssert( nDataSize > (int) strlen(pszPEString) + 10 );
1261 :
1262 46 : pabyData += 14;
1263 46 : iOffset += 14;
1264 :
1265 : /* -------------------------------------------------------------------- */
1266 : /* Set the size and offset of the mifobject. */
1267 : /* -------------------------------------------------------------------- */
1268 46 : iOffset += 8;
1269 :
1270 46 : nSize = strlen(pszPEString) + 9;
1271 :
1272 : HFAStandard( 4, &nSize );
1273 46 : memcpy( pabyData, &nSize, 4 );
1274 46 : pabyData += 4;
1275 :
1276 : HFAStandard( 4, &iOffset );
1277 46 : memcpy( pabyData, &iOffset, 4 );
1278 46 : pabyData += 4;
1279 :
1280 : /* -------------------------------------------------------------------- */
1281 : /* Set the size and offset of the string value. */
1282 : /* -------------------------------------------------------------------- */
1283 46 : nSize = strlen(pszPEString) + 1;
1284 :
1285 : HFAStandard( 4, &nSize );
1286 46 : memcpy( pabyData, &nSize, 4 );
1287 46 : pabyData += 4;
1288 :
1289 46 : iOffset = 8;
1290 : HFAStandard( 4, &iOffset );
1291 46 : memcpy( pabyData, &iOffset, 4 );
1292 46 : pabyData += 4;
1293 :
1294 : /* -------------------------------------------------------------------- */
1295 : /* Place the string itself. */
1296 : /* -------------------------------------------------------------------- */
1297 46 : memcpy( pabyData, pszPEString, strlen(pszPEString)+1 );
1298 :
1299 46 : poProX->SetStringField( "title.string", "PE" );
1300 : }
1301 :
1302 36 : return CE_None;
1303 : }
1304 :
1305 : /************************************************************************/
1306 : /* HFAGetProParameters() */
1307 : /************************************************************************/
1308 :
1309 229 : const Eprj_ProParameters *HFAGetProParameters( HFAHandle hHFA )
1310 :
1311 : {
1312 : HFAEntry *poMIEntry;
1313 : Eprj_ProParameters *psProParms;
1314 : int i;
1315 :
1316 229 : if( hHFA->nBands < 1 )
1317 0 : return NULL;
1318 :
1319 : /* -------------------------------------------------------------------- */
1320 : /* Do we already have it? */
1321 : /* -------------------------------------------------------------------- */
1322 229 : if( hHFA->pProParameters != NULL )
1323 0 : return( (Eprj_ProParameters *) hHFA->pProParameters );
1324 :
1325 : /* -------------------------------------------------------------------- */
1326 : /* Get the HFA node. */
1327 : /* -------------------------------------------------------------------- */
1328 229 : poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild( "Projection" );
1329 229 : if( poMIEntry == NULL )
1330 197 : return NULL;
1331 :
1332 : /* -------------------------------------------------------------------- */
1333 : /* Allocate the structure. */
1334 : /* -------------------------------------------------------------------- */
1335 32 : psProParms = (Eprj_ProParameters *)CPLCalloc(sizeof(Eprj_ProParameters),1);
1336 :
1337 : /* -------------------------------------------------------------------- */
1338 : /* Fetch the fields. */
1339 : /* -------------------------------------------------------------------- */
1340 32 : psProParms->proType = (Eprj_ProType) poMIEntry->GetIntField("proType");
1341 32 : psProParms->proNumber = poMIEntry->GetIntField("proNumber");
1342 32 : psProParms->proExeName =CPLStrdup(poMIEntry->GetStringField("proExeName"));
1343 32 : psProParms->proName = CPLStrdup(poMIEntry->GetStringField("proName"));
1344 32 : psProParms->proZone = poMIEntry->GetIntField("proZone");
1345 :
1346 512 : for( i = 0; i < 15; i++ )
1347 : {
1348 : char szFieldName[30];
1349 :
1350 480 : sprintf( szFieldName, "proParams[%d]", i );
1351 480 : psProParms->proParams[i] = poMIEntry->GetDoubleField(szFieldName);
1352 : }
1353 :
1354 : psProParms->proSpheroid.sphereName =
1355 32 : CPLStrdup(poMIEntry->GetStringField("proSpheroid.sphereName"));
1356 32 : psProParms->proSpheroid.a = poMIEntry->GetDoubleField("proSpheroid.a");
1357 32 : psProParms->proSpheroid.b = poMIEntry->GetDoubleField("proSpheroid.b");
1358 : psProParms->proSpheroid.eSquared =
1359 32 : poMIEntry->GetDoubleField("proSpheroid.eSquared");
1360 : psProParms->proSpheroid.radius =
1361 32 : poMIEntry->GetDoubleField("proSpheroid.radius");
1362 :
1363 32 : hHFA->pProParameters = (void *) psProParms;
1364 :
1365 32 : return psProParms;
1366 : }
1367 :
1368 : /************************************************************************/
1369 : /* HFASetProParameters() */
1370 : /************************************************************************/
1371 :
1372 33 : CPLErr HFASetProParameters( HFAHandle hHFA, const Eprj_ProParameters *poPro )
1373 :
1374 : {
1375 : /* -------------------------------------------------------------------- */
1376 : /* Loop over bands, setting information on each one. */
1377 : /* -------------------------------------------------------------------- */
1378 76 : for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1379 : {
1380 : HFAEntry *poMIEntry;
1381 :
1382 : /* -------------------------------------------------------------------- */
1383 : /* Create a new Projection if there isn't one present already. */
1384 : /* -------------------------------------------------------------------- */
1385 43 : poMIEntry = hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1386 43 : if( poMIEntry == NULL )
1387 : {
1388 : poMIEntry = new HFAEntry( hHFA, "Projection","Eprj_ProParameters",
1389 43 : hHFA->papoBand[iBand]->poNode );
1390 : }
1391 :
1392 43 : poMIEntry->MarkDirty();
1393 :
1394 : /* -------------------------------------------------------------------- */
1395 : /* Ensure we have enough space for all the data. */
1396 : /* -------------------------------------------------------------------- */
1397 : int nSize;
1398 : GByte *pabyData;
1399 :
1400 : nSize = 34 + 15 * 8
1401 : + 8 + strlen(poPro->proName) + 1
1402 43 : + 32 + 8 + strlen(poPro->proSpheroid.sphereName) + 1;
1403 :
1404 43 : if( poPro->proExeName != NULL )
1405 0 : nSize += strlen(poPro->proExeName) + 1;
1406 :
1407 43 : pabyData = poMIEntry->MakeData( nSize );
1408 43 : poMIEntry->SetPosition();
1409 :
1410 : /* -------------------------------------------------------------------- */
1411 : /* Write the various fields. */
1412 : /* -------------------------------------------------------------------- */
1413 43 : poMIEntry->SetIntField( "proType", poPro->proType );
1414 :
1415 43 : poMIEntry->SetIntField( "proNumber", poPro->proNumber );
1416 :
1417 43 : poMIEntry->SetStringField( "proExeName", poPro->proExeName );
1418 43 : poMIEntry->SetStringField( "proName", poPro->proName );
1419 43 : poMIEntry->SetIntField( "proZone", poPro->proZone );
1420 43 : poMIEntry->SetDoubleField( "proParams[0]", poPro->proParams[0] );
1421 43 : poMIEntry->SetDoubleField( "proParams[1]", poPro->proParams[1] );
1422 43 : poMIEntry->SetDoubleField( "proParams[2]", poPro->proParams[2] );
1423 43 : poMIEntry->SetDoubleField( "proParams[3]", poPro->proParams[3] );
1424 43 : poMIEntry->SetDoubleField( "proParams[4]", poPro->proParams[4] );
1425 43 : poMIEntry->SetDoubleField( "proParams[5]", poPro->proParams[5] );
1426 43 : poMIEntry->SetDoubleField( "proParams[6]", poPro->proParams[6] );
1427 43 : poMIEntry->SetDoubleField( "proParams[7]", poPro->proParams[7] );
1428 43 : poMIEntry->SetDoubleField( "proParams[8]", poPro->proParams[8] );
1429 43 : poMIEntry->SetDoubleField( "proParams[9]", poPro->proParams[9] );
1430 43 : poMIEntry->SetDoubleField( "proParams[10]", poPro->proParams[10] );
1431 43 : poMIEntry->SetDoubleField( "proParams[11]", poPro->proParams[11] );
1432 43 : poMIEntry->SetDoubleField( "proParams[12]", poPro->proParams[12] );
1433 43 : poMIEntry->SetDoubleField( "proParams[13]", poPro->proParams[13] );
1434 43 : poMIEntry->SetDoubleField( "proParams[14]", poPro->proParams[14] );
1435 : poMIEntry->SetStringField( "proSpheroid.sphereName",
1436 43 : poPro->proSpheroid.sphereName );
1437 : poMIEntry->SetDoubleField( "proSpheroid.a",
1438 43 : poPro->proSpheroid.a );
1439 : poMIEntry->SetDoubleField( "proSpheroid.b",
1440 43 : poPro->proSpheroid.b );
1441 : poMIEntry->SetDoubleField( "proSpheroid.eSquared",
1442 43 : poPro->proSpheroid.eSquared );
1443 : poMIEntry->SetDoubleField( "proSpheroid.radius",
1444 43 : poPro->proSpheroid.radius );
1445 : }
1446 :
1447 33 : return CE_None;
1448 : }
1449 :
1450 : /************************************************************************/
1451 : /* HFAGetDatum() */
1452 : /************************************************************************/
1453 :
1454 229 : const Eprj_Datum *HFAGetDatum( HFAHandle hHFA )
1455 :
1456 : {
1457 : HFAEntry *poMIEntry;
1458 : Eprj_Datum *psDatum;
1459 : int i;
1460 :
1461 229 : if( hHFA->nBands < 1 )
1462 0 : return NULL;
1463 :
1464 : /* -------------------------------------------------------------------- */
1465 : /* Do we already have it? */
1466 : /* -------------------------------------------------------------------- */
1467 229 : if( hHFA->pDatum != NULL )
1468 0 : return( (Eprj_Datum *) hHFA->pDatum );
1469 :
1470 : /* -------------------------------------------------------------------- */
1471 : /* Get the HFA node. */
1472 : /* -------------------------------------------------------------------- */
1473 229 : poMIEntry = hHFA->papoBand[0]->poNode->GetNamedChild( "Projection.Datum" );
1474 229 : if( poMIEntry == NULL )
1475 197 : return NULL;
1476 :
1477 : /* -------------------------------------------------------------------- */
1478 : /* Allocate the structure. */
1479 : /* -------------------------------------------------------------------- */
1480 32 : psDatum = (Eprj_Datum *) CPLCalloc(sizeof(Eprj_Datum),1);
1481 :
1482 : /* -------------------------------------------------------------------- */
1483 : /* Fetch the fields. */
1484 : /* -------------------------------------------------------------------- */
1485 32 : psDatum->datumname = CPLStrdup(poMIEntry->GetStringField("datumname"));
1486 32 : psDatum->type = (Eprj_DatumType) poMIEntry->GetIntField("type");
1487 :
1488 256 : for( i = 0; i < 7; i++ )
1489 : {
1490 : char szFieldName[30];
1491 :
1492 224 : sprintf( szFieldName, "params[%d]", i );
1493 224 : psDatum->params[i] = poMIEntry->GetDoubleField(szFieldName);
1494 : }
1495 :
1496 32 : psDatum->gridname = CPLStrdup(poMIEntry->GetStringField("gridname"));
1497 :
1498 32 : hHFA->pDatum = (void *) psDatum;
1499 :
1500 32 : return psDatum;
1501 : }
1502 :
1503 : /************************************************************************/
1504 : /* HFASetDatum() */
1505 : /************************************************************************/
1506 :
1507 33 : CPLErr HFASetDatum( HFAHandle hHFA, const Eprj_Datum *poDatum )
1508 :
1509 : {
1510 : /* -------------------------------------------------------------------- */
1511 : /* Loop over bands, setting information on each one. */
1512 : /* -------------------------------------------------------------------- */
1513 76 : for( int iBand = 0; iBand < hHFA->nBands; iBand++ )
1514 : {
1515 43 : HFAEntry *poDatumEntry=NULL, *poProParms;
1516 :
1517 : /* -------------------------------------------------------------------- */
1518 : /* Create a new Projection if there isn't one present already. */
1519 : /* -------------------------------------------------------------------- */
1520 : poProParms =
1521 43 : hHFA->papoBand[iBand]->poNode->GetNamedChild("Projection");
1522 43 : if( poProParms == NULL )
1523 : {
1524 : CPLError( CE_Failure, CPLE_AppDefined,
1525 0 : "Can't add Eprj_Datum with no Eprj_ProjParameters." );
1526 0 : return CE_Failure;
1527 : }
1528 :
1529 43 : poDatumEntry = poProParms->GetNamedChild("Datum");
1530 43 : if( poDatumEntry == NULL )
1531 : {
1532 : poDatumEntry = new HFAEntry( hHFA, "Datum","Eprj_Datum",
1533 43 : poProParms );
1534 : }
1535 :
1536 43 : poDatumEntry->MarkDirty();
1537 :
1538 : /* -------------------------------------------------------------------- */
1539 : /* Ensure we have enough space for all the data. */
1540 : /* -------------------------------------------------------------------- */
1541 : int nSize;
1542 : GByte *pabyData;
1543 :
1544 43 : nSize = 26 + strlen(poDatum->datumname) + 1 + 7*8;
1545 :
1546 43 : if( poDatum->gridname != NULL )
1547 19 : nSize += strlen(poDatum->gridname) + 1;
1548 :
1549 43 : pabyData = poDatumEntry->MakeData( nSize );
1550 43 : poDatumEntry->SetPosition();
1551 :
1552 : /* -------------------------------------------------------------------- */
1553 : /* Write the various fields. */
1554 : /* -------------------------------------------------------------------- */
1555 43 : poDatumEntry->SetStringField( "datumname", poDatum->datumname );
1556 43 : poDatumEntry->SetIntField( "type", poDatum->type );
1557 :
1558 43 : poDatumEntry->SetDoubleField( "params[0]", poDatum->params[0] );
1559 43 : poDatumEntry->SetDoubleField( "params[1]", poDatum->params[1] );
1560 43 : poDatumEntry->SetDoubleField( "params[2]", poDatum->params[2] );
1561 43 : poDatumEntry->SetDoubleField( "params[3]", poDatum->params[3] );
1562 43 : poDatumEntry->SetDoubleField( "params[4]", poDatum->params[4] );
1563 43 : poDatumEntry->SetDoubleField( "params[5]", poDatum->params[5] );
1564 43 : poDatumEntry->SetDoubleField( "params[6]", poDatum->params[6] );
1565 :
1566 43 : poDatumEntry->SetStringField( "gridname", poDatum->gridname );
1567 : }
1568 :
1569 33 : return CE_None;
1570 : }
1571 :
1572 : /************************************************************************/
1573 : /* HFAGetPCT() */
1574 : /* */
1575 : /* Read the PCT from a band, if it has one. */
1576 : /************************************************************************/
1577 :
1578 403 : CPLErr HFAGetPCT( HFAHandle hHFA, int nBand, int *pnColors,
1579 : double **ppadfRed, double **ppadfGreen,
1580 : double **ppadfBlue , double **ppadfAlpha,
1581 : double **ppadfBins )
1582 :
1583 : {
1584 403 : if( nBand < 1 || nBand > hHFA->nBands )
1585 0 : return CE_Failure;
1586 :
1587 403 : return( hHFA->papoBand[nBand-1]->GetPCT( pnColors, ppadfRed,
1588 : ppadfGreen, ppadfBlue,
1589 806 : ppadfAlpha, ppadfBins ) );
1590 : }
1591 :
1592 : /************************************************************************/
1593 : /* HFASetPCT() */
1594 : /* */
1595 : /* Set the PCT on a band. */
1596 : /************************************************************************/
1597 :
1598 3 : CPLErr HFASetPCT( HFAHandle hHFA, int nBand, int nColors,
1599 : double *padfRed, double *padfGreen, double *padfBlue,
1600 : double *padfAlpha )
1601 :
1602 : {
1603 3 : if( nBand < 1 || nBand > hHFA->nBands )
1604 0 : return CE_Failure;
1605 :
1606 3 : return( hHFA->papoBand[nBand-1]->SetPCT( nColors, padfRed,
1607 6 : padfGreen, padfBlue, padfAlpha ) );
1608 : }
1609 :
1610 : /************************************************************************/
1611 : /* HFAGetDataRange() */
1612 : /************************************************************************/
1613 :
1614 0 : CPLErr HFAGetDataRange( HFAHandle hHFA, int nBand,
1615 : double * pdfMin, double *pdfMax )
1616 :
1617 : {
1618 : HFAEntry *poBinInfo;
1619 :
1620 0 : if( nBand < 1 || nBand > hHFA->nBands )
1621 0 : return CE_Failure;
1622 :
1623 0 : poBinInfo = hHFA->papoBand[nBand-1]->poNode->GetNamedChild("Statistics" );
1624 :
1625 0 : if( poBinInfo == NULL )
1626 0 : return( CE_Failure );
1627 :
1628 0 : *pdfMin = poBinInfo->GetDoubleField( "minimum" );
1629 0 : *pdfMax = poBinInfo->GetDoubleField( "maximum" );
1630 :
1631 0 : if( *pdfMax > *pdfMin )
1632 0 : return CE_None;
1633 : else
1634 0 : return CE_Failure;
1635 : }
1636 :
1637 : /************************************************************************/
1638 : /* HFADumpNode() */
1639 : /************************************************************************/
1640 :
1641 0 : static void HFADumpNode( HFAEntry *poEntry, int nIndent, int bVerbose,
1642 : FILE * fp )
1643 :
1644 : {
1645 : static char szSpaces[256];
1646 : int i;
1647 :
1648 0 : for( i = 0; i < nIndent*2; i++ )
1649 0 : szSpaces[i] = ' ';
1650 0 : szSpaces[nIndent*2] = '\0';
1651 :
1652 : fprintf( fp, "%s%s(%s) @ %d + %d @ %d\n", szSpaces,
1653 : poEntry->GetName(), poEntry->GetType(),
1654 : poEntry->GetFilePos(),
1655 0 : poEntry->GetDataSize(), poEntry->GetDataPos() );
1656 :
1657 0 : if( bVerbose )
1658 : {
1659 0 : strcat( szSpaces, "+ " );
1660 0 : poEntry->DumpFieldValues( fp, szSpaces );
1661 0 : fprintf( fp, "\n" );
1662 : }
1663 :
1664 0 : if( poEntry->GetChild() != NULL )
1665 0 : HFADumpNode( poEntry->GetChild(), nIndent+1, bVerbose, fp );
1666 :
1667 0 : if( poEntry->GetNext() != NULL )
1668 0 : HFADumpNode( poEntry->GetNext(), nIndent, bVerbose, fp );
1669 0 : }
1670 :
1671 : /************************************************************************/
1672 : /* HFADumpTree() */
1673 : /* */
1674 : /* Dump the tree of information in a HFA file. */
1675 : /************************************************************************/
1676 :
1677 0 : void HFADumpTree( HFAHandle hHFA, FILE * fpOut )
1678 :
1679 : {
1680 0 : HFADumpNode( hHFA->poRoot, 0, TRUE, fpOut );
1681 0 : }
1682 :
1683 : /************************************************************************/
1684 : /* HFADumpDictionary() */
1685 : /* */
1686 : /* Dump the dictionary (in raw, and parsed form) to the named */
1687 : /* device. */
1688 : /************************************************************************/
1689 :
1690 0 : void HFADumpDictionary( HFAHandle hHFA, FILE * fpOut )
1691 :
1692 : {
1693 0 : fprintf( fpOut, "%s\n", hHFA->pszDictionary );
1694 :
1695 0 : hHFA->poDictionary->Dump( fpOut );
1696 0 : }
1697 :
1698 : /************************************************************************/
1699 : /* HFAStandard() */
1700 : /* */
1701 : /* Swap byte order on MSB systems. */
1702 : /************************************************************************/
1703 :
1704 : #ifdef CPL_MSB
1705 : void HFAStandard( int nBytes, void * pData )
1706 :
1707 : {
1708 : int i;
1709 : GByte *pabyData = (GByte *) pData;
1710 :
1711 : for( i = nBytes/2-1; i >= 0; i-- )
1712 : {
1713 : GByte byTemp;
1714 :
1715 : byTemp = pabyData[i];
1716 : pabyData[i] = pabyData[nBytes-i-1];
1717 : pabyData[nBytes-i-1] = byTemp;
1718 : }
1719 : }
1720 : #endif
1721 :
1722 : /* ==================================================================== */
1723 : /* Default data dictionary. Emitted verbatim into the imagine */
1724 : /* file. */
1725 : /* ==================================================================== */
1726 :
1727 : static const char *aszDefaultDD[] = {
1728 : "{1:lversion,1:LfreeList,1:LrootEntryPtr,1:sentryHeaderLength,1:LdictionaryPtr,}Ehfa_File,{1:Lnext,1:Lprev,1:Lparent,1:Lchild,1:Ldata,1:ldataSize,64:cname,32:ctype,1:tmodTime,}Ehfa_Entry,{16:clabel,1:LheaderPtr,}Ehfa_HeaderTag,{1:LfreeList,1:lfreeSize,}Ehfa_FreeListNode,{1:lsize,1:Lptr,}Ehfa_Data,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,",
1729 : "1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer,{1:lwidth,1:lheight,1:e3:thematic,athematic,fft of real-valued data,layerType,1:e13:u1,u2,u4,u8,s8,u16,s16,u32,s32,f32,f64,c64,c128,pixelType,1:lblockWidth,1:lblockHeight,}Eimg_Layer_SubSample,{1:e2:raster,vector,type,1:LdictionaryPtr,}Ehfa_Layer,{1:LspaceUsedForRasterData,}ImgFormatInfo831,{1:sfileCode,1:Loffset,1:lsize,1:e2:false,true,logvalid,",
1730 : "1:e2:no compression,ESRI GRID compression,compressionType,}Edms_VirtualBlockInfo,{1:lmin,1:lmax,}Edms_FreeIDList,{1:lnumvirtualblocks,1:lnumobjectsperblock,1:lnextobjectnum,1:e2:no compression,RLC compression,compressionType,0:poEdms_VirtualBlockInfo,blockinfo,0:poEdms_FreeIDList,freelist,1:tmodTime,}Edms_State,{0:pcstring,}Emif_String,{1:oEmif_String,fileName,2:LlayerStackValidFlagsOffset,2:LlayerStackDataOffset,1:LlayerStackCount,1:LlayerStackIndex,}ImgExternalRaster,{1:oEmif_String,algorithm,0:poEmif_String,nameList,}Eimg_RRDNamesList,{1:oEmif_String,projection,1:oEmif_String,units,}Eimg_MapInformation,",
1731 : "{1:oEmif_String,dependent,}Eimg_DependentFile,{1:oEmif_String,ImageLayerName,}Eimg_DependentLayerName,{1:lnumrows,1:lnumcolumns,1:e13:EGDA_TYPE_U1,EGDA_TYPE_U2,EGDA_TYPE_U4,EGDA_TYPE_U8,EGDA_TYPE_S8,EGDA_TYPE_U16,EGDA_TYPE_S16,EGDA_TYPE_U32,EGDA_TYPE_S32,EGDA_TYPE_F32,EGDA_TYPE_F64,EGDA_TYPE_C64,EGDA_TYPE_C128,datatype,1:e4:EGDA_SCALAR_OBJECT,EGDA_TABLE_OBJECT,EGDA_MATRIX_OBJECT,EGDA_RASTER_OBJECT,objecttype,}Egda_BaseData,{1:*bvalueBD,}Eimg_NonInitializedValue,{1:dx,1:dy,}Eprj_Coordinate,{1:dwidth,1:dheight,}Eprj_Size,{0:pcproName,1:*oEprj_Coordinate,upperLeftCenter,",
1732 : "1:*oEprj_Coordinate,lowerRightCenter,1:*oEprj_Size,pixelSize,0:pcunits,}Eprj_MapInfo,{0:pcdatumname,1:e3:EPRJ_DATUM_PARAMETRIC,EPRJ_DATUM_GRID,EPRJ_DATUM_REGRESSION,type,0:pdparams,0:pcgridname,}Eprj_Datum,{0:pcsphereName,1:da,1:db,1:deSquared,1:dradius,}Eprj_Spheroid,{1:e2:EPRJ_INTERNAL,EPRJ_EXTERNAL,proType,1:lproNumber,0:pcproExeName,0:pcproName,1:lproZone,0:pdproParams,1:*oEprj_Spheroid,proSpheroid,}Eprj_ProParameters,{1:dminimum,1:dmaximum,1:dmean,1:dmedian,1:dmode,1:dstddev,}Esta_Statistics,{1:lnumBins,1:e4:direct,linear,logarithmic,explicit,binFunctionType,1:dminLimit,1:dmaxLimit,1:*bbinLimits,}Edsc_BinFunction,{0:poEmif_String,LayerNames,1:*bExcludedValues,1:oEmif_String,AOIname,",
1733 : "1:lSkipFactorX,1:lSkipFactorY,1:*oEdsc_BinFunction,BinFunction,}Eimg_StatisticsParameters830,{1:lnumrows,}Edsc_Table,{1:lnumRows,1:LcolumnDataPtr,1:e4:integer,real,complex,string,dataType,1:lmaxNumChars,}Edsc_Column,{1:lposition,0:pcname,1:e2:EMSC_FALSE,EMSC_TRUE,editable,1:e3:LEFT,CENTER,RIGHT,alignment,0:pcformat,1:e3:DEFAULT,APPLY,AUTO-APPLY,formulamode,0:pcformula,1:dcolumnwidth,0:pcunits,1:e5:NO_COLOR,RED,GREEN,BLUE,COLOR,colorflag,0:pcgreenname,0:pcbluename,}Eded_ColumnAttributes_1,{1:lversion,1:lnumobjects,1:e2:EAOI_UNION,EAOI_INTERSECTION,operation,}Eaoi_AreaOfInterest,",
1734 : "{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,",
1735 : "{1:x{1:x{0:pcstring,}Emif_String,type,1:x{0:pcstring,}Emif_String,MIFDictionary,0:pCMIFObject,}Emif_MIFObject,projection,1:x{0:pcstring,}Emif_String,title,}Eprj_MapProjection842,",
1736 : "{0:poEmif_String,titleList,}Exfr_GenericXFormHeader,{1:lorder,1:lnumdimtransform,1:lnumdimpolynomial,1:ltermcount,0:plexponentlist,1:*bpolycoefmtx,1:*bpolycoefvector,}Efga_Polynomial,",
1737 : ".",
1738 : NULL
1739 : };
1740 :
1741 :
1742 :
1743 : /************************************************************************/
1744 : /* HFACreateLL() */
1745 : /* */
1746 : /* Low level creation of an Imagine file. Writes out the */
1747 : /* Ehfa_HeaderTag, dictionary and Ehfa_File. */
1748 : /************************************************************************/
1749 :
1750 96 : HFAHandle HFACreateLL( const char * pszFilename )
1751 :
1752 : {
1753 : FILE *fp;
1754 : HFAInfo_t *psInfo;
1755 :
1756 : /* -------------------------------------------------------------------- */
1757 : /* Create the file in the file system. */
1758 : /* -------------------------------------------------------------------- */
1759 96 : fp = VSIFOpenL( pszFilename, "w+b" );
1760 96 : if( fp == NULL )
1761 : {
1762 : CPLError( CE_Failure, CPLE_OpenFailed,
1763 : "Creation of file %s failed.",
1764 0 : pszFilename );
1765 0 : return NULL;
1766 : }
1767 :
1768 : /* -------------------------------------------------------------------- */
1769 : /* Create the HFAInfo_t */
1770 : /* -------------------------------------------------------------------- */
1771 96 : psInfo = (HFAInfo_t *) CPLCalloc(sizeof(HFAInfo_t),1);
1772 :
1773 96 : psInfo->fp = fp;
1774 96 : psInfo->eAccess = HFA_Update;
1775 96 : psInfo->nXSize = 0;
1776 96 : psInfo->nYSize = 0;
1777 96 : psInfo->nBands = 0;
1778 96 : psInfo->papoBand = NULL;
1779 96 : psInfo->pMapInfo = NULL;
1780 96 : psInfo->pDatum = NULL;
1781 96 : psInfo->pProParameters = NULL;
1782 96 : psInfo->bTreeDirty = FALSE;
1783 96 : psInfo->pszFilename = CPLStrdup(CPLGetFilename(pszFilename));
1784 96 : psInfo->pszPath = CPLStrdup(CPLGetPath(pszFilename));
1785 :
1786 : /* -------------------------------------------------------------------- */
1787 : /* Write out the Ehfa_HeaderTag */
1788 : /* -------------------------------------------------------------------- */
1789 : GInt32 nHeaderPos;
1790 :
1791 96 : VSIFWriteL( (void *) "EHFA_HEADER_TAG", 1, 16, fp );
1792 :
1793 96 : nHeaderPos = 20;
1794 : HFAStandard( 4, &nHeaderPos );
1795 96 : VSIFWriteL( &nHeaderPos, 4, 1, fp );
1796 :
1797 : /* -------------------------------------------------------------------- */
1798 : /* Write the Ehfa_File node, locked in at offset 20. */
1799 : /* -------------------------------------------------------------------- */
1800 96 : GInt32 nVersion = 1, nFreeList = 0, nRootEntry = 0;
1801 96 : GInt16 nEntryHeaderLength = 128;
1802 96 : GInt32 nDictionaryPtr = 38;
1803 :
1804 96 : psInfo->nEntryHeaderLength = nEntryHeaderLength;
1805 96 : psInfo->nRootPos = 0;
1806 96 : psInfo->nDictionaryPos = nDictionaryPtr;
1807 96 : psInfo->nVersion = nVersion;
1808 :
1809 : HFAStandard( 4, &nVersion );
1810 : HFAStandard( 4, &nFreeList );
1811 : HFAStandard( 4, &nRootEntry );
1812 : HFAStandard( 2, &nEntryHeaderLength );
1813 : HFAStandard( 4, &nDictionaryPtr );
1814 :
1815 96 : VSIFWriteL( &nVersion, 4, 1, fp );
1816 96 : VSIFWriteL( &nFreeList, 4, 1, fp );
1817 96 : VSIFWriteL( &nRootEntry, 4, 1, fp );
1818 96 : VSIFWriteL( &nEntryHeaderLength, 2, 1, fp );
1819 96 : VSIFWriteL( &nDictionaryPtr, 4, 1, fp );
1820 :
1821 : /* -------------------------------------------------------------------- */
1822 : /* Write the dictionary, locked in at location 38. Note that */
1823 : /* we jump through a bunch of hoops to operate on the */
1824 : /* dictionary in chunks because some compiles (such as VC++) */
1825 : /* don't allow particularly large static strings. */
1826 : /* -------------------------------------------------------------------- */
1827 96 : int nDictLen = 0, iChunk;
1828 :
1829 1056 : for( iChunk = 0; aszDefaultDD[iChunk] != NULL; iChunk++ )
1830 960 : nDictLen += strlen(aszDefaultDD[iChunk]);
1831 :
1832 96 : psInfo->pszDictionary = (char *) CPLMalloc(nDictLen+1);
1833 96 : psInfo->pszDictionary[0] = '\0';
1834 :
1835 1056 : for( iChunk = 0; aszDefaultDD[iChunk] != NULL; iChunk++ )
1836 960 : strcat( psInfo->pszDictionary, aszDefaultDD[iChunk] );
1837 :
1838 : VSIFWriteL( (void *) psInfo->pszDictionary, 1,
1839 96 : strlen(psInfo->pszDictionary)+1, fp );
1840 :
1841 96 : psInfo->poDictionary = new HFADictionary( psInfo->pszDictionary );
1842 :
1843 96 : psInfo->nEndOfFile = (GUInt32) VSIFTellL( fp );
1844 :
1845 : /* -------------------------------------------------------------------- */
1846 : /* Create a root entry. */
1847 : /* -------------------------------------------------------------------- */
1848 192 : psInfo->poRoot = new HFAEntry( psInfo, "root", "root", NULL );
1849 :
1850 : /* -------------------------------------------------------------------- */
1851 : /* If an .ige or .rrd file exists with the same base name, */
1852 : /* delete them. (#1784) */
1853 : /* -------------------------------------------------------------------- */
1854 96 : CPLString osExtension = CPLGetExtension(pszFilename);
1855 192 : if( !EQUAL(osExtension,"rrd") && !EQUAL(osExtension,"aux") )
1856 : {
1857 93 : CPLString osPath = CPLGetPath( pszFilename );
1858 93 : CPLString osBasename = CPLGetBasename( pszFilename );
1859 : VSIStatBufL sStatBuf;
1860 93 : CPLString osSupFile = CPLFormCIFilename( osPath, osBasename, "rrd" );
1861 :
1862 93 : if( VSIStatL( osSupFile, &sStatBuf ) == 0 )
1863 0 : VSIUnlink( osSupFile );
1864 :
1865 93 : osSupFile = CPLFormCIFilename( osPath, osBasename, "ige" );
1866 :
1867 93 : if( VSIStatL( osSupFile, &sStatBuf ) == 0 )
1868 7 : VSIUnlink( osSupFile );
1869 : }
1870 :
1871 96 : return psInfo;
1872 : }
1873 :
1874 : /************************************************************************/
1875 : /* HFAAllocateSpace() */
1876 : /* */
1877 : /* Return an area in the file to the caller to write the */
1878 : /* requested number of bytes. Currently this is always at the */
1879 : /* end of the file, but eventually we might actually keep track */
1880 : /* of free space. The HFAInfo_t's concept of file size is */
1881 : /* updated, even if nothing ever gets written to this region. */
1882 : /* */
1883 : /* Returns the offset to the requested space, or zero one */
1884 : /* failure. */
1885 : /************************************************************************/
1886 :
1887 1490 : GUInt32 HFAAllocateSpace( HFAInfo_t *psInfo, GUInt32 nBytes )
1888 :
1889 : {
1890 : /* should check if this will wrap over 2GB limit */
1891 :
1892 1490 : psInfo->nEndOfFile += nBytes;
1893 1490 : return psInfo->nEndOfFile - nBytes;
1894 : }
1895 :
1896 : /************************************************************************/
1897 : /* HFAFlush() */
1898 : /* */
1899 : /* Write out any dirty tree information to disk, putting the */
1900 : /* disk file in a consistent state. */
1901 : /************************************************************************/
1902 :
1903 179 : CPLErr HFAFlush( HFAHandle hHFA )
1904 :
1905 : {
1906 : CPLErr eErr;
1907 :
1908 179 : if( !hHFA->bTreeDirty && !hHFA->poDictionary->bDictionaryTextDirty )
1909 0 : return CE_None;
1910 :
1911 : CPLAssert( hHFA->poRoot != NULL );
1912 :
1913 : /* -------------------------------------------------------------------- */
1914 : /* Flush HFAEntry tree to disk. */
1915 : /* -------------------------------------------------------------------- */
1916 179 : if( hHFA->bTreeDirty )
1917 : {
1918 178 : eErr = hHFA->poRoot->FlushToDisk();
1919 178 : if( eErr != CE_None )
1920 0 : return eErr;
1921 :
1922 178 : hHFA->bTreeDirty = FALSE;
1923 : }
1924 :
1925 : /* -------------------------------------------------------------------- */
1926 : /* Flush Dictionary to disk. */
1927 : /* -------------------------------------------------------------------- */
1928 179 : GUInt32 nNewDictionaryPos = hHFA->nDictionaryPos;
1929 :
1930 179 : if( hHFA->poDictionary->bDictionaryTextDirty )
1931 : {
1932 1 : VSIFSeekL( hHFA->fp, 0, SEEK_END );
1933 1 : nNewDictionaryPos = (GUInt32) VSIFTellL( hHFA->fp );
1934 : VSIFWriteL( hHFA->poDictionary->osDictionaryText.c_str(),
1935 : strlen(hHFA->poDictionary->osDictionaryText.c_str()) + 1,
1936 1 : 1, hHFA->fp );
1937 1 : hHFA->poDictionary->bDictionaryTextDirty = FALSE;
1938 : }
1939 :
1940 : /* -------------------------------------------------------------------- */
1941 : /* do we need to update the Ehfa_File pointer to the root node? */
1942 : /* -------------------------------------------------------------------- */
1943 179 : if( hHFA->nRootPos != hHFA->poRoot->GetFilePos()
1944 : || nNewDictionaryPos != hHFA->nDictionaryPos )
1945 : {
1946 : GUInt32 nOffset;
1947 : GUInt32 nHeaderPos;
1948 :
1949 97 : VSIFSeekL( hHFA->fp, 16, SEEK_SET );
1950 97 : VSIFReadL( &nHeaderPos, sizeof(GInt32), 1, hHFA->fp );
1951 : HFAStandard( 4, &nHeaderPos );
1952 :
1953 97 : nOffset = hHFA->nRootPos = hHFA->poRoot->GetFilePos();
1954 : HFAStandard( 4, &nOffset );
1955 97 : VSIFSeekL( hHFA->fp, nHeaderPos+8, SEEK_SET );
1956 97 : VSIFWriteL( &nOffset, 4, 1, hHFA->fp );
1957 :
1958 97 : nOffset = hHFA->nDictionaryPos = nNewDictionaryPos;
1959 : HFAStandard( 4, &nOffset );
1960 97 : VSIFSeekL( hHFA->fp, nHeaderPos+14, SEEK_SET );
1961 97 : VSIFWriteL( &nOffset, 4, 1, hHFA->fp );
1962 : }
1963 :
1964 179 : return CE_None;
1965 : }
1966 :
1967 : /************************************************************************/
1968 : /* HFACreateLayer() */
1969 : /* */
1970 : /* Create a layer object, and corresponding RasterDMS. */
1971 : /* Suitable for use with primary layers, and overviews. */
1972 : /************************************************************************/
1973 :
1974 : int
1975 150 : HFACreateLayer( HFAHandle psInfo, HFAEntry *poParent,
1976 : const char *pszLayerName,
1977 : int bOverview, int nBlockSize,
1978 : int bCreateCompressed, int bCreateLargeRaster,
1979 : int bDependentLayer,
1980 : int nXSize, int nYSize, int nDataType,
1981 : char **papszOptions,
1982 :
1983 : // these are only related to external (large) files
1984 : GIntBig nStackValidFlagsOffset,
1985 : GIntBig nStackDataOffset,
1986 : int nStackCount, int nStackIndex )
1987 :
1988 : {
1989 :
1990 : HFAEntry *poEimg_Layer;
1991 : const char *pszLayerType;
1992 :
1993 150 : if( bOverview )
1994 5 : pszLayerType = "Eimg_Layer_SubSample";
1995 : else
1996 145 : pszLayerType = "Eimg_Layer";
1997 :
1998 150 : if (nBlockSize <= 0)
1999 : {
2000 0 : CPLError(CE_Failure, CPLE_IllegalArg, "HFACreateLayer : nBlockXSize < 0");
2001 0 : return FALSE;
2002 : }
2003 :
2004 : /* -------------------------------------------------------------------- */
2005 : /* Work out some details about the tiling scheme. */
2006 : /* -------------------------------------------------------------------- */
2007 : int nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock;
2008 :
2009 150 : nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2010 150 : nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2011 150 : nBlocks = nBlocksPerRow * nBlocksPerColumn;
2012 : nBytesPerBlock = (nBlockSize * nBlockSize
2013 150 : * HFAGetDataTypeBits(nDataType) + 7) / 8;
2014 :
2015 : /* -------------------------------------------------------------------- */
2016 : /* Create the Eimg_Layer for the band. */
2017 : /* -------------------------------------------------------------------- */
2018 : poEimg_Layer =
2019 150 : new HFAEntry( psInfo, pszLayerName, pszLayerType, poParent );
2020 :
2021 150 : poEimg_Layer->SetIntField( "width", nXSize );
2022 150 : poEimg_Layer->SetIntField( "height", nYSize );
2023 150 : poEimg_Layer->SetStringField( "layerType", "athematic" );
2024 150 : poEimg_Layer->SetIntField( "pixelType", nDataType );
2025 150 : poEimg_Layer->SetIntField( "blockWidth", nBlockSize );
2026 150 : poEimg_Layer->SetIntField( "blockHeight", nBlockSize );
2027 :
2028 : /* -------------------------------------------------------------------- */
2029 : /* Create the RasterDMS (block list). This is a complex type */
2030 : /* with pointers, and variable size. We set the superstructure */
2031 : /* ourselves rather than trying to have the HFA type management */
2032 : /* system do it for us (since this would be hard to implement). */
2033 : /* -------------------------------------------------------------------- */
2034 424 : if ( !bCreateLargeRaster && !bDependentLayer )
2035 : {
2036 : int nDmsSize;
2037 : HFAEntry *poEdms_State;
2038 : GByte *pabyData;
2039 :
2040 : poEdms_State =
2041 137 : new HFAEntry( psInfo, "RasterDMS", "Edms_State", poEimg_Layer );
2042 :
2043 137 : nDmsSize = 14 * nBlocks + 38;
2044 137 : pabyData = poEdms_State->MakeData( nDmsSize );
2045 :
2046 : /* set some simple values */
2047 137 : poEdms_State->SetIntField( "numvirtualblocks", nBlocks );
2048 : poEdms_State->SetIntField( "numobjectsperblock",
2049 137 : nBlockSize*nBlockSize );
2050 : poEdms_State->SetIntField( "nextobjectnum",
2051 137 : nBlockSize*nBlockSize*nBlocks );
2052 :
2053 : /* Is file compressed or not? */
2054 137 : if( bCreateCompressed )
2055 : {
2056 14 : poEdms_State->SetStringField( "compressionType", "RLC compression" );
2057 : }
2058 : else
2059 : {
2060 123 : poEdms_State->SetStringField( "compressionType", "no compression" );
2061 : }
2062 :
2063 : /* we need to hardcode file offset into the data, so locate it now */
2064 137 : poEdms_State->SetPosition();
2065 :
2066 : /* Set block info headers */
2067 : GUInt32 nValue;
2068 :
2069 : /* blockinfo count */
2070 137 : nValue = nBlocks;
2071 : HFAStandard( 4, &nValue );
2072 137 : memcpy( pabyData + 14, &nValue, 4 );
2073 :
2074 : /* blockinfo position */
2075 137 : nValue = poEdms_State->GetDataPos() + 22;
2076 : HFAStandard( 4, &nValue );
2077 137 : memcpy( pabyData + 18, &nValue, 4 );
2078 :
2079 : /* Set each blockinfo */
2080 438 : for( int iBlock = 0; iBlock < nBlocks; iBlock++ )
2081 : {
2082 : GInt16 nValue16;
2083 301 : int nOffset = 22 + 14 * iBlock;
2084 :
2085 : /* fileCode */
2086 301 : nValue16 = 0;
2087 : HFAStandard( 2, &nValue16 );
2088 301 : memcpy( pabyData + nOffset, &nValue16, 2 );
2089 :
2090 : /* offset */
2091 301 : if( bCreateCompressed )
2092 : {
2093 : /* flag it with zero offset - will allocate space when we compress it */
2094 14 : nValue = 0;
2095 : }
2096 : else
2097 : {
2098 287 : nValue = HFAAllocateSpace( psInfo, nBytesPerBlock );
2099 : }
2100 : HFAStandard( 4, &nValue );
2101 301 : memcpy( pabyData + nOffset + 2, &nValue, 4 );
2102 :
2103 : /* size */
2104 301 : if( bCreateCompressed )
2105 : {
2106 : /* flag it with zero size - don't know until we compress it */
2107 14 : nValue = 0;
2108 : }
2109 : else
2110 : {
2111 287 : nValue = nBytesPerBlock;
2112 : }
2113 : HFAStandard( 4, &nValue );
2114 301 : memcpy( pabyData + nOffset + 6, &nValue, 4 );
2115 :
2116 : /* logValid (false) */
2117 301 : nValue16 = 0;
2118 : HFAStandard( 2, &nValue16 );
2119 301 : memcpy( pabyData + nOffset + 10, &nValue16, 2 );
2120 :
2121 : /* compressionType */
2122 301 : if( bCreateCompressed )
2123 14 : nValue16 = 1;
2124 : else
2125 287 : nValue16 = 0;
2126 :
2127 : HFAStandard( 2, &nValue16 );
2128 301 : memcpy( pabyData + nOffset + 12, &nValue16, 2 );
2129 : }
2130 :
2131 : }
2132 : /* -------------------------------------------------------------------- */
2133 : /* Create ExternalRasterDMS object. */
2134 : /* -------------------------------------------------------------------- */
2135 13 : else if( bCreateLargeRaster )
2136 : {
2137 : HFAEntry *poEdms_State;
2138 :
2139 : poEdms_State =
2140 : new HFAEntry( psInfo, "ExternalRasterDMS",
2141 13 : "ImgExternalRaster", poEimg_Layer );
2142 13 : poEdms_State->MakeData( 8 + strlen(psInfo->pszIGEFilename) + 1 + 6 * 4 );
2143 :
2144 : poEdms_State->SetStringField( "fileName.string",
2145 13 : psInfo->pszIGEFilename );
2146 :
2147 : poEdms_State->SetIntField( "layerStackValidFlagsOffset[0]",
2148 13 : (int) (nStackValidFlagsOffset & 0xFFFFFFFF));
2149 : poEdms_State->SetIntField( "layerStackValidFlagsOffset[1]",
2150 13 : (int) (nStackValidFlagsOffset >> 32) );
2151 :
2152 : poEdms_State->SetIntField( "layerStackDataOffset[0]",
2153 13 : (int) (nStackDataOffset & 0xFFFFFFFF) );
2154 : poEdms_State->SetIntField( "layerStackDataOffset[1]",
2155 13 : (int) (nStackDataOffset >> 32 ) );
2156 13 : poEdms_State->SetIntField( "layerStackCount", nStackCount );
2157 13 : poEdms_State->SetIntField( "layerStackIndex", nStackIndex );
2158 : }
2159 :
2160 : /* -------------------------------------------------------------------- */
2161 : /* Dependent... */
2162 : /* -------------------------------------------------------------------- */
2163 0 : else if( bDependentLayer )
2164 : {
2165 : HFAEntry *poDepLayerName;
2166 :
2167 : poDepLayerName =
2168 : new HFAEntry( psInfo, "DependentLayerName",
2169 0 : "Eimg_DependentLayerName", poEimg_Layer );
2170 0 : poDepLayerName->MakeData( 8 + strlen(pszLayerName) + 2 );
2171 :
2172 : poDepLayerName->SetStringField( "ImageLayerName.string",
2173 0 : pszLayerName );
2174 : }
2175 :
2176 : /* -------------------------------------------------------------------- */
2177 : /* Create the Ehfa_Layer. */
2178 : /* -------------------------------------------------------------------- */
2179 : HFAEntry *poEhfa_Layer;
2180 : GUInt32 nLDict;
2181 : char szLDict[128], chBandType;
2182 :
2183 150 : if( nDataType == EPT_u1 )
2184 2 : chBandType = '1';
2185 148 : else if( nDataType == EPT_u2 )
2186 0 : chBandType = '2';
2187 148 : else if( nDataType == EPT_u4 )
2188 0 : chBandType = '4';
2189 148 : else if( nDataType == EPT_u8 )
2190 69 : chBandType = 'c';
2191 79 : else if( nDataType == EPT_s8 )
2192 0 : chBandType = 'C';
2193 79 : else if( nDataType == EPT_u16 )
2194 16 : chBandType = 's';
2195 63 : else if( nDataType == EPT_s16 )
2196 8 : chBandType = 'S';
2197 55 : else if( nDataType == EPT_u32 )
2198 : // for some reason erdas imagine expects an L for unsinged 32 bit ints
2199 : // otherwise it gives strange "out of memory errors"
2200 8 : chBandType = 'L';
2201 47 : else if( nDataType == EPT_s32 )
2202 8 : chBandType = 'L';
2203 39 : else if( nDataType == EPT_f32 )
2204 8 : chBandType = 'f';
2205 31 : else if( nDataType == EPT_f64 )
2206 15 : chBandType = 'd';
2207 16 : else if( nDataType == EPT_c64 )
2208 8 : chBandType = 'm';
2209 8 : else if( nDataType == EPT_c128 )
2210 8 : chBandType = 'M';
2211 : else
2212 : {
2213 : CPLAssert( FALSE );
2214 0 : chBandType = 'c';
2215 : }
2216 :
2217 : // the first value in the entry below gives the number of pixels within a block
2218 150 : sprintf( szLDict, "{%d:%cdata,}RasterDMS,.", nBlockSize*nBlockSize, chBandType );
2219 :
2220 : poEhfa_Layer = new HFAEntry( psInfo, "Ehfa_Layer", "Ehfa_Layer",
2221 150 : poEimg_Layer );
2222 150 : poEhfa_Layer->MakeData();
2223 150 : poEhfa_Layer->SetPosition();
2224 150 : nLDict = HFAAllocateSpace( psInfo, strlen(szLDict) + 1 );
2225 :
2226 150 : poEhfa_Layer->SetStringField( "type", "raster" );
2227 150 : poEhfa_Layer->SetIntField( "dictionaryPtr", nLDict );
2228 :
2229 150 : VSIFSeekL( psInfo->fp, nLDict, SEEK_SET );
2230 150 : VSIFWriteL( (void *) szLDict, strlen(szLDict) + 1, 1, psInfo->fp );
2231 :
2232 150 : return TRUE;
2233 : }
2234 :
2235 :
2236 : /************************************************************************/
2237 : /* HFACreate() */
2238 : /************************************************************************/
2239 :
2240 95 : HFAHandle HFACreate( const char * pszFilename,
2241 : int nXSize, int nYSize, int nBands,
2242 : int nDataType, char ** papszOptions )
2243 :
2244 : {
2245 : HFAHandle psInfo;
2246 95 : int nBlockSize = 64;
2247 95 : const char * pszValue = CSLFetchNameValue( papszOptions, "BLOCKSIZE" );
2248 :
2249 95 : if ( pszValue != NULL )
2250 : {
2251 0 : nBlockSize = atoi( pszValue );
2252 : // check for sane values
2253 0 : if ( ( nBlockSize < 32 ) || (nBlockSize > 2048) )
2254 : {
2255 0 : nBlockSize = 64;
2256 : }
2257 : }
2258 : int bCreateLargeRaster = CSLFetchBoolean(papszOptions,"USE_SPILL",
2259 95 : FALSE);
2260 : int bCreateCompressed =
2261 : CSLFetchBoolean(papszOptions,"COMPRESS", FALSE)
2262 95 : || CSLFetchBoolean(papszOptions,"COMPRESSED", FALSE);
2263 95 : int bCreateAux = CSLFetchBoolean(papszOptions,"AUX", FALSE);
2264 :
2265 95 : char *pszFullFilename = NULL, *pszRawFilename = NULL;
2266 :
2267 : /* -------------------------------------------------------------------- */
2268 : /* Create the low level structure. */
2269 : /* -------------------------------------------------------------------- */
2270 95 : psInfo = HFACreateLL( pszFilename );
2271 95 : if( psInfo == NULL )
2272 0 : return NULL;
2273 :
2274 : /* -------------------------------------------------------------------- */
2275 : /* Create the DependentFile node if requested. */
2276 : /* -------------------------------------------------------------------- */
2277 : const char *pszDependentFile =
2278 95 : CSLFetchNameValue( papszOptions, "DEPENDENT_FILE" );
2279 :
2280 95 : if( pszDependentFile != NULL )
2281 : {
2282 : HFAEntry *poDF = new HFAEntry( psInfo, "DependentFile",
2283 2 : "Eimg_DependentFile", psInfo->poRoot );
2284 :
2285 2 : poDF->MakeData( strlen(pszDependentFile) + 50 );
2286 2 : poDF->SetPosition();
2287 2 : poDF->SetStringField( "dependent.string", pszDependentFile );
2288 : }
2289 :
2290 : /* -------------------------------------------------------------------- */
2291 : /* Work out some details about the tiling scheme. */
2292 : /* -------------------------------------------------------------------- */
2293 : int nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock;
2294 :
2295 95 : nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2296 95 : nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2297 95 : nBlocks = nBlocksPerRow * nBlocksPerColumn;
2298 : nBytesPerBlock = (nBlockSize * nBlockSize
2299 95 : * HFAGetDataTypeBits(nDataType) + 7) / 8;
2300 :
2301 : CPLDebug( "HFACreate", "Blocks per row %d, blocks per column %d, "
2302 : "total number of blocks %d, bytes per block %d.",
2303 95 : nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock );
2304 :
2305 : /* -------------------------------------------------------------------- */
2306 : /* Check whether we should create external large file with */
2307 : /* image. We create a spill file if the amount of imagery is */
2308 : /* close to 2GB. We don't check the amount of auxilary */
2309 : /* information, so in theory if there were an awful lot of */
2310 : /* non-imagery data our approximate size could be smaller than */
2311 : /* the file will actually we be. We leave room for 10MB of */
2312 : /* auxilary data. */
2313 : /* We can also force spill file creation using option */
2314 : /* SPILL_FILE=YES. */
2315 : /* -------------------------------------------------------------------- */
2316 : double dfApproxSize = (double)nBytesPerBlock * (double)nBlocks *
2317 95 : (double)nBands + 10000000.0;
2318 :
2319 95 : if( dfApproxSize > 2147483648.0 && !bCreateAux )
2320 0 : bCreateLargeRaster = TRUE;
2321 :
2322 : // erdas imagine creates this entry even if an external spill file is used
2323 95 : if( !bCreateAux )
2324 : {
2325 : HFAEntry *poImgFormat;
2326 : poImgFormat = new HFAEntry( psInfo, "IMGFormatInfo",
2327 95 : "ImgFormatInfo831", psInfo->poRoot );
2328 95 : poImgFormat->MakeData();
2329 95 : if ( bCreateLargeRaster )
2330 : {
2331 7 : poImgFormat->SetIntField( "spaceUsedForRasterData", 0 );
2332 7 : bCreateCompressed = FALSE; // Can't be compressed if we are creating a spillfile
2333 : }
2334 : else
2335 : {
2336 : poImgFormat->SetIntField( "spaceUsedForRasterData",
2337 88 : nBytesPerBlock*nBlocks*nBands );
2338 : }
2339 : }
2340 :
2341 : /* -------------------------------------------------------------------- */
2342 : /* Create external file and write its header. */
2343 : /* -------------------------------------------------------------------- */
2344 95 : GIntBig nValidFlagsOffset = 0, nDataOffset = 0;
2345 :
2346 95 : if( bCreateLargeRaster )
2347 : {
2348 7 : if( !HFACreateSpillStack( psInfo, nXSize, nYSize, nBands,
2349 : nBlockSize, nDataType,
2350 : &nValidFlagsOffset, &nDataOffset ) )
2351 : {
2352 0 : CPLFree( pszRawFilename );
2353 0 : CPLFree( pszFullFilename );
2354 0 : return NULL;
2355 : }
2356 : }
2357 :
2358 : /* ==================================================================== */
2359 : /* Create each band (layer) */
2360 : /* ==================================================================== */
2361 : int iBand;
2362 :
2363 240 : for( iBand = 0; iBand < nBands; iBand++ )
2364 : {
2365 : char szName[128];
2366 :
2367 145 : sprintf( szName, "Layer_%d", iBand + 1 );
2368 :
2369 145 : if( !HFACreateLayer( psInfo, psInfo->poRoot, szName, FALSE, nBlockSize,
2370 : bCreateCompressed, bCreateLargeRaster, bCreateAux,
2371 : nXSize, nYSize, nDataType, papszOptions,
2372 : nValidFlagsOffset, nDataOffset,
2373 : nBands, iBand ) )
2374 : {
2375 0 : HFAClose( psInfo );
2376 0 : return NULL;
2377 : }
2378 : }
2379 :
2380 : /* -------------------------------------------------------------------- */
2381 : /* Initialize the band information. */
2382 : /* -------------------------------------------------------------------- */
2383 95 : HFAParseBandInfo( psInfo );
2384 :
2385 95 : return psInfo;
2386 : }
2387 :
2388 : /************************************************************************/
2389 : /* HFACreateOverview() */
2390 : /* */
2391 : /* Create an overview layer object for a band. */
2392 : /************************************************************************/
2393 :
2394 5 : int HFACreateOverview( HFAHandle hHFA, int nBand, int nOverviewLevel )
2395 :
2396 : {
2397 5 : if( nBand < 1 || nBand > hHFA->nBands )
2398 0 : return -1;
2399 : else
2400 : {
2401 5 : HFABand *poBand = hHFA->papoBand[nBand-1];
2402 5 : return poBand->CreateOverview( nOverviewLevel );
2403 : }
2404 : }
2405 :
2406 : /************************************************************************/
2407 : /* HFAGetMetadata() */
2408 : /* */
2409 : /* Read metadata structured in a table called GDAL_MetaData. */
2410 : /************************************************************************/
2411 :
2412 678 : char ** HFAGetMetadata( HFAHandle hHFA, int nBand )
2413 :
2414 : {
2415 : HFAEntry *poTable;
2416 :
2417 1081 : if( nBand > 0 && nBand <= hHFA->nBands )
2418 403 : poTable = hHFA->papoBand[nBand - 1]->poNode->GetChild();
2419 275 : else if( nBand == 0 )
2420 275 : poTable = hHFA->poRoot->GetChild();
2421 : else
2422 0 : return NULL;
2423 :
2424 678 : for( ; poTable != NULL && !EQUAL(poTable->GetName(),"GDAL_MetaData");
2425 : poTable = poTable->GetNext() ) {}
2426 :
2427 678 : if( poTable == NULL || !EQUAL(poTable->GetType(),"Edsc_Table") )
2428 608 : return NULL;
2429 :
2430 70 : if( poTable->GetIntField( "numRows" ) != 1 )
2431 : {
2432 : CPLDebug( "HFADataset", "GDAL_MetaData.numRows = %d, expected 1!",
2433 0 : poTable->GetIntField( "numRows" ) );
2434 0 : return NULL;
2435 : }
2436 :
2437 : /* -------------------------------------------------------------------- */
2438 : /* Loop over each column. Each column will be one metadata */
2439 : /* entry, with the title being the key, and the row value being */
2440 : /* the value. There is only ever one row in GDAL_MetaData */
2441 : /* tables. */
2442 : /* -------------------------------------------------------------------- */
2443 : HFAEntry *poColumn;
2444 70 : char **papszMD = NULL;
2445 :
2446 212 : for( poColumn = poTable->GetChild();
2447 : poColumn != NULL;
2448 : poColumn = poColumn->GetNext() )
2449 : {
2450 : const char *pszValue;
2451 : int columnDataPtr;
2452 :
2453 : // Skip the #Bin_Function# entry.
2454 142 : if( EQUALN(poColumn->GetName(),"#",1) )
2455 70 : continue;
2456 :
2457 72 : pszValue = poColumn->GetStringField( "dataType" );
2458 72 : if( pszValue == NULL || !EQUAL(pszValue,"string") )
2459 0 : continue;
2460 :
2461 72 : columnDataPtr = poColumn->GetIntField( "columnDataPtr" );
2462 72 : if( columnDataPtr == 0 )
2463 0 : continue;
2464 :
2465 : /* -------------------------------------------------------------------- */
2466 : /* read up to nMaxNumChars bytes from the indicated location. */
2467 : /* allocate required space temporarily */
2468 : /* nMaxNumChars should have been set by GDAL orginally so we should*/
2469 : /* trust it, but who knows... */
2470 : /* -------------------------------------------------------------------- */
2471 72 : int nMaxNumChars = poColumn->GetIntField( "maxNumChars" );
2472 :
2473 72 : if( nMaxNumChars == 0 )
2474 : {
2475 0 : papszMD = CSLSetNameValue( papszMD, poColumn->GetName(), "" );
2476 : }
2477 : else
2478 : {
2479 72 : char *pszMDValue = (char*) VSIMalloc(nMaxNumChars);
2480 72 : if (pszMDValue == NULL)
2481 : {
2482 : CPLError(CE_Failure, CPLE_OutOfMemory,
2483 0 : "HFAGetMetadata : Out of memory while allocating %d bytes", nMaxNumChars);
2484 0 : continue;
2485 : }
2486 :
2487 72 : if( VSIFSeekL( hHFA->fp, columnDataPtr, SEEK_SET ) != 0 )
2488 0 : continue;
2489 :
2490 72 : int nMDBytes = VSIFReadL( pszMDValue, 1, nMaxNumChars, hHFA->fp );
2491 72 : if( nMDBytes == 0 )
2492 : {
2493 0 : CPLFree( pszMDValue );
2494 0 : continue;
2495 : }
2496 :
2497 72 : pszMDValue[nMaxNumChars-1] = '\0';
2498 :
2499 : papszMD = CSLSetNameValue( papszMD, poColumn->GetName(),
2500 72 : pszMDValue );
2501 72 : CPLFree( pszMDValue );
2502 : }
2503 : }
2504 :
2505 70 : return papszMD;
2506 : }
2507 :
2508 : /************************************************************************/
2509 : /* HFASetGDALMetadata() */
2510 : /* */
2511 : /* This function is used to set metadata in a table called */
2512 : /* GDAL_MetaData. It is called by HFASetMetadata() for all */
2513 : /* metadata items that aren't some specific supported */
2514 : /* information (like histogram or stats info). */
2515 : /************************************************************************/
2516 :
2517 : static CPLErr
2518 44 : HFASetGDALMetadata( HFAHandle hHFA, int nBand, char **papszMD )
2519 :
2520 : {
2521 44 : if( papszMD == NULL )
2522 0 : return CE_None;
2523 :
2524 : HFAEntry *poNode;
2525 :
2526 46 : if( nBand > 0 && nBand <= hHFA->nBands )
2527 2 : poNode = hHFA->papoBand[nBand - 1]->poNode;
2528 42 : else if( nBand == 0 )
2529 42 : poNode = hHFA->poRoot;
2530 : else
2531 0 : return CE_Failure;
2532 :
2533 : /* -------------------------------------------------------------------- */
2534 : /* Create the Descriptor table. */
2535 : /* Check we have no table with this name already */
2536 : /* -------------------------------------------------------------------- */
2537 44 : HFAEntry *poEdsc_Table = poNode->GetNamedChild( "GDAL_MetaData" );
2538 :
2539 44 : if( poEdsc_Table == NULL || !EQUAL(poEdsc_Table->GetType(),"Edsc_Table") )
2540 : poEdsc_Table = new HFAEntry( hHFA, "GDAL_MetaData", "Edsc_Table",
2541 43 : poNode );
2542 :
2543 44 : poEdsc_Table->SetIntField( "numrows", 1 );
2544 :
2545 : /* -------------------------------------------------------------------- */
2546 : /* Create the Binning function node. I am not sure that we */
2547 : /* really need this though. */
2548 : /* Check it doesn't exist already */
2549 : /* -------------------------------------------------------------------- */
2550 : HFAEntry *poEdsc_BinFunction =
2551 44 : poEdsc_Table->GetNamedChild( "#Bin_Function#" );
2552 :
2553 44 : if( poEdsc_BinFunction == NULL
2554 : || !EQUAL(poEdsc_BinFunction->GetType(),"Edsc_BinFunction") )
2555 : poEdsc_BinFunction = new HFAEntry( hHFA, "#Bin_Function#",
2556 43 : "Edsc_BinFunction", poEdsc_Table );
2557 :
2558 : // Because of the BaseData we have to hardcode the size.
2559 44 : poEdsc_BinFunction->MakeData( 30 );
2560 :
2561 44 : poEdsc_BinFunction->SetIntField( "numBins", 1 );
2562 44 : poEdsc_BinFunction->SetStringField( "binFunction", "direct" );
2563 44 : poEdsc_BinFunction->SetDoubleField( "minLimit", 0.0 );
2564 44 : poEdsc_BinFunction->SetDoubleField( "maxLimit", 0.0 );
2565 :
2566 : /* -------------------------------------------------------------------- */
2567 : /* Process each metadata item as a separate column. */
2568 : /* -------------------------------------------------------------------- */
2569 89 : for( int iColumn = 0; papszMD[iColumn] != NULL; iColumn++ )
2570 : {
2571 : HFAEntry *poEdsc_Column;
2572 45 : char *pszKey = NULL;
2573 : const char *pszValue;
2574 :
2575 45 : pszValue = CPLParseNameValue( papszMD[iColumn], &pszKey );
2576 45 : if( pszValue == NULL )
2577 0 : continue;
2578 :
2579 : /* -------------------------------------------------------------------- */
2580 : /* Create the Edsc_Column. */
2581 : /* Check it doesn't exist already */
2582 : /* -------------------------------------------------------------------- */
2583 45 : poEdsc_Column = poEdsc_Table->GetNamedChild(pszKey);
2584 :
2585 45 : if( poEdsc_Column == NULL
2586 : || !EQUAL(poEdsc_Column->GetType(),"Edsc_Column") )
2587 : poEdsc_Column = new HFAEntry( hHFA, pszKey, "Edsc_Column",
2588 44 : poEdsc_Table );
2589 :
2590 45 : poEdsc_Column->SetIntField( "numRows", 1 );
2591 45 : poEdsc_Column->SetStringField( "dataType", "string" );
2592 45 : poEdsc_Column->SetIntField( "maxNumChars", strlen(pszValue)+1 );
2593 :
2594 : /* -------------------------------------------------------------------- */
2595 : /* Write the data out. */
2596 : /* -------------------------------------------------------------------- */
2597 45 : int nOffset = HFAAllocateSpace( hHFA, strlen(pszValue)+1);
2598 :
2599 45 : poEdsc_Column->SetIntField( "columnDataPtr", nOffset );
2600 :
2601 45 : VSIFSeekL( hHFA->fp, nOffset, SEEK_SET );
2602 45 : VSIFWriteL( (void *) pszValue, 1, strlen(pszValue)+1, hHFA->fp );
2603 :
2604 45 : CPLFree( pszKey );
2605 : }
2606 :
2607 44 : return CE_Failure;
2608 : }
2609 :
2610 : /************************************************************************/
2611 : /* HFASetMetadata() */
2612 : /************************************************************************/
2613 :
2614 46 : CPLErr HFASetMetadata( HFAHandle hHFA, int nBand, char **papszMD )
2615 :
2616 : {
2617 46 : char **papszGDALMD = NULL;
2618 :
2619 46 : if( CSLCount(papszMD) == 0 )
2620 0 : return CE_None;
2621 :
2622 : HFAEntry *poNode;
2623 :
2624 50 : if( nBand > 0 && nBand <= hHFA->nBands )
2625 4 : poNode = hHFA->papoBand[nBand - 1]->poNode;
2626 42 : else if( nBand == 0 )
2627 42 : poNode = hHFA->poRoot;
2628 : else
2629 0 : return CE_Failure;
2630 :
2631 : /* -------------------------------------------------------------------- */
2632 : /* Check if the Metadata is an "known" entity which should be */
2633 : /* stored in a better place. */
2634 : /* -------------------------------------------------------------------- */
2635 46 : char * pszBinValues = NULL;
2636 46 : int bCreatedHistogramParameters = FALSE;
2637 46 : int bCreatedStatistics = FALSE;
2638 46 : const char ** pszAuxMetaData = GetHFAAuxMetaDataList();
2639 : // check each metadata item
2640 120 : for( int iColumn = 0; papszMD[iColumn] != NULL; iColumn++ )
2641 : {
2642 74 : char *pszKey = NULL;
2643 : const char *pszValue;
2644 :
2645 74 : pszValue = CPLParseNameValue( papszMD[iColumn], &pszKey );
2646 74 : if( pszValue == NULL )
2647 0 : continue;
2648 :
2649 : // know look if its known
2650 : int i;
2651 853 : for( i = 0; pszAuxMetaData[i] != NULL; i += 4 )
2652 : {
2653 806 : if ( EQUALN( pszAuxMetaData[i + 2], pszKey, strlen(pszKey) ) )
2654 27 : break;
2655 : }
2656 74 : if ( pszAuxMetaData[i] != NULL )
2657 : {
2658 : // found one, get the right entry
2659 : HFAEntry *poEntry;
2660 :
2661 27 : if( strlen(pszAuxMetaData[i]) > 0 )
2662 24 : poEntry = poNode->GetNamedChild( pszAuxMetaData[i] );
2663 : else
2664 3 : poEntry = poNode;
2665 :
2666 27 : if( poEntry == NULL && strlen(pszAuxMetaData[i+3]) > 0 )
2667 : {
2668 : // child does not yet exist --> create it
2669 : poEntry = new HFAEntry( hHFA, pszAuxMetaData[i], pszAuxMetaData[i+3],
2670 4 : poNode );
2671 :
2672 4 : if ( EQUALN( "Statistics", pszAuxMetaData[i], 10 ) )
2673 2 : bCreatedStatistics = TRUE;
2674 :
2675 4 : if( EQUALN( "HistogramParameters", pszAuxMetaData[i], 19 ) )
2676 : {
2677 : // this is a bit nasty I need to set the string field for the object
2678 : // first because the SetStringField sets the count for the object
2679 : // BinFunction to the length of the string
2680 2 : poEntry->MakeData( 70 );
2681 2 : poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
2682 :
2683 2 : bCreatedHistogramParameters = TRUE;
2684 : }
2685 : }
2686 27 : if ( poEntry == NULL )
2687 : {
2688 6 : CPLFree( pszKey );
2689 6 : continue;
2690 : }
2691 :
2692 21 : const char *pszFieldName = pszAuxMetaData[i+1] + 1;
2693 21 : switch( pszAuxMetaData[i+1][0] )
2694 : {
2695 : case 'd':
2696 : {
2697 16 : double dfValue = atof( pszValue );
2698 16 : poEntry->SetDoubleField( pszFieldName, dfValue );
2699 : }
2700 16 : break;
2701 : case 'i':
2702 : case 'l':
2703 : {
2704 2 : int nValue = atoi( pszValue );
2705 2 : poEntry->SetIntField( pszFieldName, nValue );
2706 : }
2707 2 : break;
2708 : case 's':
2709 : case 'e':
2710 : {
2711 3 : poEntry->SetStringField( pszFieldName, pszValue );
2712 : }
2713 : break;
2714 : default:
2715 : CPLAssert( FALSE );
2716 : }
2717 : }
2718 47 : else if ( EQUALN( "STATISTICS_HISTOBINVALUES", pszKey, strlen(pszKey) ) )
2719 : {
2720 2 : pszBinValues = strdup( pszValue );
2721 : }
2722 : else
2723 45 : papszGDALMD = CSLAddString( papszGDALMD, papszMD[iColumn] );
2724 :
2725 68 : CPLFree( pszKey );
2726 : }
2727 :
2728 : /* -------------------------------------------------------------------- */
2729 : /* Special case to write out the histogram. */
2730 : /* -------------------------------------------------------------------- */
2731 46 : if ( pszBinValues != NULL )
2732 : {
2733 2 : HFAEntry * poEntry = poNode->GetNamedChild( "HistogramParameters" );
2734 2 : if ( poEntry != NULL && bCreatedHistogramParameters )
2735 : {
2736 : // if this node exists we have added Histogram data -- complete with some defaults
2737 2 : poEntry->SetIntField( "SkipFactorX", 1 );
2738 2 : poEntry->SetIntField( "SkipFactorY", 1 );
2739 :
2740 2 : int nNumBins = poEntry->GetIntField( "BinFunction.numBins" );
2741 2 : double dMinLimit = poEntry->GetDoubleField( "BinFunction.minLimit" );
2742 2 : double dMaxLimit = poEntry->GetDoubleField( "BinFunction.maxLimit" );
2743 :
2744 : // fill the descriptor table - check it isn't there already
2745 2 : poEntry = poNode->GetNamedChild( "Descriptor_Table" );
2746 2 : if( poEntry == NULL || !EQUAL(poEntry->GetType(),"Edsc_Table") )
2747 2 : poEntry = new HFAEntry( hHFA, "Descriptor_Table", "Edsc_Table", poNode );
2748 :
2749 2 : poEntry->SetIntField( "numRows", nNumBins );
2750 :
2751 : // bin function
2752 2 : HFAEntry * poBinFunc = poEntry->GetNamedChild( "#Bin_Function#" );
2753 2 : if( poBinFunc == NULL || !EQUAL(poBinFunc->GetType(),"Edsc_BinFunction") )
2754 2 : poBinFunc = new HFAEntry( hHFA, "#Bin_Function#", "Edsc_BinFunction", poEntry );
2755 :
2756 2 : poBinFunc->MakeData( 30 );
2757 2 : poBinFunc->SetIntField( "numBins", nNumBins );
2758 2 : poBinFunc->SetDoubleField( "minLimit", dMinLimit );
2759 2 : poBinFunc->SetDoubleField( "maxLimit", dMaxLimit );
2760 2 : poBinFunc->SetStringField( "binFunctionType", "linear" ); // we use always a linear
2761 :
2762 : // we need a child named histogram
2763 2 : HFAEntry * poHisto = poEntry->GetNamedChild( "Histogram" );
2764 2 : if( poHisto == NULL || !EQUAL(poHisto->GetType(),"Edsc_Column") )
2765 2 : poHisto = new HFAEntry( hHFA, "Histogram", "Edsc_Column", poEntry );
2766 :
2767 2 : poHisto->SetIntField( "numRows", nNumBins );
2768 : // allocate space for the bin values
2769 2 : GUInt32 nOffset = HFAAllocateSpace( hHFA, nNumBins*4 );
2770 2 : poHisto->SetIntField( "columnDataPtr", nOffset );
2771 2 : poHisto->SetStringField( "dataType", "integer" );
2772 2 : poHisto->SetIntField( "maxNumChars", 0 );
2773 : // write out histogram data
2774 2 : char * pszWork = pszBinValues;
2775 503 : for ( int nBin = 0; nBin < nNumBins; ++nBin )
2776 : {
2777 501 : char * pszEnd = strchr( pszWork, '|' );
2778 501 : if ( pszEnd != NULL )
2779 : {
2780 501 : *pszEnd = 0;
2781 501 : VSIFSeekL( hHFA->fp, nOffset + 4*nBin, SEEK_SET );
2782 501 : int nValue = atoi( pszWork );
2783 : HFAStandard( 4, &nValue );
2784 :
2785 501 : VSIFWriteL( (void *)&nValue, 1, 4, hHFA->fp );
2786 501 : pszWork = pszEnd + 1;
2787 : }
2788 : }
2789 : }
2790 2 : free( pszBinValues );
2791 : }
2792 :
2793 : /* -------------------------------------------------------------------- */
2794 : /* If we created a statistics node then try to create a */
2795 : /* StatisticsParameters node too. */
2796 : /* -------------------------------------------------------------------- */
2797 46 : if( bCreatedStatistics )
2798 : {
2799 : HFAEntry *poEntry =
2800 : new HFAEntry( hHFA, "StatisticsParameters",
2801 2 : "Eimg_StatisticsParameters830", poNode );
2802 :
2803 2 : poEntry->MakeData( 70 );
2804 : //poEntry->SetStringField( "BinFunction.binFunctionType", "linear" );
2805 :
2806 2 : poEntry->SetIntField( "SkipFactorX", 1 );
2807 2 : poEntry->SetIntField( "SkipFactorY", 1 );
2808 : }
2809 :
2810 : /* -------------------------------------------------------------------- */
2811 : /* Write out metadata items without a special place. */
2812 : /* -------------------------------------------------------------------- */
2813 46 : if( CSLCount( papszGDALMD) != 0 )
2814 : {
2815 44 : CPLErr eErr = HFASetGDALMetadata( hHFA, nBand, papszGDALMD );
2816 :
2817 44 : CSLDestroy( papszGDALMD );
2818 44 : return eErr;
2819 : }
2820 : else
2821 2 : return CE_Failure;
2822 : }
2823 :
2824 : /************************************************************************/
2825 : /* HFACreateSpillStack() */
2826 : /* */
2827 : /* Create a new stack of raster layers in the spill (.ige) */
2828 : /* file. Create the spill file if it didn't exist before. */
2829 : /************************************************************************/
2830 :
2831 7 : int HFACreateSpillStack( HFAInfo_t *psInfo, int nXSize, int nYSize,
2832 : int nLayers, int nBlockSize, int nDataType,
2833 : GIntBig *pnValidFlagsOffset,
2834 : GIntBig *pnDataOffset )
2835 :
2836 : {
2837 : /* -------------------------------------------------------------------- */
2838 : /* Form .ige filename. */
2839 : /* -------------------------------------------------------------------- */
2840 : char *pszFullFilename;
2841 :
2842 7 : if (nBlockSize <= 0)
2843 : {
2844 0 : CPLError(CE_Failure, CPLE_IllegalArg, "HFACreateSpillStack : nBlockXSize < 0");
2845 0 : return FALSE;
2846 : }
2847 :
2848 7 : if( psInfo->pszIGEFilename == NULL )
2849 : psInfo->pszIGEFilename =
2850 7 : CPLStrdup( CPLResetExtension( psInfo->pszFilename, "ige" ) );
2851 :
2852 : pszFullFilename =
2853 7 : CPLStrdup( CPLFormFilename( psInfo->pszPath, psInfo->pszIGEFilename, NULL ) );
2854 :
2855 : /* -------------------------------------------------------------------- */
2856 : /* Try and open it. If we fail, create it and write the magic */
2857 : /* header. */
2858 : /* -------------------------------------------------------------------- */
2859 : static const char *pszMagick = "ERDAS_IMG_EXTERNAL_RASTER";
2860 : FILE *fpVSIL;
2861 :
2862 7 : fpVSIL = VSIFOpenL( pszFullFilename, "r+b" );
2863 7 : if( fpVSIL == NULL )
2864 : {
2865 7 : fpVSIL = VSIFOpenL( pszFullFilename, "w+" );
2866 7 : if( fpVSIL == NULL )
2867 : {
2868 : CPLError( CE_Failure, CPLE_OpenFailed,
2869 : "Failed to create spill file %s.\n%s",
2870 0 : psInfo->pszIGEFilename, VSIStrerror( errno ) );
2871 0 : return FALSE;
2872 : }
2873 :
2874 7 : VSIFWriteL( (void *) pszMagick, 1, strlen(pszMagick)+1, fpVSIL );
2875 : }
2876 :
2877 7 : CPLFree( pszFullFilename );
2878 :
2879 : /* -------------------------------------------------------------------- */
2880 : /* Work out some details about the tiling scheme. */
2881 : /* -------------------------------------------------------------------- */
2882 : int nBlocksPerRow, nBlocksPerColumn, nBlocks, nBytesPerBlock;
2883 : int nBytesPerRow, nBlockMapSize, iFlagsSize;
2884 :
2885 7 : nBlocksPerRow = (nXSize + nBlockSize - 1) / nBlockSize;
2886 7 : nBlocksPerColumn = (nYSize + nBlockSize - 1) / nBlockSize;
2887 7 : nBlocks = nBlocksPerRow * nBlocksPerColumn;
2888 : nBytesPerBlock = (nBlockSize * nBlockSize
2889 7 : * HFAGetDataTypeBits(nDataType) + 7) / 8;
2890 :
2891 7 : nBytesPerRow = ( nBlocksPerRow + 7 ) / 8;
2892 7 : nBlockMapSize = nBytesPerRow * nBlocksPerColumn;
2893 7 : iFlagsSize = nBlockMapSize + 20;
2894 :
2895 : /* -------------------------------------------------------------------- */
2896 : /* Write stack prefix information. */
2897 : /* -------------------------------------------------------------------- */
2898 : GByte bUnknown;
2899 : GInt32 nValue32;
2900 :
2901 7 : VSIFSeekL( fpVSIL, 0, SEEK_END );
2902 :
2903 7 : bUnknown = 1;
2904 7 : VSIFWriteL( &bUnknown, 1, 1, fpVSIL );
2905 7 : nValue32 = nLayers;
2906 : HFAStandard( 4, &nValue32 );
2907 7 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2908 7 : nValue32 = nXSize;
2909 : HFAStandard( 4, &nValue32 );
2910 7 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2911 7 : nValue32 = nYSize;
2912 : HFAStandard( 4, &nValue32 );
2913 7 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2914 7 : nValue32 = nBlockSize;
2915 : HFAStandard( 4, &nValue32 );
2916 7 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2917 7 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2918 7 : bUnknown = 3;
2919 7 : VSIFWriteL( &bUnknown, 1, 1, fpVSIL );
2920 7 : bUnknown = 0;
2921 7 : VSIFWriteL( &bUnknown, 1, 1, fpVSIL );
2922 :
2923 : /* -------------------------------------------------------------------- */
2924 : /* Write out ValidFlags section(s). */
2925 : /* -------------------------------------------------------------------- */
2926 : unsigned char *pabyBlockMap;
2927 : int iBand;
2928 :
2929 7 : *pnValidFlagsOffset = VSIFTellL( fpVSIL );
2930 :
2931 7 : pabyBlockMap = (unsigned char *) VSIMalloc( nBlockMapSize );
2932 7 : if (pabyBlockMap == NULL)
2933 : {
2934 0 : CPLError(CE_Failure, CPLE_OutOfMemory, "HFACreateSpillStack : Out of memory");
2935 0 : VSIFCloseL( fpVSIL );
2936 0 : return FALSE;
2937 : }
2938 :
2939 7 : memset( pabyBlockMap, 0xff, nBlockMapSize );
2940 20 : for ( iBand = 0; iBand < nLayers; iBand++ )
2941 : {
2942 : int i, iRemainder;
2943 :
2944 13 : nValue32 = 1; // Unknown
2945 : HFAStandard( 4, &nValue32 );
2946 13 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2947 13 : nValue32 = 0; // Unknown
2948 13 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2949 13 : nValue32 = nBlocksPerColumn;
2950 : HFAStandard( 4, &nValue32 );
2951 13 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2952 13 : nValue32 = nBlocksPerRow;
2953 : HFAStandard( 4, &nValue32 );
2954 13 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2955 13 : nValue32 = 0x30000; // Unknown
2956 : HFAStandard( 4, &nValue32 );
2957 13 : VSIFWriteL( &nValue32, 4, 1, fpVSIL );
2958 :
2959 13 : iRemainder = nBlocksPerRow % 8;
2960 : CPLDebug( "HFACreate",
2961 : "Block map size %d, bytes per row %d, remainder %d.",
2962 13 : nBlockMapSize, nBytesPerRow, iRemainder );
2963 13 : if ( iRemainder )
2964 : {
2965 26 : for ( i = nBytesPerRow - 1; i < nBlockMapSize; i+=nBytesPerRow )
2966 13 : pabyBlockMap[i] = (GByte) ((1<<iRemainder) - 1);
2967 : }
2968 :
2969 13 : VSIFWriteL( pabyBlockMap, 1, nBlockMapSize, fpVSIL );
2970 : }
2971 7 : CPLFree(pabyBlockMap);
2972 7 : pabyBlockMap = NULL;
2973 :
2974 : /* -------------------------------------------------------------------- */
2975 : /* Extend the file to account for all the imagery space. */
2976 : /* -------------------------------------------------------------------- */
2977 : GIntBig nTileDataSize = ((GIntBig) nBytesPerBlock)
2978 7 : * nBlocksPerRow * nBlocksPerColumn * nLayers;
2979 :
2980 7 : *pnDataOffset = VSIFTellL( fpVSIL );
2981 :
2982 7 : if( VSIFSeekL( fpVSIL, nTileDataSize - 1 + *pnDataOffset, SEEK_SET ) != 0
2983 : || VSIFWriteL( (void *) "", 1, 1, fpVSIL ) != 1 )
2984 : {
2985 : CPLError( CE_Failure, CPLE_FileIO,
2986 : "Failed to extend %s to full size (%g bytes),\n"
2987 : "likely out of disk space.\n%s",
2988 : psInfo->pszIGEFilename,
2989 : (double) nTileDataSize - 1 + *pnDataOffset,
2990 0 : VSIStrerror( errno ) );
2991 :
2992 0 : VSIFCloseL( fpVSIL );
2993 0 : return FALSE;
2994 : }
2995 :
2996 7 : VSIFCloseL( fpVSIL );
2997 :
2998 7 : return TRUE;
2999 : }
3000 :
3001 : /************************************************************************/
3002 : /* HFAReadAndValidatePoly() */
3003 : /************************************************************************/
3004 :
3005 3 : static int HFAReadAndValidatePoly( HFAEntry *poTarget,
3006 : const char *pszName,
3007 : Efga_Polynomial *psRetPoly )
3008 :
3009 : {
3010 3 : CPLString osFldName;
3011 :
3012 3 : memset( psRetPoly, 0, sizeof(Efga_Polynomial) );
3013 :
3014 3 : osFldName.Printf( "%sorder", pszName );
3015 3 : psRetPoly->order = poTarget->GetIntField(osFldName);
3016 :
3017 3 : if( psRetPoly->order < 1 || psRetPoly->order > 3 )
3018 0 : return FALSE;
3019 :
3020 : /* -------------------------------------------------------------------- */
3021 : /* Validate that things are in a "well known" form. */
3022 : /* -------------------------------------------------------------------- */
3023 : int numdimtransform, numdimpolynomial, termcount;
3024 :
3025 3 : osFldName.Printf( "%snumdimtransform", pszName );
3026 3 : numdimtransform = poTarget->GetIntField(osFldName);
3027 :
3028 3 : osFldName.Printf( "%snumdimpolynomial", pszName );
3029 3 : numdimpolynomial = poTarget->GetIntField(osFldName);
3030 :
3031 3 : osFldName.Printf( "%stermcount", pszName );
3032 3 : termcount = poTarget->GetIntField(osFldName);
3033 :
3034 3 : if( numdimtransform != 2 || numdimpolynomial != 2 )
3035 0 : return FALSE;
3036 :
3037 3 : if( (psRetPoly->order == 1 && termcount != 3)
3038 : || (psRetPoly->order == 2 && termcount != 6)
3039 : || (psRetPoly->order == 3 && termcount != 10) )
3040 0 : return FALSE;
3041 :
3042 : // we don't check the exponent organization for now. Hopefully
3043 : // it is always standard.
3044 :
3045 : /* -------------------------------------------------------------------- */
3046 : /* Get coefficients. */
3047 : /* -------------------------------------------------------------------- */
3048 : int i;
3049 :
3050 43 : for( i = 0; i < termcount*2 - 2; i++ )
3051 : {
3052 40 : osFldName.Printf( "%spolycoefmtx[%d]", pszName, i );
3053 40 : psRetPoly->polycoefmtx[i] = poTarget->GetDoubleField(osFldName);
3054 : }
3055 :
3056 9 : for( i = 0; i < 2; i++ )
3057 : {
3058 6 : osFldName.Printf( "%spolycoefvector[%d]", pszName, i );
3059 6 : psRetPoly->polycoefvector[i] = poTarget->GetDoubleField(osFldName);
3060 : }
3061 :
3062 3 : return TRUE;
3063 : }
3064 :
3065 : /************************************************************************/
3066 : /* HFAReadXFormStack() */
3067 : /************************************************************************/
3068 :
3069 :
3070 175 : int HFAReadXFormStack( HFAHandle hHFA,
3071 : Efga_Polynomial **ppasPolyListForward,
3072 : Efga_Polynomial **ppasPolyListReverse )
3073 :
3074 : {
3075 175 : if( hHFA->nBands == 0 )
3076 0 : return 0;
3077 :
3078 : /* -------------------------------------------------------------------- */
3079 : /* Get the HFA node. */
3080 : /* -------------------------------------------------------------------- */
3081 : HFAEntry *poXFormHeader;
3082 :
3083 175 : poXFormHeader = hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm" );
3084 175 : if( poXFormHeader == NULL )
3085 174 : return 0;
3086 :
3087 : /* -------------------------------------------------------------------- */
3088 : /* Loop over children, collecting XForms. */
3089 : /* -------------------------------------------------------------------- */
3090 : HFAEntry *poXForm;
3091 1 : int nStepCount = 0;
3092 1 : *ppasPolyListForward = NULL;
3093 1 : *ppasPolyListReverse = NULL;
3094 :
3095 3 : for( poXForm = poXFormHeader->GetChild();
3096 : poXForm != NULL;
3097 : poXForm = poXForm->GetNext() )
3098 : {
3099 2 : int bSuccess = FALSE;
3100 : Efga_Polynomial sForward, sReverse;
3101 :
3102 2 : if( EQUAL(poXForm->GetType(),"Efga_Polynomial") )
3103 : {
3104 : bSuccess =
3105 1 : HFAReadAndValidatePoly( poXForm, "", &sForward );
3106 :
3107 1 : if( bSuccess )
3108 : {
3109 : double adfGT[6], adfInvGT[6];
3110 :
3111 1 : adfGT[0] = sForward.polycoefvector[0];
3112 1 : adfGT[1] = sForward.polycoefmtx[0];
3113 1 : adfGT[2] = sForward.polycoefmtx[2];
3114 1 : adfGT[3] = sForward.polycoefvector[1];
3115 1 : adfGT[4] = sForward.polycoefmtx[1];
3116 1 : adfGT[5] = sForward.polycoefmtx[3];
3117 :
3118 1 : bSuccess = HFAInvGeoTransform( adfGT, adfInvGT );
3119 :
3120 1 : memset( &sReverse, 0, sizeof(sReverse) );
3121 :
3122 1 : sReverse.order = sForward.order;
3123 1 : sReverse.polycoefvector[0] = adfInvGT[0];
3124 1 : sReverse.polycoefmtx[0] = adfInvGT[1];
3125 1 : sReverse.polycoefmtx[2] = adfInvGT[2];
3126 1 : sReverse.polycoefvector[1] = adfInvGT[3];
3127 1 : sReverse.polycoefmtx[1] = adfInvGT[4];
3128 1 : sReverse.polycoefmtx[3] = adfInvGT[5];
3129 : }
3130 : }
3131 1 : else if( EQUAL(poXForm->GetType(),"GM_PolyPair") )
3132 : {
3133 : bSuccess =
3134 1 : HFAReadAndValidatePoly( poXForm, "forward.", &sForward );
3135 : bSuccess = bSuccess &&
3136 1 : HFAReadAndValidatePoly( poXForm, "reverse.", &sReverse );
3137 : }
3138 :
3139 2 : if( bSuccess )
3140 : {
3141 2 : nStepCount++;
3142 : *ppasPolyListForward = (Efga_Polynomial *)
3143 : CPLRealloc( *ppasPolyListForward,
3144 2 : sizeof(Efga_Polynomial) * nStepCount);
3145 : memcpy( *ppasPolyListForward + nStepCount - 1,
3146 2 : &sForward, sizeof(sForward) );
3147 :
3148 : *ppasPolyListReverse = (Efga_Polynomial *)
3149 : CPLRealloc( *ppasPolyListReverse,
3150 2 : sizeof(Efga_Polynomial) * nStepCount);
3151 : memcpy( *ppasPolyListReverse + nStepCount - 1,
3152 2 : &sReverse, sizeof(sReverse) );
3153 : }
3154 : }
3155 :
3156 1 : return nStepCount;
3157 : }
3158 :
3159 : /************************************************************************/
3160 : /* HFAEvaluateXFormStack() */
3161 : /************************************************************************/
3162 :
3163 36 : int HFAEvaluateXFormStack( int nStepCount, int bForward,
3164 : Efga_Polynomial *pasPolyList,
3165 : double *pdfX, double *pdfY )
3166 :
3167 : {
3168 : int iStep;
3169 :
3170 108 : for( iStep = 0; iStep < nStepCount; iStep++ )
3171 : {
3172 : double dfXOut, dfYOut;
3173 : Efga_Polynomial *psStep;
3174 :
3175 72 : if( bForward )
3176 0 : psStep = pasPolyList + iStep;
3177 : else
3178 72 : psStep = pasPolyList + nStepCount - iStep - 1;
3179 :
3180 72 : if( psStep->order == 1 )
3181 : {
3182 36 : dfXOut = psStep->polycoefvector[0]
3183 36 : + psStep->polycoefmtx[0] * *pdfX
3184 36 : + psStep->polycoefmtx[2] * *pdfY;
3185 :
3186 36 : dfYOut = psStep->polycoefvector[1]
3187 36 : + psStep->polycoefmtx[1] * *pdfX
3188 36 : + psStep->polycoefmtx[3] * *pdfY;
3189 :
3190 36 : *pdfX = dfXOut;
3191 36 : *pdfY = dfYOut;
3192 : }
3193 36 : else if( psStep->order == 2 )
3194 : {
3195 0 : dfXOut = psStep->polycoefvector[0]
3196 0 : + psStep->polycoefmtx[0] * *pdfX
3197 0 : + psStep->polycoefmtx[2] * *pdfY
3198 0 : + psStep->polycoefmtx[4] * *pdfX * *pdfX
3199 0 : + psStep->polycoefmtx[6] * *pdfX * *pdfY
3200 0 : + psStep->polycoefmtx[8] * *pdfY * *pdfY;
3201 0 : dfYOut = psStep->polycoefvector[1]
3202 0 : + psStep->polycoefmtx[1] * *pdfX
3203 0 : + psStep->polycoefmtx[3] * *pdfY
3204 0 : + psStep->polycoefmtx[5] * *pdfX * *pdfX
3205 0 : + psStep->polycoefmtx[7] * *pdfX * *pdfY
3206 0 : + psStep->polycoefmtx[9] * *pdfY * *pdfY;
3207 :
3208 0 : *pdfX = dfXOut;
3209 0 : *pdfY = dfYOut;
3210 : }
3211 36 : else if( psStep->order == 3 )
3212 : {
3213 36 : dfXOut = psStep->polycoefvector[0]
3214 36 : + psStep->polycoefmtx[ 0] * *pdfX
3215 36 : + psStep->polycoefmtx[ 2] * *pdfY
3216 36 : + psStep->polycoefmtx[ 4] * *pdfX * *pdfX
3217 36 : + psStep->polycoefmtx[ 6] * *pdfX * *pdfY
3218 36 : + psStep->polycoefmtx[ 8] * *pdfY * *pdfY
3219 36 : + psStep->polycoefmtx[10] * *pdfX * *pdfX * *pdfX
3220 36 : + psStep->polycoefmtx[12] * *pdfX * *pdfX * *pdfY
3221 36 : + psStep->polycoefmtx[14] * *pdfX * *pdfY * *pdfY
3222 288 : + psStep->polycoefmtx[16] * *pdfY * *pdfY * *pdfY;
3223 36 : dfYOut = psStep->polycoefvector[1]
3224 36 : + psStep->polycoefmtx[ 1] * *pdfX
3225 36 : + psStep->polycoefmtx[ 3] * *pdfY
3226 36 : + psStep->polycoefmtx[ 5] * *pdfX * *pdfX
3227 36 : + psStep->polycoefmtx[ 7] * *pdfX * *pdfY
3228 36 : + psStep->polycoefmtx[ 9] * *pdfY * *pdfY
3229 36 : + psStep->polycoefmtx[11] * *pdfX * *pdfX * *pdfX
3230 36 : + psStep->polycoefmtx[13] * *pdfX * *pdfX * *pdfY
3231 36 : + psStep->polycoefmtx[15] * *pdfX * *pdfY * *pdfY
3232 288 : + psStep->polycoefmtx[17] * *pdfY * *pdfY * *pdfY;
3233 :
3234 36 : *pdfX = dfXOut;
3235 36 : *pdfY = dfYOut;
3236 : }
3237 : else
3238 0 : return FALSE;
3239 : }
3240 :
3241 36 : return TRUE;
3242 : }
3243 :
3244 : /************************************************************************/
3245 : /* HFAWriteXFormStack() */
3246 : /************************************************************************/
3247 :
3248 2 : CPLErr HFAWriteXFormStack( HFAHandle hHFA, int nBand, int nXFormCount,
3249 : Efga_Polynomial **ppasPolyListForward,
3250 : Efga_Polynomial **ppasPolyListReverse )
3251 :
3252 : {
3253 2 : if( nXFormCount == 0 )
3254 0 : return CE_None;
3255 :
3256 2 : if( ppasPolyListForward[0]->order != 1 )
3257 : {
3258 : CPLError( CE_Failure, CPLE_AppDefined,
3259 0 : "For now HFAWriteXFormStack() only supports order 1 polynomials" );
3260 0 : return CE_Failure;
3261 : }
3262 :
3263 2 : if( nBand < 0 || nBand > hHFA->nBands )
3264 0 : return CE_Failure;
3265 :
3266 : /* -------------------------------------------------------------------- */
3267 : /* If no band number is provided, operate on all bands. */
3268 : /* -------------------------------------------------------------------- */
3269 2 : if( nBand == 0 )
3270 : {
3271 1 : CPLErr eErr = CE_None;
3272 :
3273 2 : for( nBand = 1; nBand <= hHFA->nBands; nBand++ )
3274 : {
3275 : eErr = HFAWriteXFormStack( hHFA, nBand, nXFormCount,
3276 : ppasPolyListForward,
3277 1 : ppasPolyListReverse );
3278 1 : if( eErr != CE_None )
3279 0 : return eErr;
3280 : }
3281 :
3282 1 : return eErr;
3283 : }
3284 :
3285 : /* -------------------------------------------------------------------- */
3286 : /* Fetch our band node. */
3287 : /* -------------------------------------------------------------------- */
3288 1 : HFAEntry *poBandNode = hHFA->papoBand[nBand-1]->poNode;
3289 : HFAEntry *poXFormHeader;
3290 :
3291 1 : poXFormHeader = poBandNode->GetNamedChild( "MapToPixelXForm" );
3292 1 : if( poXFormHeader == NULL )
3293 : {
3294 : poXFormHeader = new HFAEntry( hHFA, "MapToPixelXForm",
3295 1 : "Exfr_GenericXFormHeader", poBandNode );
3296 1 : poXFormHeader->MakeData( 23 );
3297 1 : poXFormHeader->SetPosition();
3298 1 : poXFormHeader->SetStringField( "titleList.string", "Affine" );
3299 : }
3300 :
3301 : /* -------------------------------------------------------------------- */
3302 : /* Loop over XForms. */
3303 : /* -------------------------------------------------------------------- */
3304 2 : for( int iXForm = 0; iXForm < nXFormCount; iXForm++ )
3305 : {
3306 1 : Efga_Polynomial *psForward = *ppasPolyListForward + iXForm;
3307 1 : CPLString osXFormName;
3308 1 : osXFormName.Printf( "XForm%d", iXForm );
3309 :
3310 1 : HFAEntry *poXForm = poXFormHeader->GetNamedChild( osXFormName );
3311 :
3312 1 : if( poXForm == NULL )
3313 : {
3314 : poXForm = new HFAEntry( hHFA, osXFormName, "Efga_Polynomial",
3315 1 : poXFormHeader );
3316 1 : poXForm->MakeData( 136 );
3317 1 : poXForm->SetPosition();
3318 : }
3319 :
3320 1 : poXForm->SetIntField( "order", 1 );
3321 1 : poXForm->SetIntField( "numdimtransform", 2 );
3322 1 : poXForm->SetIntField( "numdimpolynomial", 2 );
3323 1 : poXForm->SetIntField( "termcount", 3 );
3324 1 : poXForm->SetIntField( "exponentlist[0]", 0 );
3325 1 : poXForm->SetIntField( "exponentlist[1]", 0 );
3326 1 : poXForm->SetIntField( "exponentlist[2]", 1 );
3327 1 : poXForm->SetIntField( "exponentlist[3]", 0 );
3328 1 : poXForm->SetIntField( "exponentlist[4]", 0 );
3329 1 : poXForm->SetIntField( "exponentlist[5]", 1 );
3330 :
3331 1 : poXForm->SetIntField( "polycoefmtx[-3]", EPT_f64 );
3332 1 : poXForm->SetIntField( "polycoefmtx[-2]", 2 );
3333 1 : poXForm->SetIntField( "polycoefmtx[-1]", 2 );
3334 : poXForm->SetDoubleField( "polycoefmtx[0]",
3335 1 : psForward->polycoefmtx[0] );
3336 : poXForm->SetDoubleField( "polycoefmtx[1]",
3337 1 : psForward->polycoefmtx[1] );
3338 : poXForm->SetDoubleField( "polycoefmtx[2]",
3339 1 : psForward->polycoefmtx[2] );
3340 : poXForm->SetDoubleField( "polycoefmtx[3]",
3341 1 : psForward->polycoefmtx[3] );
3342 :
3343 1 : poXForm->SetIntField( "polycoefvector[-3]", EPT_f64 );
3344 1 : poXForm->SetIntField( "polycoefvector[-2]", 1 );
3345 1 : poXForm->SetIntField( "polycoefvector[-1]", 2 );
3346 : poXForm->SetDoubleField( "polycoefvector[0]",
3347 1 : psForward->polycoefvector[0] );
3348 : poXForm->SetDoubleField( "polycoefvector[1]",
3349 1 : psForward->polycoefvector[1] );
3350 : }
3351 :
3352 1 : return CE_None;
3353 : }
3354 :
3355 : /************************************************************************/
3356 : /* HFAReadCameraModel() */
3357 : /************************************************************************/
3358 :
3359 275 : char **HFAReadCameraModel( HFAHandle hHFA )
3360 :
3361 : {
3362 275 : if( hHFA->nBands == 0 )
3363 0 : return NULL;
3364 :
3365 : /* -------------------------------------------------------------------- */
3366 : /* Get the camera model node, and confirm it's type. */
3367 : /* -------------------------------------------------------------------- */
3368 : HFAEntry *poXForm;
3369 :
3370 : poXForm =
3371 275 : hHFA->papoBand[0]->poNode->GetNamedChild( "MapToPixelXForm.XForm0" );
3372 275 : if( poXForm == NULL )
3373 270 : return NULL;
3374 :
3375 5 : if( !EQUAL(poXForm->GetType(),"Camera_ModelX") )
3376 5 : return NULL;
3377 :
3378 : /* -------------------------------------------------------------------- */
3379 : /* Convert the values to metadata. */
3380 : /* -------------------------------------------------------------------- */
3381 : const char *pszValue;
3382 : int i;
3383 0 : char **papszMD = NULL;
3384 : static const char *apszFields[] = {
3385 : "direction", "refType", "demsource", "PhotoDirection", "RotationSystem",
3386 : "demfilename", "demzunits",
3387 : "forSrcAffine[0]", "forSrcAffine[1]", "forSrcAffine[2]",
3388 : "forSrcAffine[3]", "forSrcAffine[4]", "forSrcAffine[5]",
3389 : "forDstAffine[0]", "forDstAffine[1]", "forDstAffine[2]",
3390 : "forDstAffine[3]", "forDstAffine[4]", "forDstAffine[5]",
3391 : "invSrcAffine[0]", "invSrcAffine[1]", "invSrcAffine[2]",
3392 : "invSrcAffine[3]", "invSrcAffine[4]", "invSrcAffine[5]",
3393 : "invDstAffine[0]", "invDstAffine[1]", "invDstAffine[2]",
3394 : "invDstAffine[3]", "invDstAffine[4]", "invDstAffine[5]",
3395 : "z_mean", "lat0", "lon0",
3396 : "coeffs[0]", "coeffs[1]", "coeffs[2]",
3397 : "coeffs[3]", "coeffs[4]", "coeffs[5]",
3398 : "coeffs[6]", "coeffs[7]", "coeffs[8]",
3399 : "LensDistortion[0]", "LensDistortion[1]", "LensDistortion[2]",
3400 : NULL };
3401 :
3402 0 : for( i = 0; apszFields[i] != NULL; i++ )
3403 : {
3404 0 : pszValue = poXForm->GetStringField( apszFields[i] );
3405 0 : if( pszValue == NULL )
3406 0 : pszValue = "";
3407 :
3408 0 : papszMD = CSLSetNameValue( papszMD, apszFields[i], pszValue );
3409 : }
3410 :
3411 : /* -------------------------------------------------------------------- */
3412 : /* Create a pseudo-entry for the MIFObject with the */
3413 : /* outputProjection. */
3414 : /* -------------------------------------------------------------------- */
3415 0 : HFAEntry *poProjInfo = new HFAEntry( poXForm, "outputProjection" );
3416 :
3417 : /* -------------------------------------------------------------------- */
3418 : /* Fetch the datum. */
3419 : /* -------------------------------------------------------------------- */
3420 : Eprj_Datum sDatum;
3421 :
3422 0 : memset( &sDatum, 0, sizeof(sDatum));
3423 :
3424 : sDatum.datumname =
3425 0 : (char *) poProjInfo->GetStringField("earthModel.datum.datumname");
3426 : sDatum.type = (Eprj_DatumType) poProjInfo->GetIntField(
3427 0 : "earthModel.datum.type");
3428 :
3429 0 : for( i = 0; i < 7; i++ )
3430 : {
3431 : char szFieldName[60];
3432 :
3433 0 : sprintf( szFieldName, "earthModel.datum.params[%d]", i );
3434 0 : sDatum.params[i] = poProjInfo->GetDoubleField(szFieldName);
3435 : }
3436 :
3437 : sDatum.gridname = (char *)
3438 0 : poProjInfo->GetStringField("earthModel.datum.gridname");
3439 :
3440 : /* -------------------------------------------------------------------- */
3441 : /* Fetch the projection parameters. */
3442 : /* -------------------------------------------------------------------- */
3443 : Eprj_ProParameters sPro;
3444 :
3445 0 : memset( &sPro, 0, sizeof(sPro) );
3446 :
3447 0 : sPro.proType = (Eprj_ProType) poProjInfo->GetIntField("projectionObject.proType");
3448 0 : sPro.proNumber = poProjInfo->GetIntField("projectionObject.proNumber");
3449 0 : sPro.proExeName = (char *) poProjInfo->GetStringField("projectionObject.proExeName");
3450 0 : sPro.proName = (char *) poProjInfo->GetStringField("projectionObject.proName");
3451 0 : sPro.proZone = poProjInfo->GetIntField("projectionObject.proZone");
3452 :
3453 0 : for( i = 0; i < 15; i++ )
3454 : {
3455 : char szFieldName[30];
3456 :
3457 0 : sprintf( szFieldName, "projectionObject.proParams[%d]", i );
3458 0 : sPro.proParams[i] = poProjInfo->GetDoubleField(szFieldName);
3459 : }
3460 :
3461 : /* -------------------------------------------------------------------- */
3462 : /* Fetch the spheroid. */
3463 : /* -------------------------------------------------------------------- */
3464 : sPro.proSpheroid.sphereName = (char *)
3465 0 : poProjInfo->GetStringField("earthModel.proSpheroid.sphereName");
3466 0 : sPro.proSpheroid.a = poProjInfo->GetDoubleField("earthModel.proSpheroid.a");
3467 0 : sPro.proSpheroid.b = poProjInfo->GetDoubleField("earthModel.proSpheroid.b");
3468 : sPro.proSpheroid.eSquared =
3469 0 : poProjInfo->GetDoubleField("earthModel.proSpheroid.eSquared");
3470 : sPro.proSpheroid.radius =
3471 0 : poProjInfo->GetDoubleField("earthModel.proSpheroid.radius");
3472 :
3473 : /* -------------------------------------------------------------------- */
3474 : /* Fetch the projection info. */
3475 : /* -------------------------------------------------------------------- */
3476 : char *pszProjection;
3477 :
3478 : // poProjInfo->DumpFieldValues( stdout, "" );
3479 :
3480 0 : pszProjection = HFAPCSStructToWKT( &sDatum, &sPro, NULL, NULL );
3481 :
3482 0 : if( pszProjection )
3483 : {
3484 : papszMD =
3485 0 : CSLSetNameValue( papszMD, "outputProjection", pszProjection );
3486 0 : CPLFree( pszProjection );
3487 : }
3488 :
3489 0 : delete poProjInfo;
3490 :
3491 : /* -------------------------------------------------------------------- */
3492 : /* Fetch the horizontal units. */
3493 : /* -------------------------------------------------------------------- */
3494 0 : pszValue = poXForm->GetStringField( "outputHorizontalUnits.string" );
3495 0 : if( pszValue == NULL )
3496 0 : pszValue = "";
3497 :
3498 0 : papszMD = CSLSetNameValue( papszMD, "outputHorizontalUnits", pszValue );
3499 :
3500 : /* -------------------------------------------------------------------- */
3501 : /* Fetch the elevationinfo. */
3502 : /* -------------------------------------------------------------------- */
3503 0 : HFAEntry *poElevInfo = new HFAEntry( poXForm, "outputElevationInfo" );
3504 : //poElevInfo->DumpFieldValues( stdout, "" );
3505 :
3506 0 : if( poElevInfo->GetDataSize() != 0 )
3507 : {
3508 : static const char *apszEFields[] = {
3509 : "verticalDatum.datumname",
3510 : "verticalDatum.type",
3511 : "elevationUnit",
3512 : "elevationType",
3513 : NULL };
3514 :
3515 0 : for( i = 0; apszEFields[i] != NULL; i++ )
3516 : {
3517 0 : pszValue = poElevInfo->GetStringField( apszEFields[i] );
3518 0 : if( pszValue == NULL )
3519 0 : pszValue = "";
3520 :
3521 0 : papszMD = CSLSetNameValue( papszMD, apszEFields[i], pszValue );
3522 : }
3523 : }
3524 :
3525 0 : delete poElevInfo;
3526 :
3527 0 : return papszMD;
3528 : }
3529 :
3530 : /************************************************************************/
3531 : /* HFASetGeoTransform() */
3532 : /* */
3533 : /* Set a MapInformation and XForm block. Allows for rotated */
3534 : /* and shared geotransforms. */
3535 : /************************************************************************/
3536 :
3537 1 : CPLErr HFASetGeoTransform( HFAHandle hHFA,
3538 : const char *pszProName,
3539 : const char *pszUnits,
3540 : double *padfGeoTransform )
3541 :
3542 : {
3543 : /* -------------------------------------------------------------------- */
3544 : /* Write MapInformation. */
3545 : /* -------------------------------------------------------------------- */
3546 : int nBand;
3547 :
3548 2 : for( nBand = 1; nBand <= hHFA->nBands; nBand++ )
3549 : {
3550 1 : HFAEntry *poBandNode = hHFA->papoBand[nBand-1]->poNode;
3551 :
3552 1 : HFAEntry *poMI = poBandNode->GetNamedChild( "MapInformation" );
3553 1 : if( poMI == NULL )
3554 : {
3555 : poMI = new HFAEntry( hHFA, "MapInformation",
3556 1 : "Eimg_MapInformation", poBandNode );
3557 1 : poMI->MakeData( 18 + strlen(pszProName) + strlen(pszUnits) );
3558 1 : poMI->SetPosition();
3559 : }
3560 :
3561 1 : poMI->SetStringField( "projection.string", pszProName );
3562 1 : poMI->SetStringField( "units.string", pszUnits );
3563 : }
3564 :
3565 : /* -------------------------------------------------------------------- */
3566 : /* Write XForm. */
3567 : /* -------------------------------------------------------------------- */
3568 : Efga_Polynomial sForward, sReverse;
3569 : double adfAdjTransform[6], adfRevTransform[6];
3570 :
3571 : // Offset by half pixel.
3572 :
3573 1 : memcpy( adfAdjTransform, padfGeoTransform, sizeof(double) * 6 );
3574 1 : adfAdjTransform[0] += adfAdjTransform[1] * 0.5;
3575 1 : adfAdjTransform[0] += adfAdjTransform[2] * 0.5;
3576 1 : adfAdjTransform[3] += adfAdjTransform[4] * 0.5;
3577 1 : adfAdjTransform[3] += adfAdjTransform[5] * 0.5;
3578 :
3579 : // Invert
3580 1 : HFAInvGeoTransform( adfAdjTransform, adfRevTransform );
3581 :
3582 : // Assign to polynomial object.
3583 :
3584 1 : sForward.order = 1;
3585 1 : sForward.polycoefvector[0] = adfRevTransform[0];
3586 1 : sForward.polycoefmtx[0] = adfRevTransform[1];
3587 1 : sForward.polycoefmtx[1] = adfRevTransform[4];
3588 1 : sForward.polycoefvector[1] = adfRevTransform[3];
3589 1 : sForward.polycoefmtx[2] = adfRevTransform[2];
3590 1 : sForward.polycoefmtx[3] = adfRevTransform[5];
3591 :
3592 1 : sReverse = sForward;
3593 1 : Efga_Polynomial *psForward=&sForward, *psReverse=&sReverse;
3594 :
3595 1 : return HFAWriteXFormStack( hHFA, 0, 1, &psForward, &psReverse );
3596 : }
|
