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