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