Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Dataset storage of geolocation array and backmap
5 : * Author: Even Rouault, <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2022, Planet Labs
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalcachedpixelaccessor.h"
14 :
15 : /*! @cond Doxygen_Suppress */
16 :
17 : /************************************************************************/
18 : /* GDALGeoLocDatasetAccessors */
19 : /************************************************************************/
20 :
21 : class GDALGeoLocDatasetAccessors
22 : {
23 : typedef class GDALGeoLocDatasetAccessors AccessorType;
24 :
25 : GDALGeoLocTransformInfo *m_psTransform;
26 :
27 : CPLStringList m_aosGTiffCreationOptions{};
28 :
29 : GDALDataset *m_poGeolocTmpDataset = nullptr;
30 : GDALDataset *m_poBackmapTmpDataset = nullptr;
31 : GDALDataset *m_poBackmapWeightsTmpDataset = nullptr;
32 :
33 : GDALGeoLocDatasetAccessors(const GDALGeoLocDatasetAccessors &) = delete;
34 : GDALGeoLocDatasetAccessors &
35 : operator=(const GDALGeoLocDatasetAccessors &) = delete;
36 :
37 : bool LoadGeoloc(bool bIsRegularGrid);
38 :
39 : public:
40 : static constexpr int TILE_SIZE = 256;
41 : static constexpr int TILE_COUNT = 64;
42 :
43 : GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocXAccessor;
44 : GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocYAccessor;
45 : GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapXAccessor;
46 : GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapYAccessor;
47 : GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapWeightAccessor;
48 :
49 2 : explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo *psTransform)
50 2 : : m_psTransform(psTransform), geolocXAccessor(nullptr),
51 : geolocYAccessor(nullptr), backMapXAccessor(nullptr),
52 2 : backMapYAccessor(nullptr), backMapWeightAccessor(nullptr)
53 : {
54 2 : m_aosGTiffCreationOptions.SetNameValue("TILED", "YES");
55 2 : m_aosGTiffCreationOptions.SetNameValue("INTERLEAVE", "BAND");
56 : m_aosGTiffCreationOptions.SetNameValue("BLOCKXSIZE",
57 2 : CPLSPrintf("%d", TILE_SIZE));
58 : m_aosGTiffCreationOptions.SetNameValue("BLOCKYSIZE",
59 2 : CPLSPrintf("%d", TILE_SIZE));
60 2 : }
61 :
62 : ~GDALGeoLocDatasetAccessors();
63 :
64 : bool Load(bool bIsRegularGrid, bool bUseQuadtree);
65 :
66 : bool AllocateBackMap();
67 :
68 : GDALDataset *GetBackmapDataset();
69 : void FlushBackmapCaches();
70 :
71 2 : static void ReleaseBackmapDataset(GDALDataset *)
72 : {
73 2 : }
74 :
75 : void FreeWghtsBackMap();
76 : };
77 :
78 : /************************************************************************/
79 : /* ~GDALGeoLocDatasetAccessors() */
80 : /************************************************************************/
81 :
82 2 : GDALGeoLocDatasetAccessors::~GDALGeoLocDatasetAccessors()
83 : {
84 2 : geolocXAccessor.ResetModifiedFlag();
85 2 : geolocYAccessor.ResetModifiedFlag();
86 2 : backMapXAccessor.ResetModifiedFlag();
87 2 : backMapYAccessor.ResetModifiedFlag();
88 :
89 2 : FreeWghtsBackMap();
90 :
91 2 : delete m_poGeolocTmpDataset;
92 2 : delete m_poBackmapTmpDataset;
93 2 : }
94 :
95 : /************************************************************************/
96 : /* AllocateBackMap() */
97 : /************************************************************************/
98 :
99 2 : bool GDALGeoLocDatasetAccessors::AllocateBackMap()
100 : {
101 2 : auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
102 2 : if (poDriver == nullptr)
103 0 : return false;
104 :
105 : // CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
106 : // so store them in a long-lived std::string
107 : const std::string osBackmapTmpFilename =
108 4 : CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif");
109 2 : m_poBackmapTmpDataset = poDriver->Create(
110 2 : osBackmapTmpFilename.c_str(), m_psTransform->nBackMapWidth,
111 2 : m_psTransform->nBackMapHeight, 2, GDT_Float32,
112 2 : m_aosGTiffCreationOptions.List());
113 2 : if (m_poBackmapTmpDataset == nullptr)
114 : {
115 0 : return false;
116 : }
117 2 : m_poBackmapTmpDataset->MarkSuppressOnClose();
118 2 : VSIUnlink(m_poBackmapTmpDataset->GetDescription());
119 2 : auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
120 2 : auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
121 :
122 2 : backMapXAccessor.SetBand(poBandX);
123 2 : backMapYAccessor.SetBand(poBandY);
124 :
125 : // CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
126 : // so store them in a long-lived std::string
127 : const std::string osBackmapWeightsTmpFilename =
128 4 : CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif");
129 2 : m_poBackmapWeightsTmpDataset = poDriver->Create(
130 2 : osBackmapWeightsTmpFilename.c_str(), m_psTransform->nBackMapWidth,
131 2 : m_psTransform->nBackMapHeight, 1, GDT_Float32,
132 2 : m_aosGTiffCreationOptions.List());
133 2 : if (m_poBackmapWeightsTmpDataset == nullptr)
134 : {
135 0 : return false;
136 : }
137 2 : m_poBackmapWeightsTmpDataset->MarkSuppressOnClose();
138 2 : VSIUnlink(m_poBackmapWeightsTmpDataset->GetDescription());
139 2 : backMapWeightAccessor.SetBand(
140 2 : m_poBackmapWeightsTmpDataset->GetRasterBand(1));
141 :
142 2 : return true;
143 : }
144 :
145 : /************************************************************************/
146 : /* FreeWghtsBackMap() */
147 : /************************************************************************/
148 :
149 4 : void GDALGeoLocDatasetAccessors::FreeWghtsBackMap()
150 : {
151 4 : if (m_poBackmapWeightsTmpDataset)
152 : {
153 2 : backMapWeightAccessor.ResetModifiedFlag();
154 2 : delete m_poBackmapWeightsTmpDataset;
155 2 : m_poBackmapWeightsTmpDataset = nullptr;
156 : }
157 4 : }
158 :
159 : /************************************************************************/
160 : /* GetBackmapDataset() */
161 : /************************************************************************/
162 :
163 2 : GDALDataset *GDALGeoLocDatasetAccessors::GetBackmapDataset()
164 : {
165 2 : auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
166 2 : auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
167 2 : poBandX->SetNoDataValue(INVALID_BMXY);
168 2 : poBandY->SetNoDataValue(INVALID_BMXY);
169 2 : return m_poBackmapTmpDataset;
170 : }
171 :
172 : /************************************************************************/
173 : /* FlushBackmapCaches() */
174 : /************************************************************************/
175 :
176 2 : void GDALGeoLocDatasetAccessors::FlushBackmapCaches()
177 : {
178 2 : backMapXAccessor.FlushCache();
179 2 : backMapYAccessor.FlushCache();
180 2 : }
181 :
182 : /************************************************************************/
183 : /* Load() */
184 : /************************************************************************/
185 :
186 2 : bool GDALGeoLocDatasetAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
187 : {
188 4 : return LoadGeoloc(bIsRegularGrid) &&
189 0 : ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
190 2 : (!bUseQuadtree &&
191 4 : GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
192 : }
193 :
194 : /************************************************************************/
195 : /* LoadGeoloc() */
196 : /************************************************************************/
197 :
198 2 : bool GDALGeoLocDatasetAccessors::LoadGeoloc(bool bIsRegularGrid)
199 :
200 : {
201 2 : if (bIsRegularGrid)
202 : {
203 1 : const int nXSize = m_psTransform->nGeoLocXSize;
204 1 : const int nYSize = m_psTransform->nGeoLocYSize;
205 :
206 1 : auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
207 1 : if (poDriver == nullptr)
208 0 : return false;
209 :
210 : // CPLResetExtension / CPLGenerateTempFilename generate short-lived
211 : // strings, so store them in a long-lived std::string
212 : const std::string osGeolocTmpFilename =
213 1 : CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif");
214 1 : m_poGeolocTmpDataset =
215 1 : poDriver->Create(osGeolocTmpFilename.c_str(), nXSize, nYSize, 2,
216 1 : GDT_Float64, m_aosGTiffCreationOptions.List());
217 1 : if (m_poGeolocTmpDataset == nullptr)
218 : {
219 0 : return false;
220 : }
221 1 : m_poGeolocTmpDataset->MarkSuppressOnClose();
222 1 : VSIUnlink(m_poGeolocTmpDataset->GetDescription());
223 :
224 1 : auto poXBand = m_poGeolocTmpDataset->GetRasterBand(1);
225 1 : auto poYBand = m_poGeolocTmpDataset->GetRasterBand(2);
226 :
227 : // Case of regular grid.
228 : // The XBAND contains the x coordinates for all lines.
229 : // The YBAND contains the y coordinates for all columns.
230 :
231 : double *padfTempX =
232 1 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
233 : double *padfTempY =
234 1 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
235 1 : if (padfTempX == nullptr || padfTempY == nullptr)
236 : {
237 0 : CPLFree(padfTempX);
238 0 : CPLFree(padfTempY);
239 0 : return false;
240 : }
241 :
242 : CPLErr eErr =
243 1 : GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
244 : padfTempX, nXSize, 1, GDT_Float64, 0, 0);
245 :
246 81 : for (int j = 0; j < nYSize; j++)
247 : {
248 80 : if (poXBand->RasterIO(GF_Write, 0, j, nXSize, 1, padfTempX, nXSize,
249 80 : 1, GDT_Float64, 0, 0, nullptr) != CE_None)
250 : {
251 0 : eErr = CE_Failure;
252 0 : break;
253 : }
254 : }
255 :
256 1 : if (eErr == CE_None)
257 : {
258 1 : eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
259 : 1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
260 :
261 361 : for (int i = 0; i < nXSize; i++)
262 : {
263 360 : if (poYBand->RasterIO(GF_Write, i, 0, 1, nYSize, padfTempY, 1,
264 : nYSize, GDT_Float64, 0, 0,
265 360 : nullptr) != CE_None)
266 : {
267 0 : eErr = CE_Failure;
268 0 : break;
269 : }
270 : }
271 : }
272 :
273 1 : CPLFree(padfTempX);
274 1 : CPLFree(padfTempY);
275 :
276 1 : if (eErr != CE_None)
277 0 : return false;
278 :
279 1 : geolocXAccessor.SetBand(poXBand);
280 1 : geolocYAccessor.SetBand(poYBand);
281 : }
282 : else
283 : {
284 1 : geolocXAccessor.SetBand(
285 1 : GDALRasterBand::FromHandle(m_psTransform->hBand_X));
286 1 : geolocYAccessor.SetBand(
287 1 : GDALRasterBand::FromHandle(m_psTransform->hBand_Y));
288 : }
289 :
290 2 : GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(m_psTransform);
291 2 : return true;
292 : }
293 :
294 : /*! @endcond */
|
