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