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