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