1 /******************************************************************************
2 *
3 * Project: Polarimetric Workstation
4 * Purpose: Convair PolGASP data (.img/.hdr format).
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
9 * Copyright (c) 2009, Even Rouault <even dot rouault at spatialys.com>
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included
19 * in all copies or substantial portions of the Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include "cpl_string.h"
31 #include "gdal_frmts.h"
32 #include "ogr_spatialref.h"
33 #include "rawdataset.h"
34
35 #include <vector>
36
37 CPL_CVSID("$Id: cpgdataset.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
38
39 enum Interleave { BSQ, BIL, BIP };
40
41 /************************************************************************/
42 /* ==================================================================== */
43 /* CPGDataset */
44 /* ==================================================================== */
45 /************************************************************************/
46
47 class SIRC_QSLCRasterBand;
48 class CPG_STOKESRasterBand;
49
50 class CPGDataset final: public RawDataset
51 {
52 friend class SIRC_QSLCRasterBand;
53 friend class CPG_STOKESRasterBand;
54
55 VSILFILE *afpImage[4];
56 std::vector<CPLString> aosImageFilenames{};
57
58 int nGCPCount;
59 GDAL_GCP *pasGCPList;
60 char *pszGCPProjection{};
61
62 double adfGeoTransform[6];
63 char *pszProjection{};
64
65 int nLoadedStokesLine;
66 float *padfStokesMatrix;
67
68 int nInterleave;
69 static int AdjustFilename( char **, const char *, const char * );
70 static int FindType1( const char *pszWorkname );
71 static int FindType2( const char *pszWorkname );
72 #ifdef notdef
73 static int FindType3( const char *pszWorkname );
74 #endif
75 static GDALDataset *InitializeType1Or2Dataset( const char *pszWorkname );
76 #ifdef notdef
77 static GDALDataset *InitializeType3Dataset( const char *pszWorkname );
78 #endif
79 CPLErr LoadStokesLine( int iLine, int bNativeOrder );
80
81 CPL_DISALLOW_COPY_ASSIGN(CPGDataset)
82
83 public:
84 CPGDataset();
85 ~CPGDataset() override;
86
87 int GetGCPCount() override;
88 const char *_GetGCPProjection() override;
GetGCPSpatialRef() const89 const OGRSpatialReference* GetGCPSpatialRef() const override {
90 return GetGCPSpatialRefFromOldGetGCPProjection();
91 }
92 const GDAL_GCP *GetGCPs() override;
93
94 const char *_GetProjectionRef() override;
GetSpatialRef() const95 const OGRSpatialReference* GetSpatialRef() const override {
96 return GetSpatialRefFromOldGetProjectionRef();
97 }
98 CPLErr GetGeoTransform( double * ) override;
99
100 char **GetFileList() override;
101
102 static GDALDataset *Open( GDALOpenInfo * );
103 };
104
105 /************************************************************************/
106 /* CPGDataset() */
107 /************************************************************************/
108
CPGDataset()109 CPGDataset::CPGDataset() :
110 nGCPCount(0),
111 pasGCPList(nullptr),
112 nLoadedStokesLine(-1),
113 padfStokesMatrix(nullptr),
114 nInterleave(0)
115 {
116 pszProjection = CPLStrdup("");
117 pszGCPProjection = CPLStrdup("");
118 adfGeoTransform[0] = 0.0;
119 adfGeoTransform[1] = 1.0;
120 adfGeoTransform[2] = 0.0;
121 adfGeoTransform[3] = 0.0;
122 adfGeoTransform[4] = 0.0;
123 adfGeoTransform[5] = 1.0;
124
125 for( int iBand = 0; iBand < 4; iBand++ )
126 afpImage[iBand] = nullptr;
127 }
128
129 /************************************************************************/
130 /* ~CPGDataset() */
131 /************************************************************************/
132
~CPGDataset()133 CPGDataset::~CPGDataset()
134
135 {
136 FlushCache();
137
138 for( int iBand = 0; iBand < 4; iBand++ )
139 {
140 if( afpImage[iBand] != nullptr )
141 VSIFCloseL( afpImage[iBand] );
142 }
143
144 if( nGCPCount > 0 )
145 {
146 GDALDeinitGCPs( nGCPCount, pasGCPList );
147 CPLFree( pasGCPList );
148 }
149
150 CPLFree( pszProjection );
151 CPLFree( pszGCPProjection );
152 CPLFree( padfStokesMatrix );
153 }
154
155 /************************************************************************/
156 /* GetFileList() */
157 /************************************************************************/
158
GetFileList()159 char **CPGDataset::GetFileList()
160
161 {
162 char **papszFileList = RawDataset::GetFileList();
163 for( size_t i = 0; i < aosImageFilenames.size(); ++i )
164 papszFileList = CSLAddString(papszFileList, aosImageFilenames[i]);
165 return papszFileList;
166 }
167
168 /************************************************************************/
169 /* ==================================================================== */
170 /* SIRC_QSLCPRasterBand */
171 /* ==================================================================== */
172 /************************************************************************/
173
174 class SIRC_QSLCRasterBand final: public GDALRasterBand
175 {
176 friend class CPGDataset;
177
178 public:
179 SIRC_QSLCRasterBand( CPGDataset *, int, GDALDataType );
~SIRC_QSLCRasterBand()180 ~SIRC_QSLCRasterBand() override {}
181
182 CPLErr IReadBlock( int, int, void * ) override;
183 };
184
185 constexpr int M11 = 0;
186 //constexpr int M12 = 1;
187 constexpr int M13 = 2;
188 constexpr int M14 = 3;
189 //constexpr int M21 = 4;
190 constexpr int M22 = 5;
191 constexpr int M23 = 6;
192 constexpr int M24 = 7;
193 constexpr int M31 = 8;
194 constexpr int M32 = 9;
195 constexpr int M33 = 10;
196 constexpr int M34 = 11;
197 constexpr int M41 = 12;
198 constexpr int M42 = 13;
199 constexpr int M43 = 14;
200 constexpr int M44 = 15;
201
202 /************************************************************************/
203 /* ==================================================================== */
204 /* CPG_STOKESRasterBand */
205 /* ==================================================================== */
206 /************************************************************************/
207
208 class CPG_STOKESRasterBand final: public GDALRasterBand
209 {
210 friend class CPGDataset;
211
212 int bNativeOrder;
213
214 public:
215 CPG_STOKESRasterBand( GDALDataset *poDS,
216 GDALDataType eType,
217 int bNativeOrder );
~CPG_STOKESRasterBand()218 ~CPG_STOKESRasterBand() override {}
219
220 CPLErr IReadBlock( int, int, void * ) override;
221 };
222
223 /************************************************************************/
224 /* AdjustFilename() */
225 /* */
226 /* Try to find the file with the request polarization and */
227 /* extension and update the passed filename accordingly. */
228 /* */
229 /* Return TRUE if file found otherwise FALSE. */
230 /************************************************************************/
231
AdjustFilename(char ** pszFilename,const char * pszPolarization,const char * pszExtension)232 int CPGDataset::AdjustFilename( char **pszFilename,
233 const char *pszPolarization,
234 const char *pszExtension )
235
236 {
237 // TODO: Eventually we should handle upper/lower case.
238 if ( EQUAL(pszPolarization,"stokes") )
239 {
240 const char *pszNewName =
241 CPLResetExtension(*pszFilename,
242 pszExtension);
243 CPLFree(*pszFilename);
244 *pszFilename = CPLStrdup(pszNewName);
245 }
246 else if (strlen(pszPolarization) == 2)
247 {
248 char *subptr = strstr(*pszFilename,"hh");
249 if (subptr == nullptr)
250 subptr = strstr(*pszFilename,"hv");
251 if (subptr == nullptr)
252 subptr = strstr(*pszFilename,"vv");
253 if (subptr == nullptr)
254 subptr = strstr(*pszFilename,"vh");
255 if (subptr == nullptr)
256 return FALSE;
257
258 strncpy( subptr, pszPolarization, 2);
259 const char *pszNewName =
260 CPLResetExtension(*pszFilename,
261 pszExtension);
262 CPLFree(*pszFilename);
263 *pszFilename = CPLStrdup(pszNewName);
264 }
265 else
266 {
267 const char *pszNewName =
268 CPLResetExtension(*pszFilename,
269 pszExtension);
270 CPLFree(*pszFilename);
271 *pszFilename = CPLStrdup(pszNewName);
272 }
273 VSIStatBufL sStatBuf;
274 return VSIStatL( *pszFilename, &sStatBuf ) == 0;
275 }
276
277 /************************************************************************/
278 /* Search for the various types of Convair filesets */
279 /* Return TRUE for a match, FALSE for no match */
280 /************************************************************************/
FindType1(const char * pszFilename)281 int CPGDataset::FindType1( const char *pszFilename )
282 {
283 const int nNameLen = static_cast<int>(strlen(pszFilename));
284
285 if ((strstr(pszFilename,"sso") == nullptr) &&
286 (strstr(pszFilename,"polgasp") == nullptr))
287 return FALSE;
288
289 if (( strlen(pszFilename) < 5) ||
290 (!EQUAL(pszFilename+nNameLen-4,".hdr")
291 && !EQUAL(pszFilename+nNameLen-4,".img")))
292 return FALSE;
293
294 /* Expect all bands and headers to be present */
295 char* pszTemp = CPLStrdup(pszFilename);
296
297 const bool bNotFound = !AdjustFilename( &pszTemp, "hh", "img" )
298 || !AdjustFilename( &pszTemp, "hh", "hdr" )
299 || !AdjustFilename( &pszTemp, "hv", "img" )
300 || !AdjustFilename( &pszTemp, "hv", "hdr" )
301 || !AdjustFilename( &pszTemp, "vh", "img" )
302 || !AdjustFilename( &pszTemp, "vh", "hdr" )
303 || !AdjustFilename( &pszTemp, "vv", "img" )
304 || !AdjustFilename( &pszTemp, "vv", "hdr" );
305
306 CPLFree(pszTemp);
307
308 return !bNotFound;
309 }
310
FindType2(const char * pszFilename)311 int CPGDataset::FindType2( const char *pszFilename )
312 {
313 const int nNameLen = static_cast<int>(strlen( pszFilename ));
314
315 if (( strlen(pszFilename) < 9) ||
316 (!EQUAL(pszFilename+nNameLen-8,"SIRC.hdr")
317 && !EQUAL(pszFilename+nNameLen-8,"SIRC.img")))
318 return FALSE;
319
320 char* pszTemp = CPLStrdup(pszFilename);
321 const bool bNotFound = !AdjustFilename( &pszTemp, "", "img" )
322 || !AdjustFilename( &pszTemp, "", "hdr" );
323 CPLFree(pszTemp);
324
325 return !bNotFound;
326 }
327
328 #ifdef notdef
FindType3(const char * pszFilename)329 int CPGDataset::FindType3( const char *pszFilename )
330 {
331 const int nNameLen = static_cast<int>(strlen( pszFilename ));
332
333 if ((strstr(pszFilename,"sso") == NULL) &&
334 (strstr(pszFilename,"polgasp") == NULL))
335 return FALSE;
336
337 if (( strlen(pszFilename) < 9) ||
338 (!EQUAL(pszFilename+nNameLen-4,".img")
339 && !EQUAL(pszFilename+nNameLen-8,".img_def")))
340 return FALSE;
341
342 char* pszTemp = CPLStrdup(pszFilename);
343 const bool bNotFound = !AdjustFilename( &pszTemp, "stokes", "img" )
344 || !AdjustFilename( &pszTemp, "stokes", "img_def" );
345 CPLFree(pszTemp);
346
347 return !bNotFound;
348 }
349 #endif
350
351 /************************************************************************/
352 /* LoadStokesLine() */
353 /************************************************************************/
354
LoadStokesLine(int iLine,int bNativeOrder)355 CPLErr CPGDataset::LoadStokesLine( int iLine, int bNativeOrder )
356
357 {
358 if( iLine == nLoadedStokesLine )
359 return CE_None;
360
361 const int nDataSize = GDALGetDataTypeSize(GDT_Float32) / 8;
362
363 /* -------------------------------------------------------------------- */
364 /* allocate working buffers if we don't have them already. */
365 /* -------------------------------------------------------------------- */
366 if( padfStokesMatrix == nullptr )
367 {
368 padfStokesMatrix = reinterpret_cast<float *>(
369 CPLMalloc( sizeof(float) * nRasterXSize * 16 ) );
370 }
371
372 /* -------------------------------------------------------------------- */
373 /* Load all the pixel data associated with this scanline. */
374 /* Retains same interleaving as original dataset. */
375 /* -------------------------------------------------------------------- */
376 if ( nInterleave == BIP )
377 {
378 const int offset = nRasterXSize * iLine * nDataSize * 16;
379 const int nBytesToRead = nDataSize * nRasterXSize*16;
380 if (( VSIFSeekL( afpImage[0], offset, SEEK_SET ) != 0 ) ||
381 static_cast<int>( VSIFReadL(
382 reinterpret_cast<GByte *>( padfStokesMatrix ),
383 1, nBytesToRead, afpImage[0] ) ) != nBytesToRead )
384 {
385 CPLError( CE_Failure, CPLE_FileIO,
386 "Error reading %d bytes of Stokes Convair at offset %d.\n"
387 "Reading file %s failed.",
388 nBytesToRead, offset, GetDescription() );
389 CPLFree( padfStokesMatrix );
390 padfStokesMatrix = nullptr;
391 nLoadedStokesLine = -1;
392 return CE_Failure;
393 }
394 }
395 else if ( nInterleave == BIL )
396 {
397 for ( int band_index = 0; band_index < 16; band_index++)
398 {
399 const int offset = nDataSize * (nRasterXSize * iLine +
400 nRasterXSize*band_index);
401 const int nBytesToRead = nDataSize * nRasterXSize;
402 if (( VSIFSeekL( afpImage[0], offset, SEEK_SET ) != 0 ) ||
403 static_cast<int>( VSIFReadL(
404 reinterpret_cast<GByte *>(
405 padfStokesMatrix + nBytesToRead*band_index ),
406 1, nBytesToRead, afpImage[0] ) ) != nBytesToRead )
407 {
408 CPLError( CE_Failure, CPLE_FileIO,
409 "Error reading %d bytes of Stokes Convair at offset %d.\n"
410 "Reading file %s failed.",
411 nBytesToRead, offset, GetDescription() );
412 CPLFree( padfStokesMatrix );
413 padfStokesMatrix = nullptr;
414 nLoadedStokesLine = -1;
415 return CE_Failure;
416 }
417 }
418 }
419 else
420 {
421 for ( int band_index = 0; band_index < 16; band_index++)
422 {
423 const int offset =
424 nDataSize * ( nRasterXSize * iLine +
425 nRasterXSize * nRasterYSize * band_index );
426 const int nBytesToRead = nDataSize * nRasterXSize;
427 if (( VSIFSeekL( afpImage[0], offset, SEEK_SET ) != 0 ) ||
428 static_cast<int>( VSIFReadL(
429 reinterpret_cast<GByte *>(
430 padfStokesMatrix + nBytesToRead * band_index ),
431 1, nBytesToRead, afpImage[0] ) ) != nBytesToRead )
432 {
433 CPLError( CE_Failure, CPLE_FileIO,
434 "Error reading %d bytes of Stokes Convair at offset %d.\n"
435 "Reading file %s failed.",
436 nBytesToRead, offset, GetDescription() );
437 CPLFree( padfStokesMatrix );
438 padfStokesMatrix = nullptr;
439 nLoadedStokesLine = -1;
440 return CE_Failure;
441 }
442 }
443 }
444
445 if (!bNativeOrder)
446 GDALSwapWords( padfStokesMatrix,nDataSize,nRasterXSize*16, nDataSize );
447
448 nLoadedStokesLine = iLine;
449
450 return CE_None;
451 }
452
453 /************************************************************************/
454 /* Parse header information and initialize dataset for the */
455 /* appropriate Convair dataset style. */
456 /* Returns dataset if successful; NULL if there was a problem. */
457 /************************************************************************/
458
InitializeType1Or2Dataset(const char * pszFilename)459 GDALDataset* CPGDataset::InitializeType1Or2Dataset( const char *pszFilename )
460 {
461
462 /* -------------------------------------------------------------------- */
463 /* Read the .hdr file (the hh one for the .sso and polgasp cases) */
464 /* and parse it. */
465 /* -------------------------------------------------------------------- */
466 int nLines = 0;
467 int nSamples = 0;
468 int nError = 0;
469
470 /* Parameters required for pseudo-geocoding. GCPs map */
471 /* slant range to ground range at 16 points. */
472 int iGeoParamsFound = 0;
473 int itransposed = 0;
474 double dfaltitude = 0.0;
475 double dfnear_srd = 0.0;
476 double dfsample_size = 0.0;
477 double dfsample_size_az = 0.0;
478
479 /* Parameters in geogratis geocoded images */
480 int iUTMParamsFound = 0;
481 int iUTMZone = 0;
482 // int iCorner = 0;
483 double dfnorth = 0.0;
484 double dfeast = 0.0;
485
486 char* pszWorkname = CPLStrdup(pszFilename);
487 AdjustFilename( &pszWorkname, "hh", "hdr" );
488 char **papszHdrLines = CSLLoad( pszWorkname );
489
490 for( int iLine = 0; papszHdrLines && papszHdrLines[iLine] != nullptr; iLine++ )
491 {
492 char **papszTokens = CSLTokenizeString( papszHdrLines[iLine] );
493
494 /* Note: some cv580 file seem to have comments with #, hence the >=
495 * instead of = for token checking, and the equalN for the corner.
496 */
497
498 if( CSLCount( papszTokens ) < 2 )
499 {
500 /* ignore */;
501 }
502 else if ( ( CSLCount( papszTokens ) >= 3 ) &&
503 EQUAL(papszTokens[0],"reference") &&
504 EQUAL(papszTokens[1],"north") )
505 {
506 dfnorth = CPLAtof(papszTokens[2]);
507 iUTMParamsFound++;
508 }
509 else if ( ( CSLCount( papszTokens ) >= 3 ) &&
510 EQUAL(papszTokens[0],"reference") &&
511 EQUAL(papszTokens[1],"east") )
512 {
513 dfeast = CPLAtof(papszTokens[2]);
514 iUTMParamsFound++;
515 }
516 else if ( ( CSLCount( papszTokens ) >= 5 ) &&
517 EQUAL(papszTokens[0],"reference") &&
518 EQUAL(papszTokens[1],"projection") &&
519 EQUAL(papszTokens[2],"UTM") &&
520 EQUAL(papszTokens[3],"zone") )
521 {
522 iUTMZone = atoi(papszTokens[4]);
523 iUTMParamsFound++;
524 }
525 else if ( ( CSLCount( papszTokens ) >= 3 ) &&
526 EQUAL(papszTokens[0],"reference") &&
527 EQUAL(papszTokens[1],"corner") &&
528 STARTS_WITH_CI(papszTokens[2], "Upper_Left") )
529 {
530 /* iCorner = 0; */
531 iUTMParamsFound++;
532 }
533 else if( EQUAL(papszTokens[0],"number_lines") )
534 nLines = atoi(papszTokens[1]);
535
536 else if( EQUAL(papszTokens[0],"number_samples") )
537 nSamples = atoi(papszTokens[1]);
538
539 else if( (EQUAL(papszTokens[0],"header_offset")
540 && atoi(papszTokens[1]) != 0)
541 || (EQUAL(papszTokens[0],"number_channels")
542 && (atoi(papszTokens[1]) != 1)
543 && (atoi(papszTokens[1]) != 10))
544 || (EQUAL(papszTokens[0],"datatype")
545 && atoi(papszTokens[1]) != 1)
546 || (EQUAL(papszTokens[0],"number_format")
547 && !EQUAL(papszTokens[1],"float32")
548 && !EQUAL(papszTokens[1],"int8")))
549 {
550 CPLError( CE_Failure, CPLE_AppDefined,
551 "Keyword %s has value %s which does not match CPG driver expectation.",
552 papszTokens[0], papszTokens[1] );
553 nError = 1;
554 }
555 else if( EQUAL(papszTokens[0],"altitude") )
556 {
557 dfaltitude = CPLAtof(papszTokens[1]);
558 iGeoParamsFound++;
559 }
560 else if( EQUAL(papszTokens[0],"near_srd") )
561 {
562 dfnear_srd = CPLAtof(papszTokens[1]);
563 iGeoParamsFound++;
564 }
565
566 else if( EQUAL(papszTokens[0],"sample_size") )
567 {
568 dfsample_size = CPLAtof(papszTokens[1]);
569 iGeoParamsFound++;
570 iUTMParamsFound++;
571 }
572 else if( EQUAL(papszTokens[0],"sample_size_az") )
573 {
574 dfsample_size_az = CPLAtof(papszTokens[1]);
575 iGeoParamsFound++;
576 iUTMParamsFound++;
577 }
578 else if( EQUAL(papszTokens[0],"transposed") )
579 {
580 itransposed = atoi(papszTokens[1]);
581 iGeoParamsFound++;
582 iUTMParamsFound++;
583 }
584
585 CSLDestroy( papszTokens );
586 }
587 CSLDestroy( papszHdrLines );
588 /* -------------------------------------------------------------------- */
589 /* Check for successful completion. */
590 /* -------------------------------------------------------------------- */
591 if( nError )
592 {
593 CPLFree(pszWorkname);
594 return nullptr;
595 }
596
597 if( nLines <= 0 || nSamples <= 0 )
598 {
599 CPLError( CE_Failure, CPLE_AppDefined,
600 "Did not find valid number_lines or number_samples keywords in %s.",
601 pszWorkname );
602 CPLFree(pszWorkname);
603 return nullptr;
604 }
605
606 /* -------------------------------------------------------------------- */
607 /* Initialize dataset. */
608 /* -------------------------------------------------------------------- */
609 CPGDataset *poDS = new CPGDataset();
610
611 poDS->nRasterXSize = nSamples;
612 poDS->nRasterYSize = nLines;
613
614 /* -------------------------------------------------------------------- */
615 /* Open the four bands. */
616 /* -------------------------------------------------------------------- */
617 static const char * const apszPolarizations[4] = { "hh", "hv", "vv", "vh" };
618
619 const int nNameLen = static_cast<int>(strlen(pszWorkname));
620
621 if ( EQUAL(pszWorkname+nNameLen-7,"IRC.hdr") ||
622 EQUAL(pszWorkname+nNameLen-7,"IRC.img") )
623 {
624
625 AdjustFilename( &pszWorkname, "" , "img" );
626 poDS->afpImage[0] = VSIFOpenL( pszWorkname, "rb" );
627 if( poDS->afpImage[0] == nullptr )
628 {
629 CPLError( CE_Failure, CPLE_OpenFailed,
630 "Failed to open .img file: %s",
631 pszWorkname );
632 CPLFree(pszWorkname);
633 delete poDS;
634 return nullptr;
635 }
636 poDS->aosImageFilenames.push_back(pszWorkname);
637 for( int iBand = 0; iBand < 4; iBand++ )
638 {
639 SIRC_QSLCRasterBand *poBand =
640 new SIRC_QSLCRasterBand( poDS, iBand+1, GDT_CFloat32 );
641 poDS->SetBand( iBand+1, poBand );
642 poBand->SetMetadataItem( "POLARIMETRIC_INTERP",
643 apszPolarizations[iBand] );
644 }
645 }
646 else
647 {
648 for( int iBand = 0; iBand < 4; iBand++ )
649 {
650 AdjustFilename( &pszWorkname, apszPolarizations[iBand], "img" );
651
652 poDS->afpImage[iBand] = VSIFOpenL( pszWorkname, "rb" );
653 if( poDS->afpImage[iBand] == nullptr )
654 {
655 CPLError( CE_Failure, CPLE_OpenFailed,
656 "Failed to open .img file: %s",
657 pszWorkname );
658 CPLFree(pszWorkname);
659 delete poDS;
660 return nullptr;
661 }
662 poDS->aosImageFilenames.push_back(pszWorkname);
663
664 RawRasterBand *poBand
665 = new RawRasterBand( poDS, iBand+1, poDS->afpImage[iBand],
666 0, 8, 8*nSamples,
667 GDT_CFloat32, !CPL_IS_LSB,
668 RawRasterBand::OwnFP::NO );
669 poDS->SetBand( iBand+1, poBand );
670
671 poBand->SetMetadataItem( "POLARIMETRIC_INTERP",
672 apszPolarizations[iBand] );
673 }
674 }
675
676 /* Set an appropriate matrix representation metadata item for the set */
677 if ( poDS->GetRasterCount() == 4 ) {
678 poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "SCATTERING" );
679 }
680
681 /* ------------------------------------------------------------------------- */
682 /* Add georeferencing or pseudo-geocoding, if enough information found. */
683 /* ------------------------------------------------------------------------- */
684 if (iUTMParamsFound == 7)
685 {
686 poDS->adfGeoTransform[1] = 0.0;
687 poDS->adfGeoTransform[2] = 0.0;
688 poDS->adfGeoTransform[4] = 0.0;
689 poDS->adfGeoTransform[5] = 0.0;
690
691 double dfnorth_center;
692 if (itransposed == 1)
693 {
694 CPLError(CE_Warning, CPLE_AppDefined,
695 "Did not have a convair SIRC-style test dataset\n"
696 "with transposed=1 for testing. Georeferencing may be "
697 "wrong.\n" );
698 dfnorth_center = dfnorth - nSamples*dfsample_size/2.0;
699 poDS->adfGeoTransform[0] = dfeast;
700 poDS->adfGeoTransform[2] = dfsample_size_az;
701 poDS->adfGeoTransform[3] = dfnorth;
702 poDS->adfGeoTransform[4] = -1*dfsample_size;
703 }
704 else
705 {
706 dfnorth_center = dfnorth - nLines*dfsample_size/2.0;
707 poDS->adfGeoTransform[0] = dfeast;
708 poDS->adfGeoTransform[1] = dfsample_size_az;
709 poDS->adfGeoTransform[3] = dfnorth;
710 poDS->adfGeoTransform[5] = -1*dfsample_size;
711 }
712
713 OGRSpatialReference oUTM;
714 if (dfnorth_center < 0)
715 oUTM.SetUTM(iUTMZone, 0);
716 else
717 oUTM.SetUTM(iUTMZone, 1);
718
719 /* Assuming WGS84 */
720 oUTM.SetWellKnownGeogCS( "WGS84" );
721 CPLFree( poDS->pszProjection );
722 poDS->pszProjection = nullptr;
723 oUTM.exportToWkt( &(poDS->pszProjection) );
724 }
725 else if (iGeoParamsFound == 5)
726 {
727 double dfgcpLine, dfgcpPixel, dfgcpX, dfgcpY, dftemp;
728
729 poDS->nGCPCount = 16;
730 poDS->pasGCPList = reinterpret_cast<GDAL_GCP *>(
731 CPLCalloc( sizeof(GDAL_GCP), poDS->nGCPCount ) );
732 GDALInitGCPs(poDS->nGCPCount, poDS->pasGCPList);
733
734 for( int ngcp = 0; ngcp < 16; ngcp ++ )
735 {
736 char szID[32];
737
738 snprintf( szID, sizeof(szID), "%d",ngcp+1);
739 if (itransposed == 1)
740 {
741 if (ngcp < 4)
742 dfgcpPixel = 0.0;
743 else if (ngcp < 8)
744 dfgcpPixel = nSamples/3.0;
745 else if (ngcp < 12)
746 dfgcpPixel = 2.0*nSamples/3.0;
747 else
748 dfgcpPixel = nSamples;
749
750 dfgcpLine = nLines*( ngcp % 4 )/3.0;
751
752 dftemp = dfnear_srd + (dfsample_size*dfgcpLine);
753 /* -1 so that 0,0 maps to largest Y */
754 dfgcpY = -1*sqrt( dftemp*dftemp - dfaltitude*dfaltitude );
755 dfgcpX = dfgcpPixel*dfsample_size_az;
756 }
757 else
758 {
759 if (ngcp < 4)
760 dfgcpLine = 0.0;
761 else if (ngcp < 8)
762 dfgcpLine = nLines/3.0;
763 else if (ngcp < 12)
764 dfgcpLine = 2.0*nLines/3.0;
765 else
766 dfgcpLine = nLines;
767
768 dfgcpPixel = nSamples*( ngcp % 4 )/3.0;
769
770 dftemp = dfnear_srd + (dfsample_size*dfgcpPixel);
771 dfgcpX = sqrt( dftemp*dftemp - dfaltitude*dfaltitude );
772 dfgcpY = (nLines - dfgcpLine)*dfsample_size_az;
773 }
774 poDS->pasGCPList[ngcp].dfGCPX = dfgcpX;
775 poDS->pasGCPList[ngcp].dfGCPY = dfgcpY;
776 poDS->pasGCPList[ngcp].dfGCPZ = 0.0;
777
778 poDS->pasGCPList[ngcp].dfGCPPixel = dfgcpPixel;
779 poDS->pasGCPList[ngcp].dfGCPLine = dfgcpLine;
780
781 CPLFree(poDS->pasGCPList[ngcp].pszId);
782 poDS->pasGCPList[ngcp].pszId = CPLStrdup( szID );
783 }
784
785 CPLFree(poDS->pszGCPProjection);
786 poDS->pszGCPProjection = CPLStrdup(
787 "LOCAL_CS[\"Ground range view / unreferenced meters\","
788 "UNIT[\"Meter\",1.0]]");
789 }
790
791 CPLFree(pszWorkname);
792
793 return poDS;
794 }
795
796 #ifdef notdef
InitializeType3Dataset(const char * pszFilename)797 GDALDataset *CPGDataset::InitializeType3Dataset( const char *pszFilename )
798 {
799 int iBytesPerPixel = 0;
800 int iInterleave = -1;
801 int nLines = 0;
802 int nSamples = 0;
803 int nBands = 0;
804 int nError = 0;
805
806 /* Parameters in geogratis geocoded images */
807 int iUTMParamsFound = 0;
808 int iUTMZone = 0;
809 double dfnorth = 0.0;
810 double dfeast = 0.0;
811 double dfOffsetX = 0.0;
812 double dfOffsetY = 0.0;
813 double dfxsize = 0.0;
814 double dfysize = 0.0;
815
816 char* pszWorkname = CPLStrdup(pszFilename);
817 AdjustFilename( &pszWorkname, "stokes", "img_def" );
818 char **papszHdrLines = CSLLoad( pszWorkname );
819
820 for( int iLine = 0; papszHdrLines && papszHdrLines[iLine] != NULL; iLine++ )
821 {
822 char **papszTokens
823 = CSLTokenizeString2( papszHdrLines[iLine], " \t",
824 CSLT_HONOURSTRINGS & CSLT_ALLOWEMPTYTOKENS );
825
826 /* Note: some cv580 file seem to have comments with #, hence the >=
827 * instead of = for token checking, and the equalN for the corner.
828 */
829
830 if ( ( CSLCount( papszTokens ) >= 3 ) &&
831 EQUAL(papszTokens[0],"data") &&
832 EQUAL(papszTokens[1],"organization:"))
833 {
834
835 if( STARTS_WITH_CI(papszTokens[2], "BSQ") )
836 iInterleave = BSQ;
837 else if( STARTS_WITH_CI(papszTokens[2], "BIL") )
838 iInterleave = BIL;
839 else if( STARTS_WITH_CI(papszTokens[2], "BIP") )
840 iInterleave = BIP;
841 else
842 {
843 CPLError( CE_Failure, CPLE_AppDefined,
844 "The interleaving type of the file (%s) is not supported.",
845 papszTokens[2] );
846 nError = 1;
847 }
848 }
849 else if ( ( CSLCount( papszTokens ) >= 3 ) &&
850 EQUAL(papszTokens[0],"data") &&
851 EQUAL(papszTokens[1],"state:") )
852 {
853
854 if( !STARTS_WITH_CI(papszTokens[2], "RAW") &&
855 !STARTS_WITH_CI(papszTokens[2], "GEO") )
856 {
857 CPLError( CE_Failure, CPLE_AppDefined,
858 "The data state of the file (%s) is not "
859 "supported.\n. Only RAW and GEO are currently "
860 "recognized.",
861 papszTokens[2] );
862 nError = 1;
863 }
864 }
865 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
866 EQUAL(papszTokens[0],"data") &&
867 EQUAL(papszTokens[1],"origin") &&
868 EQUAL(papszTokens[2],"point:") )
869 {
870 if (!STARTS_WITH_CI(papszTokens[3], "Upper_Left"))
871 {
872 CPLError( CE_Failure, CPLE_AppDefined,
873 "Unexpected value (%s) for data origin point- expect Upper_Left.",
874 papszTokens[3] );
875 nError = 1;
876 }
877 iUTMParamsFound++;
878 }
879 else if ( ( CSLCount( papszTokens ) >= 5 ) &&
880 EQUAL(papszTokens[0],"map") &&
881 EQUAL(papszTokens[1],"projection:") &&
882 EQUAL(papszTokens[2],"UTM") &&
883 EQUAL(papszTokens[3],"zone") )
884 {
885 iUTMZone = atoi(papszTokens[4]);
886 iUTMParamsFound++;
887 }
888 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
889 EQUAL(papszTokens[0],"project") &&
890 EQUAL(papszTokens[1],"origin:") )
891 {
892 dfeast = CPLAtof(papszTokens[2]);
893 dfnorth = CPLAtof(papszTokens[3]);
894 iUTMParamsFound+=2;
895 }
896 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
897 EQUAL(papszTokens[0],"file") &&
898 EQUAL(papszTokens[1],"start:"))
899 {
900 dfOffsetX = CPLAtof(papszTokens[2]);
901 dfOffsetY = CPLAtof(papszTokens[3]);
902 iUTMParamsFound+=2;
903 }
904 else if ( ( CSLCount( papszTokens ) >= 6 ) &&
905 EQUAL(papszTokens[0],"pixel") &&
906 EQUAL(papszTokens[1],"size") &&
907 EQUAL(papszTokens[2],"on") &&
908 EQUAL(papszTokens[3],"ground:"))
909 {
910 dfxsize = CPLAtof(papszTokens[4]);
911 dfysize = CPLAtof(papszTokens[5]);
912 iUTMParamsFound+=2;
913 }
914 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
915 EQUAL(papszTokens[0],"number") &&
916 EQUAL(papszTokens[1],"of") &&
917 EQUAL(papszTokens[2],"pixels:"))
918 {
919 nSamples = atoi(papszTokens[3]);
920 }
921 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
922 EQUAL(papszTokens[0],"number") &&
923 EQUAL(papszTokens[1],"of") &&
924 EQUAL(papszTokens[2],"lines:"))
925 {
926 nLines = atoi(papszTokens[3]);
927 }
928 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
929 EQUAL(papszTokens[0],"number") &&
930 EQUAL(papszTokens[1],"of") &&
931 EQUAL(papszTokens[2],"bands:"))
932 {
933 nBands = atoi(papszTokens[3]);
934 if ( nBands != 16)
935 {
936 CPLError( CE_Failure, CPLE_AppDefined,
937 "Number of bands has a value %s which does not match "
938 "CPG driver\nexpectation (expect a value of 16).",
939 papszTokens[3] );
940 nError = 1;
941 }
942 }
943 else if ( ( CSLCount( papszTokens ) >= 4 ) &&
944 EQUAL(papszTokens[0],"bytes") &&
945 EQUAL(papszTokens[1],"per") &&
946 EQUAL(papszTokens[2],"pixel:"))
947 {
948 iBytesPerPixel = atoi(papszTokens[3]);
949 if (iBytesPerPixel != 4)
950 {
951 CPLError( CE_Failure, CPLE_AppDefined,
952 "Bytes per pixel has a value %s which does not match "
953 "CPG driver\nexpectation (expect a value of 4).",
954 papszTokens[1] );
955 nError = 1;
956 }
957 }
958 CSLDestroy( papszTokens );
959 }
960
961 CSLDestroy( papszHdrLines );
962
963 /* -------------------------------------------------------------------- */
964 /* Check for successful completion. */
965 /* -------------------------------------------------------------------- */
966 if( nError )
967 {
968 CPLFree(pszWorkname);
969 return NULL;
970 }
971
972 if (!GDALCheckDatasetDimensions(nSamples, nLines) ||
973 !GDALCheckBandCount(nBands, FALSE) || iInterleave == -1 ||
974 iBytesPerPixel == 0 )
975 {
976 CPLError( CE_Failure, CPLE_AppDefined,
977 "%s is missing a required parameter (number of pixels, "
978 "number of lines,\number of bands, bytes per pixel, or "
979 "data organization).",
980 pszWorkname );
981 CPLFree(pszWorkname);
982 return NULL;
983 }
984
985 /* -------------------------------------------------------------------- */
986 /* Initialize dataset. */
987 /* -------------------------------------------------------------------- */
988 CPGDataset *poDS = new CPGDataset();
989
990 poDS->nRasterXSize = nSamples;
991 poDS->nRasterYSize = nLines;
992
993 if( iInterleave == BSQ )
994 poDS->nInterleave = BSQ;
995 else if( iInterleave == BIL )
996 poDS->nInterleave = BIL;
997 else
998 poDS->nInterleave = BIP;
999
1000 /* -------------------------------------------------------------------- */
1001 /* Open the 16 bands. */
1002 /* -------------------------------------------------------------------- */
1003
1004 AdjustFilename( &pszWorkname, "stokes" , "img" );
1005 poDS->afpImage[0] = VSIFOpenL( pszWorkname, "rb" );
1006 if( poDS->afpImage[0] == NULL )
1007 {
1008 CPLError( CE_Failure, CPLE_OpenFailed,
1009 "Failed to open .img file: %s",
1010 pszWorkname );
1011 CPLFree(pszWorkname);
1012 delete poDS;
1013 return NULL;
1014 }
1015 aosImageFilenames.push_back(pszWorkname);
1016 for( int iBand = 0; iBand < 16; iBand++ )
1017 {
1018 CPG_STOKESRasterBand *poBand
1019 = new CPG_STOKESRasterBand( poDS, GDT_CFloat32,
1020 !CPL_IS_LSB );
1021 poDS->SetBand( iBand+1, poBand );
1022 }
1023
1024 /* -------------------------------------------------------------------- */
1025 /* Set appropriate MATRIX_REPRESENTATION. */
1026 /* -------------------------------------------------------------------- */
1027 if ( poDS->GetRasterCount() == 6 ) {
1028 poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "COVARIANCE" );
1029 }
1030
1031 /* ------------------------------------------------------------------------- */
1032 /* Add georeferencing, if enough information found. */
1033 /* ------------------------------------------------------------------------- */
1034 if (iUTMParamsFound == 8)
1035 {
1036 double dfnorth_center = dfnorth - nLines*dfysize/2.0;
1037 poDS->adfGeoTransform[0] = dfeast + dfOffsetX;
1038 poDS->adfGeoTransform[1] = dfxsize;
1039 poDS->adfGeoTransform[2] = 0.0;
1040 poDS->adfGeoTransform[3] = dfnorth + dfOffsetY;
1041 poDS->adfGeoTransform[4] = 0.0;
1042 poDS->adfGeoTransform[5] = -1*dfysize;
1043
1044 OGRSpatialReference oUTM;
1045 if (dfnorth_center < 0)
1046 oUTM.SetUTM(iUTMZone, 0);
1047 else
1048 oUTM.SetUTM(iUTMZone, 1);
1049
1050 /* Assuming WGS84 */
1051 oUTM.SetWellKnownGeogCS( "WGS84" );
1052 CPLFree( poDS->pszProjection );
1053 poDS->pszProjection = NULL;
1054 oUTM.exportToWkt( &(poDS->pszProjection) );
1055 }
1056
1057 return poDS;
1058 }
1059 #endif
1060
1061 /************************************************************************/
1062 /* Open() */
1063 /************************************************************************/
1064
Open(GDALOpenInfo * poOpenInfo)1065 GDALDataset *CPGDataset::Open( GDALOpenInfo * poOpenInfo )
1066
1067 {
1068 /* -------------------------------------------------------------------- */
1069 /* Is this a PolGASP fileset? We expect fileset to follow */
1070 /* one of these patterns: */
1071 /* 1) <stuff>hh<stuff2>.img, <stuff>hh<stuff2>.hdr, */
1072 /* <stuff>hv<stuff2>.img, <stuff>hv<stuff2>.hdr, */
1073 /* <stuff>vh<stuff2>.img, <stuff>vh<stuff2>.hdr, */
1074 /* <stuff>vv<stuff2>.img, <stuff>vv<stuff2>.hdr, */
1075 /* where <stuff> or <stuff2> should contain the */
1076 /* substring "sso" or "polgasp" */
1077 /* 2) <stuff>SIRC.hdr and <stuff>SIRC.img */
1078 /* 3) <stuff>.img and <stuff>.img_def */
1079 /* where <stuff> should contain the */
1080 /* substring "sso" or "polgasp" */
1081 /* -------------------------------------------------------------------- */
1082 int CPGType = 0;
1083 if ( FindType1( poOpenInfo->pszFilename ))
1084 CPGType = 1;
1085 else if ( FindType2( poOpenInfo->pszFilename ))
1086 CPGType = 2;
1087
1088 /* Stokes matrix convair data: not quite working yet- something
1089 * is wrong in the interpretation of the matrix elements in terms
1090 * of hh, hv, vv, vh. Data will load if the next two lines are
1091 * uncommented, but values will be incorrect. Expect C11 = hh*conj(hh),
1092 * C12 = hh*conj(hv), etc. Used geogratis data in both scattering
1093 * matrix and stokes format for comparison.
1094 */
1095 //else if ( FindType3( poOpenInfo->pszFilename ))
1096 // CPGType = 3;
1097
1098 /* Set working name back to original */
1099
1100 if ( CPGType == 0 )
1101 {
1102 int nNameLen = static_cast<int>(strlen(poOpenInfo->pszFilename));
1103 if ( (nNameLen > 8) &&
1104 ( ( strstr(poOpenInfo->pszFilename,"sso") != nullptr ) ||
1105 ( strstr(poOpenInfo->pszFilename,"polgasp") != nullptr ) ) &&
1106 ( EQUAL(poOpenInfo->pszFilename+nNameLen-4,"img") ||
1107 EQUAL(poOpenInfo->pszFilename+nNameLen-4,"hdr") ||
1108 EQUAL(poOpenInfo->pszFilename+nNameLen-7,"img_def") ) )
1109 {
1110 CPLError( CE_Failure, CPLE_OpenFailed,
1111 "Apparent attempt to open Convair PolGASP data failed as\n"
1112 "one or more of the required files is missing (eight files\n"
1113 "are expected for scattering matrix format, two for Stokes)." );
1114 }
1115 else if ( (nNameLen > 8) &&
1116 ( strstr(poOpenInfo->pszFilename,"SIRC") != nullptr ) &&
1117 ( EQUAL(poOpenInfo->pszFilename+nNameLen-4,"img") ||
1118 EQUAL(poOpenInfo->pszFilename+nNameLen-4,"hdr")))
1119 {
1120 CPLError( CE_Failure, CPLE_OpenFailed,
1121 "Apparent attempt to open SIRC Convair PolGASP data failed \n"
1122 "as one of the expected files is missing (hdr or img)!" );
1123 }
1124 return nullptr;
1125 }
1126
1127 /* -------------------------------------------------------------------- */
1128 /* Confirm the requested access is supported. */
1129 /* -------------------------------------------------------------------- */
1130 if( poOpenInfo->eAccess == GA_Update )
1131 {
1132 CPLError( CE_Failure, CPLE_NotSupported,
1133 "The CPG driver does not support update access to existing"
1134 " datasets.\n" );
1135 return nullptr;
1136 }
1137
1138 /* Read the header info and create the dataset */
1139 #ifdef notdef
1140 if ( CPGType < 3 )
1141 #endif
1142 CPGDataset* poDS = reinterpret_cast<CPGDataset *>(
1143 InitializeType1Or2Dataset( poOpenInfo->pszFilename ) );
1144 #ifdef notdef
1145 else
1146 poDS = reinterpret_cast<CPGDataset *>(
1147 InitializeType3Dataset( poOpenInfo->pszFilename ) );
1148 #endif
1149 if( poDS == nullptr )
1150 return nullptr;
1151
1152 /* -------------------------------------------------------------------- */
1153 /* Check for overviews. */
1154 /* -------------------------------------------------------------------- */
1155 // Need to think about this.
1156 // poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1157
1158 /* -------------------------------------------------------------------- */
1159 /* Initialize any PAM information. */
1160 /* -------------------------------------------------------------------- */
1161 poDS->SetDescription( poOpenInfo->pszFilename );
1162 poDS->TryLoadXML();
1163
1164 return poDS;
1165 }
1166
1167 /************************************************************************/
1168 /* GetGCPCount() */
1169 /************************************************************************/
1170
GetGCPCount()1171 int CPGDataset::GetGCPCount()
1172
1173 {
1174 return nGCPCount;
1175 }
1176
1177 /************************************************************************/
1178 /* GetGCPProjection() */
1179 /************************************************************************/
1180
_GetGCPProjection()1181 const char *CPGDataset::_GetGCPProjection()
1182
1183 {
1184 return pszGCPProjection;
1185 }
1186
1187 /************************************************************************/
1188 /* GetGCPs() */
1189 /************************************************************************/
1190
GetGCPs()1191 const GDAL_GCP *CPGDataset::GetGCPs()
1192
1193 {
1194 return pasGCPList;
1195 }
1196
1197 /************************************************************************/
1198 /* GetProjectionRef() */
1199 /************************************************************************/
1200
_GetProjectionRef()1201 const char *CPGDataset::_GetProjectionRef()
1202
1203 {
1204 return pszProjection;
1205 }
1206
1207 /************************************************************************/
1208 /* GetGeoTransform() */
1209 /************************************************************************/
1210
GetGeoTransform(double * padfTransform)1211 CPLErr CPGDataset::GetGeoTransform( double * padfTransform )
1212
1213 {
1214 memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
1215 return CE_None;
1216 }
1217
1218 /************************************************************************/
1219 /* SIRC_QSLCRasterBand() */
1220 /************************************************************************/
1221
SIRC_QSLCRasterBand(CPGDataset * poGDSIn,int nBandIn,GDALDataType eType)1222 SIRC_QSLCRasterBand::SIRC_QSLCRasterBand( CPGDataset *poGDSIn, int nBandIn,
1223 GDALDataType eType )
1224
1225 {
1226 poDS = poGDSIn;
1227 nBand = nBandIn;
1228
1229 eDataType = eType;
1230
1231 nBlockXSize = poGDSIn->nRasterXSize;
1232 nBlockYSize = 1;
1233
1234 if( nBand == 1 )
1235 SetMetadataItem( "POLARIMETRIC_INTERP", "HH" );
1236 else if( nBand == 2 )
1237 SetMetadataItem( "POLARIMETRIC_INTERP", "HV" );
1238 else if( nBand == 3 )
1239 SetMetadataItem( "POLARIMETRIC_INTERP", "VH" );
1240 else if( nBand == 4 )
1241 SetMetadataItem( "POLARIMETRIC_INTERP", "VV" );
1242 }
1243
1244 /************************************************************************/
1245 /* IReadBlock() */
1246 /************************************************************************/
1247
1248 /* From: http://southport.jpl.nasa.gov/software/dcomp/dcomp.html
1249
1250 ysca = sqrt{ [ (Byte(2) / 254 ) + 1.5] 2Byte(1) }
1251 Re(SHH) = byte(3) ysca/127
1252 Im(SHH) = byte(4) ysca/127
1253 Re(SHV) = byte(5) ysca/127
1254 Im(SHV) = byte(6) ysca/127
1255 Re(SVH) = byte(7) ysca/127
1256 Im(SVH) = byte(8) ysca/127
1257 Re(SVV) = byte(9) ysca/127
1258 Im(SVV) = byte(10) ysca/127
1259 */
1260
IReadBlock(CPL_UNUSED int nBlockXOff,int nBlockYOff,void * pImage)1261 CPLErr SIRC_QSLCRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff,
1262 int nBlockYOff,
1263 void * pImage )
1264 {
1265 const int nBytesPerSample = 10;
1266 CPGDataset *poGDS = reinterpret_cast<CPGDataset *>( poDS );
1267 const int offset = nBlockXSize* nBlockYOff*nBytesPerSample;
1268
1269 /* -------------------------------------------------------------------- */
1270 /* Load all the pixel data associated with this scanline. */
1271 /* -------------------------------------------------------------------- */
1272 const int nBytesToRead = nBytesPerSample * nBlockXSize;
1273
1274 GByte *pabyRecord = reinterpret_cast<GByte *>(
1275 CPLMalloc( nBytesToRead ) );
1276
1277 if( VSIFSeekL( poGDS->afpImage[0], offset, SEEK_SET ) != 0
1278 || static_cast<int>( VSIFReadL(
1279 pabyRecord, 1, nBytesToRead, poGDS->afpImage[0] ) )
1280 != nBytesToRead )
1281 {
1282 CPLError( CE_Failure, CPLE_FileIO,
1283 "Error reading %d bytes of SIRC Convair at offset %d.\n"
1284 "Reading file %s failed.",
1285 nBytesToRead, offset, poGDS->GetDescription() );
1286 CPLFree( pabyRecord );
1287 return CE_Failure;
1288 }
1289
1290 /* -------------------------------------------------------------------- */
1291 /* Initialize our power table if this is our first time through. */
1292 /* -------------------------------------------------------------------- */
1293 static float afPowTable[256];
1294 static bool bPowTableInitialized = false;
1295
1296 if( !bPowTableInitialized )
1297 {
1298 bPowTableInitialized = true;
1299
1300 for( int i = 0; i < 256; i++ )
1301 {
1302 afPowTable[i] = static_cast<float>( pow( 2.0, i-128 ) );
1303 }
1304 }
1305
1306 /* -------------------------------------------------------------------- */
1307 /* Copy the desired band out based on the size of the type, and */
1308 /* the interleaving mode. */
1309 /* -------------------------------------------------------------------- */
1310 for( int iX = 0; iX < nBlockXSize; iX++ )
1311 {
1312 unsigned char *pabyGroup = pabyRecord + iX * nBytesPerSample;
1313 const signed char *Byte = reinterpret_cast<signed char *>(
1314 pabyGroup - 1 ); /* A ones based alias */
1315
1316 /* coverity[tainted_data] */
1317 const double dfScale
1318 = sqrt( (static_cast<double>(Byte[2]) / 254 + 1.5) * afPowTable[Byte[1] + 128] );
1319
1320 float *pafImage = reinterpret_cast<float *>( pImage );
1321
1322 if( nBand == 1 )
1323 {
1324 const float fReSHH = static_cast<float>(Byte[3] * dfScale / 127.0);
1325 const float fImSHH = static_cast<float>(Byte[4] * dfScale / 127.0);
1326
1327 pafImage[iX*2 ] = fReSHH;
1328 pafImage[iX*2+1] = fImSHH;
1329 }
1330 else if( nBand == 2 )
1331 {
1332 const float fReSHV = static_cast<float>(Byte[5] * dfScale / 127.0);
1333 const float fImSHV = static_cast<float>(Byte[6] * dfScale / 127.0);
1334
1335 pafImage[iX*2 ] = fReSHV;
1336 pafImage[iX*2+1] = fImSHV;
1337 }
1338 else if( nBand == 3 )
1339 {
1340 const float fReSVH = static_cast<float>(Byte[7] * dfScale / 127.0);
1341 const float fImSVH = static_cast<float>(Byte[8] * dfScale / 127.0);
1342
1343 pafImage[iX*2 ] = fReSVH;
1344 pafImage[iX*2+1] = fImSVH;
1345 }
1346 else if( nBand == 4 )
1347 {
1348 const float fReSVV = static_cast<float>(Byte[9] * dfScale / 127.0);
1349 const float fImSVV = static_cast<float>(Byte[10]* dfScale / 127.0);
1350
1351 pafImage[iX*2 ] = fReSVV;
1352 pafImage[iX*2+1] = fImSVV;
1353 }
1354 }
1355
1356 CPLFree( pabyRecord );
1357
1358 return CE_None;
1359 }
1360
1361 /************************************************************************/
1362 /* CPG_STOKESRasterBand() */
1363 /************************************************************************/
1364
CPG_STOKESRasterBand(GDALDataset * poDSIn,GDALDataType eType,int bNativeOrderIn)1365 CPG_STOKESRasterBand::CPG_STOKESRasterBand( GDALDataset *poDSIn,
1366 GDALDataType eType,
1367 int bNativeOrderIn ) :
1368 bNativeOrder(bNativeOrderIn)
1369 {
1370 static const char * const apszPolarizations[16] = {
1371 "Covariance_11",
1372 "Covariance_12",
1373 "Covariance_13",
1374 "Covariance_14",
1375 "Covariance_21",
1376 "Covariance_22",
1377 "Covariance_23",
1378 "Covariance_24",
1379 "Covariance_31",
1380 "Covariance_32",
1381 "Covariance_33",
1382 "Covariance_34",
1383 "Covariance_41",
1384 "Covariance_42",
1385 "Covariance_43",
1386 "Covariance_44" };
1387
1388 poDS = poDSIn;
1389 eDataType = eType;
1390
1391 nBlockXSize = poDSIn->GetRasterXSize();
1392 nBlockYSize = 1;
1393
1394 SetMetadataItem( "POLARIMETRIC_INTERP", apszPolarizations[nBand-1] );
1395 SetDescription( apszPolarizations[nBand-1] );
1396 }
1397
1398 /************************************************************************/
1399 /* IReadBlock() */
1400 /************************************************************************/
1401
1402 /* Convert from Stokes to Covariance representation */
1403
IReadBlock(CPL_UNUSED int nBlockXOff,int nBlockYOff,void * pImage)1404 CPLErr CPG_STOKESRasterBand::IReadBlock( CPL_UNUSED int nBlockXOff,
1405 int nBlockYOff,
1406 void * pImage )
1407
1408 {
1409 CPLAssert( nBlockXOff == 0 );
1410
1411 int m11, /* m12, */ m13, m14, /* m21, */ m22, m23, m24, step;
1412 int m31, m32, m33, m34, m41, m42, m43, m44;
1413 CPGDataset *poGDS = reinterpret_cast<CPGDataset *>( poDS );
1414
1415 CPLErr eErr = poGDS->LoadStokesLine(nBlockYOff, bNativeOrder);
1416 if( eErr != CE_None )
1417 return eErr;
1418
1419 float *M = poGDS->padfStokesMatrix;
1420 float *pafLine = reinterpret_cast<float *>( pImage );
1421
1422 if ( poGDS->nInterleave == BIP)
1423 {
1424 step = 16;
1425 m11 = M11;
1426 // m12 = M12;
1427 m13 = M13;
1428 m14 = M14;
1429 // m21 = M21;
1430 m22 = M22;
1431 m23 = M23;
1432 m24 = M24;
1433 m31 = M31;
1434 m32 = M32;
1435 m33 = M33;
1436 m34 = M34;
1437 m41 = M41;
1438 m42 = M42;
1439 m43 = M43;
1440 m44 = M44;
1441 }
1442 else
1443 {
1444 step = 1;
1445 m11=0;
1446 // m12=nRasterXSize;
1447 m13=nRasterXSize*2;
1448 m14=nRasterXSize*3;
1449 // m21=nRasterXSize*4;
1450 m22=nRasterXSize*5;
1451 m23=nRasterXSize*6;
1452 m24=nRasterXSize*7;
1453 m31=nRasterXSize*8;
1454 m32=nRasterXSize*9;
1455 m33=nRasterXSize*10;
1456 m34=nRasterXSize*11;
1457 m41=nRasterXSize*12;
1458 m42=nRasterXSize*13;
1459 m43=nRasterXSize*14;
1460 m44=nRasterXSize*15;
1461 }
1462 if ( nBand == 1 ) /* C11 */
1463 {
1464 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1465 {
1466 pafLine[iPixel*2+0] = M[m11]-M[m22]-M[m33]+M[m44];
1467 pafLine[iPixel*2+1] = 0.0;
1468 m11 += step;
1469 m22 += step;
1470 m33 += step;
1471 m44 += step;
1472 }
1473 }
1474 else if ( nBand == 2 ) /* C12 */
1475 {
1476 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1477 {
1478 pafLine[iPixel*2+0] = M[m13]-M[m23];
1479 pafLine[iPixel*2+1] = M[m14]-M[m24];
1480 m13 += step;
1481 m23 += step;
1482 m14 += step;
1483 m24 += step;
1484 }
1485 }
1486 else if ( nBand == 3 ) /* C13 */
1487 {
1488 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1489 {
1490 pafLine[iPixel*2+0] = M[m33]-M[m44];
1491 pafLine[iPixel*2+1] = M[m43]+M[m34];
1492 m33 += step;
1493 m44 += step;
1494 m43 += step;
1495 m34 += step;
1496 }
1497 }
1498 else if ( nBand == 4 ) /* C14 */
1499 {
1500 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1501 {
1502 pafLine[iPixel*2+0] = M[m31]-M[m32];
1503 pafLine[iPixel*2+1] = M[m41]-M[m42];
1504 m31 += step;
1505 m32 += step;
1506 m41 += step;
1507 m42 += step;
1508 }
1509 }
1510 else if ( nBand == 5 ) /* C21 */
1511 {
1512 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1513 {
1514 pafLine[iPixel*2+0] = M[m13]-M[m23];
1515 pafLine[iPixel*2+1] = M[m24]-M[m14];
1516 m13 += step;
1517 m23 += step;
1518 m14 += step;
1519 m24 += step;
1520 }
1521 }
1522 else if ( nBand == 6 ) /* C22 */
1523 {
1524 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1525 {
1526 pafLine[iPixel*2+0] = M[m11]+M[m22]-M[m33]-M[m44];
1527 pafLine[iPixel*2+1] = 0.0;
1528 m11 += step;
1529 m22 += step;
1530 m33 += step;
1531 m44 += step;
1532 }
1533 }
1534 else if ( nBand == 7 ) /* C23 */
1535 {
1536 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1537 {
1538 pafLine[iPixel*2+0] = M[m31]+M[m32];
1539 pafLine[iPixel*2+1] = M[m41]+M[m42];
1540 m31 += step;
1541 m32 += step;
1542 m41 += step;
1543 m42 += step;
1544 }
1545 }
1546 else if ( nBand == 8 ) /* C24 */
1547 {
1548 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1549 {
1550 pafLine[iPixel*2+0] = M[m33]+M[m44];
1551 pafLine[iPixel*2+1] = M[m43]-M[m34];
1552 m33 += step;
1553 m44 += step;
1554 m43 += step;
1555 m34 += step;
1556 }
1557 }
1558 else if ( nBand == 9 ) /* C31 */
1559 {
1560 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1561 {
1562 pafLine[iPixel*2+0] = M[m33]-M[m44];
1563 pafLine[iPixel*2+1] = -1*M[m43]-M[m34];
1564 m33 += step;
1565 m44 += step;
1566 m43 += step;
1567 m34 += step;
1568 }
1569 }
1570 else if ( nBand == 10 ) /* C32 */
1571 {
1572 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1573 {
1574 pafLine[iPixel*2+0] = M[m31]+M[m32];
1575 pafLine[iPixel*2+1] = -1*M[m41]-M[m42];
1576 m31 += step;
1577 m32 += step;
1578 m41 += step;
1579 m42 += step;
1580 }
1581 }
1582 else if ( nBand == 11 ) /* C33 */
1583 {
1584 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1585 {
1586 pafLine[iPixel*2+0] = M[m11]+M[m22]+M[m33]+M[m44];
1587 pafLine[iPixel*2+1] = 0.0;
1588 m11 += step;
1589 m22 += step;
1590 m33 += step;
1591 m44 += step;
1592 }
1593 }
1594 else if ( nBand == 12 ) /* C34 */
1595 {
1596 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1597 {
1598 pafLine[iPixel*2+0] = M[m13]-M[m23];
1599 pafLine[iPixel*2+1] = -1*M[m14]-M[m24];
1600 m13 += step;
1601 m23 += step;
1602 m14 += step;
1603 m24 += step;
1604 }
1605 }
1606 else if ( nBand == 13 ) /* C41 */
1607 {
1608 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1609 {
1610 pafLine[iPixel*2+0] = M[m31]-M[m32];
1611 pafLine[iPixel*2+1] = M[m42]-M[m41];
1612 m31 += step;
1613 m32 += step;
1614 m41 += step;
1615 m42 += step;
1616 }
1617 }
1618 else if ( nBand == 14 ) /* C42 */
1619 {
1620 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1621 {
1622 pafLine[iPixel*2+0] = M[m33]+M[m44];
1623 pafLine[iPixel*2+1] = M[m34]-M[m43];
1624 m33 += step;
1625 m44 += step;
1626 m43 += step;
1627 m34 += step;
1628 }
1629 }
1630 else if ( nBand == 15 ) /* C43 */
1631 {
1632 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1633 {
1634 pafLine[iPixel*2+0] = M[m13]-M[m23];
1635 pafLine[iPixel*2+1] = M[m14]+M[m24];
1636 m13 += step;
1637 m23 += step;
1638 m14 += step;
1639 m24 += step;
1640 }
1641 }
1642 else /* C44 */
1643 {
1644 for ( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1645 {
1646 pafLine[iPixel*2+0] = M[m11]-M[m22]+M[m33]-M[m44];
1647 pafLine[iPixel*2+1] = 0.0;
1648 m11 += step;
1649 m22 += step;
1650 m33 += step;
1651 m44 += step;
1652 }
1653 }
1654
1655 return CE_None;
1656 }
1657
1658 /************************************************************************/
1659 /* GDALRegister_CPG() */
1660 /************************************************************************/
1661
GDALRegister_CPG()1662 void GDALRegister_CPG()
1663
1664 {
1665 if( GDALGetDriverByName( "CPG" ) != nullptr )
1666 return;
1667
1668 GDALDriver *poDriver = new GDALDriver();
1669
1670 poDriver->SetDescription( "CPG" );
1671 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
1672 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Convair PolGASP" );
1673 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1674
1675 poDriver->pfnOpen = CPGDataset::Open;
1676
1677 GetGDALDriverManager()->RegisterDriver( poDriver );
1678 }
1679