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& operator= (const GDALGeoLocDatasetAccessors&) = delete;
51
52 bool LoadGeoloc(bool bIsRegularGrid);
53
54 public:
55
56 static constexpr int TILE_SIZE = 1024;
57
58 GDALCachedPixelAccessor<double, TILE_SIZE> geolocXAccessor;
59 GDALCachedPixelAccessor<double, TILE_SIZE> geolocYAccessor;
60 GDALCachedPixelAccessor<float, TILE_SIZE> backMapXAccessor;
61 GDALCachedPixelAccessor<float, TILE_SIZE> backMapYAccessor;
62 GDALCachedPixelAccessor<float, TILE_SIZE> backMapWeightAccessor;
63
64 explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo* psTransform):
65 m_psTransform(psTransform),
66 geolocXAccessor(nullptr),
67 geolocYAccessor(nullptr),
68 backMapXAccessor(nullptr),
69 backMapYAccessor(nullptr),
70 backMapWeightAccessor(nullptr)
71 {
72 m_aosGTiffCreationOptions.SetNameValue("TILED", "YES");
73 m_aosGTiffCreationOptions.SetNameValue("INTERLEAVE", "BAND");
74 m_aosGTiffCreationOptions.SetNameValue("BLOCKXSIZE", CPLSPrintf("%d", TILE_SIZE));
75 m_aosGTiffCreationOptions.SetNameValue("BLOCKYSIZE", 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 static void ReleaseBackmapDataset(GDALDataset*) {}
87
88 void FreeWghtsBackMap();
89 };
90
91 /************************************************************************/
92 /* ~GDALGeoLocDatasetAccessors() */
93 /************************************************************************/
94
95 GDALGeoLocDatasetAccessors::~GDALGeoLocDatasetAccessors()
96 {
97 geolocXAccessor.ResetModifiedFlag();
98 geolocYAccessor.ResetModifiedFlag();
99 backMapXAccessor.ResetModifiedFlag();
100 backMapYAccessor.ResetModifiedFlag();
101
102 FreeWghtsBackMap();
103
104 delete m_poGeolocTmpDataset;
105 delete m_poBackmapTmpDataset;
106 }
107
108 /************************************************************************/
109 /* AllocateBackMap() */
110 /************************************************************************/
111
112 bool GDALGeoLocDatasetAccessors::AllocateBackMap()
113 {
114 auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
115 if( poDriver == nullptr )
116 return false;
117
118 m_poBackmapTmpDataset = poDriver->Create(
119 CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif"),
120 m_psTransform->nBackMapWidth,
121 m_psTransform->nBackMapHeight,
122 2,
123 GDT_Float32,
124 m_aosGTiffCreationOptions.List());
125 if( m_poBackmapTmpDataset == nullptr )
126 {
127 return false;
128 }
129 m_poBackmapTmpDataset->MarkSuppressOnClose();
130 VSIUnlink(m_poBackmapTmpDataset->GetDescription());
131 auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
132 auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
133
134 backMapXAccessor.SetBand(poBandX);
135 backMapYAccessor.SetBand(poBandY);
136
137 m_poBackmapWeightsTmpDataset = poDriver->Create(
138 CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif"),
139 m_psTransform->nBackMapWidth,
140 m_psTransform->nBackMapHeight,
141 1,
142 GDT_Float32,
143 m_aosGTiffCreationOptions.List());
144 if( m_poBackmapWeightsTmpDataset == nullptr )
145 {
146 return false;
147 }
148 m_poBackmapWeightsTmpDataset->MarkSuppressOnClose();
149 VSIUnlink(m_poBackmapWeightsTmpDataset->GetDescription());
150 backMapWeightAccessor.SetBand(m_poBackmapWeightsTmpDataset->GetRasterBand(1));
151
152 return true;
153 }
154
155 /************************************************************************/
156 /* FreeWghtsBackMap() */
157 /************************************************************************/
158
159 void GDALGeoLocDatasetAccessors::FreeWghtsBackMap()
160 {
161 if( m_poBackmapWeightsTmpDataset )
162 {
163 backMapWeightAccessor.ResetModifiedFlag();
164 delete m_poBackmapWeightsTmpDataset;
165 m_poBackmapWeightsTmpDataset = nullptr;
166 }
167 }
168
169 /************************************************************************/
170 /* GetBackmapDataset() */
171 /************************************************************************/
172
173 GDALDataset* GDALGeoLocDatasetAccessors::GetBackmapDataset()
174 {
175 auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
176 auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
177 poBandX->SetNoDataValue(INVALID_BMXY);
178 poBandY->SetNoDataValue(INVALID_BMXY);
179 return m_poBackmapTmpDataset;
180 }
181
182 /************************************************************************/
183 /* FlushBackmapCaches() */
184 /************************************************************************/
185
186 void GDALGeoLocDatasetAccessors::FlushBackmapCaches()
187 {
188 backMapXAccessor.FlushCache();
189 backMapYAccessor.FlushCache();
190 }
191
192 /************************************************************************/
193 /* Load() */
194 /************************************************************************/
195
196 bool GDALGeoLocDatasetAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
197 {
198 return LoadGeoloc(bIsRegularGrid) &&
199 ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
200 (!bUseQuadtree && GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
201 }
202
203 /************************************************************************/
204 /* LoadGeoloc() */
205 /************************************************************************/
206
207 bool GDALGeoLocDatasetAccessors::LoadGeoloc(bool bIsRegularGrid)
208
209 {
210 if( bIsRegularGrid )
211 {
212 const int nXSize = m_psTransform->nGeoLocXSize;
213 const int nYSize = m_psTransform->nGeoLocYSize;
214
215 auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
216 if( poDriver == nullptr )
217 return false;
218
219 m_poGeolocTmpDataset = poDriver->Create(
220 CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif"),
221 nXSize,
222 nYSize,
223 2,
224 GDT_Float64,
225 m_aosGTiffCreationOptions.List());
226 if( m_poGeolocTmpDataset == nullptr )
227 {
228 return false;
229 }
230 m_poGeolocTmpDataset->MarkSuppressOnClose();
231 VSIUnlink(m_poGeolocTmpDataset->GetDescription());
232
233 auto poXBand = m_poGeolocTmpDataset->GetRasterBand(1);
234 auto poYBand = m_poGeolocTmpDataset->GetRasterBand(2);
235
236 // Case of regular grid.
237 // The XBAND contains the x coordinates for all lines.
238 // The YBAND contains the y coordinates for all columns.
239
240 double* padfTempX = static_cast<double *>(
241 VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
242 double* padfTempY = static_cast<double *>(
243 VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
244 if( padfTempX == nullptr || padfTempY == nullptr )
245 {
246 CPLFree(padfTempX);
247 CPLFree(padfTempY);
248 return false;
249 }
250
251 CPLErr eErr =
252 GDALRasterIO( m_psTransform->hBand_X, GF_Read,
253 0, 0, nXSize, 1,
254 padfTempX, nXSize, 1,
255 GDT_Float64, 0, 0 );
256
257 for( int j = 0; j < nYSize; j++ )
258 {
259 if( poXBand->RasterIO(GF_Write,
260 0, j, nXSize, 1,
261 padfTempX,
262 nXSize, 1,
263 GDT_Float64,
264 0, 0, nullptr) != CE_None )
265 {
266 eErr = CE_Failure;
267 break;
268 }
269 }
270
271 if( eErr == CE_None )
272 {
273 eErr = GDALRasterIO( m_psTransform->hBand_Y, GF_Read,
274 0, 0, nYSize, 1,
275 padfTempY, nYSize, 1,
276 GDT_Float64, 0, 0 );
277
278 for( int i = 0; i < nXSize; i++ )
279 {
280 if( poYBand->RasterIO(GF_Write,
281 i, 0, 1, nYSize,
282 padfTempY,
283 1, nYSize,
284 GDT_Float64,
285 0, 0, nullptr) != CE_None )
286 {
287 eErr = CE_Failure;
288 break;
289 }
290 }
291 }
292
293 CPLFree(padfTempX);
294 CPLFree(padfTempY);
295
296 if( eErr != CE_None )
297 return false;
298
299 geolocXAccessor.SetBand( poXBand );
300 geolocYAccessor.SetBand( poYBand );
301
302 }
303 else
304 {
305 geolocXAccessor.SetBand( GDALRasterBand::FromHandle(m_psTransform->hBand_X) );
306 geolocYAccessor.SetBand( GDALRasterBand::FromHandle(m_psTransform->hBand_Y) );
307 }
308
309 return GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(m_psTransform);
310 }
311
GDALDataset::GetRasterBand
GDALRasterBand * GetRasterBand(int)
Fetch a band object for a dataset.
Definition: gdaldataset.cpp:772
GDALRasterBand::SetNoDataValue
virtual CPLErr SetNoDataValue(double dfNoData)
Set the no data value for this band.
Definition: gdalrasterband.cpp:1850
CPLGenerateTempFilename
const char * CPLGenerateTempFilename(const char *pszStem)
Generate temporary file name.
Definition: cpl_path.cpp:1119
GDALCachedPixelAccessor
Class to have reasonably fast random pixel access to a raster band, when accessing multiple pixels th...
Definition: gdalcachedpixelaccessor.h:51
GDT_Float32
@ GDT_Float32
Definition: gdal.h:72
GDALCachedPixelAccessor::ResetModifiedFlag
void ResetModifiedFlag()
Reset the modified flag for cached tiles.
Definition: gdalcachedpixelaccessor.h:296
CPLStringList
String list class designed around our use of C "char**" string lists.
Definition: cpl_string.h:429
CPLResetExtension
const char * CPLResetExtension(const char *, const char *)
Replace the extension with the provided one.
Definition: cpl_path.cpp:442
GF_Read
@ GF_Read
Definition: gdal.h:123
GDALDataset
A set of associated raster bands, usually from one file.
Definition: gdal_priv.h:342
VSI_MALLOC2_VERBOSE
#define VSI_MALLOC2_VERBOSE(nSize1, nSize2)
VSI_MALLOC2_VERBOSE.
Definition: cpl_vsi.h:318
GF_Write
@ GF_Write
Definition: gdal.h:124
GDT_Float64
@ GDT_Float64
Definition: gdal.h:73
VSIUnlink
int VSIUnlink(const char *pszFilename)
Delete a file.
Definition: cpl_vsil.cpp:422
GDALGetDriverByName
GDALDriverH GDALGetDriverByName(const char *)
Fetch a driver based on the short name.
Definition: gdaldrivermanager.cpp:620
CPLSPrintf
const char * CPLSPrintf(const char *fmt,...)
CPLSPrintf() that works with 10 static buffer.
Definition: cpl_string.cpp:993
GDALCachedPixelAccessor::FlushCache
bool FlushCache()
Flush content of modified tiles and drop caches.
Definition: gdalcachedpixelaccessor.h:276
GDALDataset::MarkSuppressOnClose
void MarkSuppressOnClose()
Set that the dataset must be deleted on close.
Definition: gdaldataset.cpp:1431
GDALDriver::FromHandle
static GDALDriver * FromHandle(GDALDriverH hDriver)
Convert a GDALDriverH to a GDALDriver*.
Definition: gdal_priv.h:1695
CPLErr
CPLErr
Error category.
Definition: cpl_error.h:52
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:400
GDALCachedPixelAccessor::SetBand
void SetBand(GDALRasterBand *poBand)
Assign the raster band if not known at construction time.
Definition: gdalcachedpixelaccessor.h:110
CPLFree
#define CPLFree
Alias of VSIFree()
Definition: cpl_conv.h:83
GDALMajorObject::GetDescription
virtual const char * GetDescription() const
Fetch object description.
Definition: gdalmajorobject.cpp:78
GDALRasterBand::FromHandle
static GDALRasterBand * FromHandle(GDALRasterBandH hBand)
Convert a GDALRasterBandH to a GDALRasterBand*.
Definition: gdal_priv.h:1410