GDAL

gdalgeoloc_dataset_accessor.h
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 * Permission is hereby granted, free of charge, to any person obtaining a
11 * copy of this software and associated documentation files (the "Software"),
12 * to deal in the Software without restriction, including without limitation
13 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 * and/or sell copies of the Software, and to permit persons to whom the
15 * Software is furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included
18 * in all copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 * DEALINGS IN THE SOFTWARE.
27 ****************************************************************************/
28
29#include "gdalcachedpixelaccessor.h"
30
33/************************************************************************/
34/* GDALGeoLocDatasetAccessors */
35/************************************************************************/
36
37 class GDALGeoLocDatasetAccessors
38{
39 typedef class GDALGeoLocDatasetAccessors AccessorType;
40
41 GDALGeoLocTransformInfo *m_psTransform;
42
43 CPLStringList m_aosGTiffCreationOptions{};
44
45 GDALDataset *m_poGeolocTmpDataset = nullptr;
46 GDALDataset *m_poBackmapTmpDataset = nullptr;
47 GDALDataset *m_poBackmapWeightsTmpDataset = nullptr;
48
49 GDALGeoLocDatasetAccessors(const GDALGeoLocDatasetAccessors &) = delete;
50 GDALGeoLocDatasetAccessors &
51 operator=(const GDALGeoLocDatasetAccessors &) = delete;
52
53 bool LoadGeoloc(bool bIsRegularGrid);
54
55 public:
56 static constexpr int TILE_SIZE = 256;
57 static constexpr int TILE_COUNT = 64;
58
59 GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocXAccessor;
60 GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocYAccessor;
61 GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapXAccessor;
62 GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapYAccessor;
63 GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapWeightAccessor;
64
65 explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo *psTransform)
66 : m_psTransform(psTransform), geolocXAccessor(nullptr),
67 geolocYAccessor(nullptr), backMapXAccessor(nullptr),
68 backMapYAccessor(nullptr), backMapWeightAccessor(nullptr)
69 {
70 m_aosGTiffCreationOptions.SetNameValue("TILED", "YES");
71 m_aosGTiffCreationOptions.SetNameValue("INTERLEAVE", "BAND");
72 m_aosGTiffCreationOptions.SetNameValue("BLOCKXSIZE",
73 CPLSPrintf("%d", TILE_SIZE));
74 m_aosGTiffCreationOptions.SetNameValue("BLOCKYSIZE",
75 CPLSPrintf("%d", TILE_SIZE));
76 }
77
78 ~GDALGeoLocDatasetAccessors();
79
80 bool Load(bool bIsRegularGrid, bool bUseQuadtree);
81
82 bool AllocateBackMap();
83
84 GDALDataset *GetBackmapDataset();
85 void FlushBackmapCaches();
86
87 static void ReleaseBackmapDataset(GDALDataset *)
88 {
89 }
90
91 void FreeWghtsBackMap();
92};
93
94/************************************************************************/
95/* ~GDALGeoLocDatasetAccessors() */
96/************************************************************************/
97
98GDALGeoLocDatasetAccessors::~GDALGeoLocDatasetAccessors()
99{
100 geolocXAccessor.ResetModifiedFlag();
101 geolocYAccessor.ResetModifiedFlag();
102 backMapXAccessor.ResetModifiedFlag();
103 backMapYAccessor.ResetModifiedFlag();
104
105 FreeWghtsBackMap();
106
107 delete m_poGeolocTmpDataset;
108 delete m_poBackmapTmpDataset;
109}
110
111/************************************************************************/
112/* AllocateBackMap() */
113/************************************************************************/
114
115 bool GDALGeoLocDatasetAccessors::AllocateBackMap()
116{
117 auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
118 if (poDriver == nullptr)
119 return false;
120
121 // CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
122 // so store them in a long-lived std::string
123 const std::string osBackmapTmpFilename =
124 CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif");
125 m_poBackmapTmpDataset = poDriver->Create(
126 osBackmapTmpFilename.c_str(), m_psTransform->nBackMapWidth,
127 m_psTransform->nBackMapHeight, 2, GDT_Float32,
128 m_aosGTiffCreationOptions.List());
129 if (m_poBackmapTmpDataset == nullptr)
130 {
131 return false;
132 }
133 m_poBackmapTmpDataset->MarkSuppressOnClose();
134 VSIUnlink(m_poBackmapTmpDataset->GetDescription());
135 auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
136 auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
137
138 backMapXAccessor.SetBand(poBandX);
139 backMapYAccessor.SetBand(poBandY);
140
141 // CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
142 // so store them in a long-lived std::string
143 const std::string osBackmapWeightsTmpFilename =
144 CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif");
145 m_poBackmapWeightsTmpDataset = poDriver->Create(
146 osBackmapWeightsTmpFilename.c_str(), m_psTransform->nBackMapWidth,
147 m_psTransform->nBackMapHeight, 1, GDT_Float32,
148 m_aosGTiffCreationOptions.List());
149 if (m_poBackmapWeightsTmpDataset == nullptr)
150 {
151 return false;
152 }
153 m_poBackmapWeightsTmpDataset->MarkSuppressOnClose();
154 VSIUnlink(m_poBackmapWeightsTmpDataset->GetDescription());
155 backMapWeightAccessor.SetBand(
156 m_poBackmapWeightsTmpDataset->GetRasterBand(1));
157
158 return true;
159}
160
161/************************************************************************/
162/* FreeWghtsBackMap() */
163/************************************************************************/
164
165 void GDALGeoLocDatasetAccessors::FreeWghtsBackMap()
166{
167 if (m_poBackmapWeightsTmpDataset)
168 {
169 backMapWeightAccessor.ResetModifiedFlag();
170 delete m_poBackmapWeightsTmpDataset;
171 m_poBackmapWeightsTmpDataset = nullptr;
172 }
173}
174
175/************************************************************************/
176/* GetBackmapDataset() */
177/************************************************************************/
178
179 GDALDataset *GDALGeoLocDatasetAccessors::GetBackmapDataset()
180{
181 auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
182 auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
183 poBandX->SetNoDataValue(INVALID_BMXY);
184 poBandY->SetNoDataValue(INVALID_BMXY);
185 return m_poBackmapTmpDataset;
186}
187
188/************************************************************************/
189/* FlushBackmapCaches() */
190/************************************************************************/
191
192 void GDALGeoLocDatasetAccessors::FlushBackmapCaches()
193{
194 backMapXAccessor.FlushCache();
195 backMapYAccessor.FlushCache();
196}
197
198/************************************************************************/
199/* Load() */
200/************************************************************************/
201
202 bool GDALGeoLocDatasetAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
203{
204 return LoadGeoloc(bIsRegularGrid) &&
205 ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
206 (!bUseQuadtree &&
207 GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
208}
209
210/************************************************************************/
211/* LoadGeoloc() */
212/************************************************************************/
213
214 bool GDALGeoLocDatasetAccessors::LoadGeoloc(bool bIsRegularGrid)
215
216{
217 if (bIsRegularGrid)
218 {
219 const int nXSize = m_psTransform->nGeoLocXSize;
220 const int nYSize = m_psTransform->nGeoLocYSize;
221
222 auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
223 if (poDriver == nullptr)
224 return false;
225
226 // CPLResetExtension / CPLGenerateTempFilename generate short-lived
227 // strings, so store them in a long-lived std::string
228 const std::string osGeolocTmpFilename =
229 CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif");
230 m_poGeolocTmpDataset =
231 poDriver->Create(osGeolocTmpFilename.c_str(), nXSize, nYSize, 2,
232 GDT_Float64, m_aosGTiffCreationOptions.List());
233 if (m_poGeolocTmpDataset == nullptr)
234 {
235 return false;
236 }
237 m_poGeolocTmpDataset->MarkSuppressOnClose();
238 VSIUnlink(m_poGeolocTmpDataset->GetDescription());
239
240 auto poXBand = m_poGeolocTmpDataset->GetRasterBand(1);
241 auto poYBand = m_poGeolocTmpDataset->GetRasterBand(2);
242
243 // Case of regular grid.
244 // The XBAND contains the x coordinates for all lines.
245 // The YBAND contains the y coordinates for all columns.
246
247 double *padfTempX =
248 static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
249 double *padfTempY =
250 static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
251 if (padfTempX == nullptr || padfTempY == nullptr)
252 {
253 CPLFree(padfTempX);
254 CPLFree(padfTempY);
255 return false;
256 }
257
258 CPLErr eErr =
259 GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
260 padfTempX, nXSize, 1, GDT_Float64, 0, 0);
261
262 for (int j = 0; j < nYSize; j++)
263 {
264 if (poXBand->RasterIO(GF_Write, 0, j, nXSize, 1, padfTempX, nXSize,
265 1, GDT_Float64, 0, 0, nullptr) != CE_None)
266 {
267 eErr = CE_Failure;
268 break;
269 }
270 }
271
272 if (eErr == CE_None)
273 {
274 eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
275 1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
276
277 for (int i = 0; i < nXSize; i++)
278 {
279 if (poYBand->RasterIO(GF_Write, i, 0, 1, nYSize, padfTempY, 1,
280 nYSize, GDT_Float64, 0, 0,
281 nullptr) != CE_None)
282 {
283 eErr = CE_Failure;
284 break;
285 }
286 }
287 }
288
289 CPLFree(padfTempX);
290 CPLFree(padfTempY);
291
292 if (eErr != CE_None)
293 return false;
294
295 geolocXAccessor.SetBand(poXBand);
296 geolocYAccessor.SetBand(poYBand);
297 }
298 else
299 {
300 geolocXAccessor.SetBand(
301 GDALRasterBand::FromHandle(m_psTransform->hBand_X));
302 geolocYAccessor.SetBand(
303 GDALRasterBand::FromHandle(m_psTransform->hBand_Y));
304 }
305
306 GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(m_psTransform);
307 return true;
308}
309
CPLStringList
String list class designed around our use of C "char**" string lists.
Definition cpl_string.h:449
GDALCachedPixelAccessor
Class to have reasonably fast random pixel access to a raster band, when accessing multiple pixels th...
Definition gdalcachedpixelaccessor.h:53
GDALCachedPixelAccessor::FlushCache
bool FlushCache()
Flush content of modified tiles and drop caches.
Definition gdalcachedpixelaccessor.h:281
GDALCachedPixelAccessor::SetBand
void SetBand(GDALRasterBand *poBand)
Assign the raster band if not known at construction time.
Definition gdalcachedpixelaccessor.h:84
GDALCachedPixelAccessor::ResetModifiedFlag
void ResetModifiedFlag()
Reset the modified flag for cached tiles.
Definition gdalcachedpixelaccessor.h:302
GDALDataset
A set of associated raster bands, usually from one file.
Definition gdal_priv.h:490
GDALDataset::MarkSuppressOnClose
void MarkSuppressOnClose()
Set that the dataset must be deleted on close.
Definition gdaldataset.cpp:1723
GDALDataset::GetRasterBand
GDALRasterBand * GetRasterBand(int)
Fetch a band object for a dataset.
Definition gdaldataset.cpp:1013
GDALDriver::FromHandle
static GDALDriver * FromHandle(GDALDriverH hDriver)
Convert a GDALDriverH to a GDALDriver*.
Definition gdal_priv.h:2173
GDALMajorObject::GetDescription
virtual const char * GetDescription() const
Fetch object description.
Definition gdalmajorobject.cpp:77
GDALRasterBand::SetNoDataValue
virtual CPLErr SetNoDataValue(double dfNoData)
Set the no data value for this band.
Definition gdalrasterband.cpp:1995
GDALRasterBand::FromHandle
static GDALRasterBand * FromHandle(GDALRasterBandH hBand)
Convert a GDALRasterBandH to a GDALRasterBand*.
Definition gdal_priv.h:1764
CPLFree
#define CPLFree
Alias of VSIFree()
Definition cpl_conv.h:98
CPLGenerateTempFilename
const char * CPLGenerateTempFilename(const char *pszStem)
Generate temporary file name.
Definition cpl_path.cpp:1110
CPLResetExtension
const char * CPLResetExtension(const char *, const char *)
Replace the extension with the provided one.
Definition cpl_path.cpp:466
CPLErr
CPLErr
Error category.
Definition cpl_error.h:53
VSIUnlink
int VSIUnlink(const char *pszFilename)
Delete a file.
Definition cpl_vsil.cpp:450
VSI_MALLOC2_VERBOSE
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition cpl_vsi.h:353
GDT_Float64
@ GDT_Float64
Definition gdal.h:75
GDT_Float32
@ GDT_Float32
Definition gdal.h:74
GF_Write
@ GF_Write
Definition gdal.h:134
GF_Read
@ GF_Read
Definition gdal.h:133
GDALGetDriverByName
GDALDriverH GDALGetDriverByName(const char *)
Fetch a driver based on the short name.
Definition gdaldrivermanager.cpp:691
GDALRasterIO
CPLErr GDALRasterIO(GDALRasterBandH hRBand, GDALRWFlag eRWFlag, int nDSXOff, int nDSYOff, int nDSXSize, int nDSYSize, void *pBuffer, int nBXSize, int nBYSize, GDALDataType eBDataType, int nPixelSpace, int nLineSpace)
Read/write a region of image data for this band.
Definition gdalrasterband.cpp:448