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