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 : * 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 :
31 : /*! @cond Doxygen_Suppress */
32 :
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 :
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 2 : explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo *psTransform)
65 2 : : m_psTransform(psTransform), geolocXAccessor(nullptr),
66 : geolocYAccessor(nullptr), backMapXAccessor(nullptr),
67 2 : backMapYAccessor(nullptr), backMapWeightAccessor(nullptr)
68 : {
69 2 : m_aosGTiffCreationOptions.SetNameValue("TILED", "YES");
70 2 : m_aosGTiffCreationOptions.SetNameValue("INTERLEAVE", "BAND");
71 : m_aosGTiffCreationOptions.SetNameValue("BLOCKXSIZE",
72 2 : CPLSPrintf("%d", TILE_SIZE));
73 : m_aosGTiffCreationOptions.SetNameValue("BLOCKYSIZE",
74 2 : CPLSPrintf("%d", TILE_SIZE));
75 2 : }
76 :
77 : ~GDALGeoLocDatasetAccessors();
78 :
79 : bool Load(bool bIsRegularGrid, bool bUseQuadtree);
80 :
81 : bool AllocateBackMap();
82 :
83 : GDALDataset *GetBackmapDataset();
84 : void FlushBackmapCaches();
85 :
86 2 : static void ReleaseBackmapDataset(GDALDataset *)
87 : {
88 2 : }
89 :
90 : void FreeWghtsBackMap();
91 : };
92 :
93 : /************************************************************************/
94 : /* ~GDALGeoLocDatasetAccessors() */
95 : /************************************************************************/
96 :
97 2 : GDALGeoLocDatasetAccessors::~GDALGeoLocDatasetAccessors()
98 : {
99 2 : geolocXAccessor.ResetModifiedFlag();
100 2 : geolocYAccessor.ResetModifiedFlag();
101 2 : backMapXAccessor.ResetModifiedFlag();
102 2 : backMapYAccessor.ResetModifiedFlag();
103 :
104 2 : FreeWghtsBackMap();
105 :
106 2 : delete m_poGeolocTmpDataset;
107 2 : delete m_poBackmapTmpDataset;
108 2 : }
109 :
110 : /************************************************************************/
111 : /* AllocateBackMap() */
112 : /************************************************************************/
113 :
114 2 : bool GDALGeoLocDatasetAccessors::AllocateBackMap()
115 : {
116 2 : auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
117 2 : if (poDriver == nullptr)
118 0 : return false;
119 :
120 2 : m_poBackmapTmpDataset = poDriver->Create(
121 : CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif"),
122 2 : m_psTransform->nBackMapWidth, m_psTransform->nBackMapHeight, 2,
123 2 : GDT_Float32, m_aosGTiffCreationOptions.List());
124 2 : if (m_poBackmapTmpDataset == nullptr)
125 : {
126 0 : return false;
127 : }
128 2 : m_poBackmapTmpDataset->MarkSuppressOnClose();
129 2 : VSIUnlink(m_poBackmapTmpDataset->GetDescription());
130 2 : auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
131 2 : auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
132 :
133 2 : backMapXAccessor.SetBand(poBandX);
134 2 : backMapYAccessor.SetBand(poBandY);
135 :
136 2 : m_poBackmapWeightsTmpDataset = poDriver->Create(
137 : CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif"),
138 2 : m_psTransform->nBackMapWidth, m_psTransform->nBackMapHeight, 1,
139 2 : GDT_Float32, m_aosGTiffCreationOptions.List());
140 2 : if (m_poBackmapWeightsTmpDataset == nullptr)
141 : {
142 0 : return false;
143 : }
144 2 : m_poBackmapWeightsTmpDataset->MarkSuppressOnClose();
145 2 : VSIUnlink(m_poBackmapWeightsTmpDataset->GetDescription());
146 2 : backMapWeightAccessor.SetBand(
147 2 : m_poBackmapWeightsTmpDataset->GetRasterBand(1));
148 :
149 2 : return true;
150 : }
151 :
152 : /************************************************************************/
153 : /* FreeWghtsBackMap() */
154 : /************************************************************************/
155 :
156 4 : void GDALGeoLocDatasetAccessors::FreeWghtsBackMap()
157 : {
158 4 : if (m_poBackmapWeightsTmpDataset)
159 : {
160 2 : backMapWeightAccessor.ResetModifiedFlag();
161 2 : delete m_poBackmapWeightsTmpDataset;
162 2 : m_poBackmapWeightsTmpDataset = nullptr;
163 : }
164 4 : }
165 :
166 : /************************************************************************/
167 : /* GetBackmapDataset() */
168 : /************************************************************************/
169 :
170 2 : GDALDataset *GDALGeoLocDatasetAccessors::GetBackmapDataset()
171 : {
172 2 : auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
173 2 : auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
174 2 : poBandX->SetNoDataValue(INVALID_BMXY);
175 2 : poBandY->SetNoDataValue(INVALID_BMXY);
176 2 : return m_poBackmapTmpDataset;
177 : }
178 :
179 : /************************************************************************/
180 : /* FlushBackmapCaches() */
181 : /************************************************************************/
182 :
183 2 : void GDALGeoLocDatasetAccessors::FlushBackmapCaches()
184 : {
185 2 : backMapXAccessor.FlushCache();
186 2 : backMapYAccessor.FlushCache();
187 2 : }
188 :
189 : /************************************************************************/
190 : /* Load() */
191 : /************************************************************************/
192 :
193 2 : bool GDALGeoLocDatasetAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
194 : {
195 4 : return LoadGeoloc(bIsRegularGrid) &&
196 0 : ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
197 2 : (!bUseQuadtree &&
198 4 : GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
199 : }
200 :
201 : /************************************************************************/
202 : /* LoadGeoloc() */
203 : /************************************************************************/
204 :
205 2 : bool GDALGeoLocDatasetAccessors::LoadGeoloc(bool bIsRegularGrid)
206 :
207 : {
208 2 : if (bIsRegularGrid)
209 : {
210 1 : const int nXSize = m_psTransform->nGeoLocXSize;
211 1 : const int nYSize = m_psTransform->nGeoLocYSize;
212 :
213 1 : auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
214 1 : if (poDriver == nullptr)
215 0 : return false;
216 :
217 1 : m_poGeolocTmpDataset = poDriver->Create(
218 : CPLResetExtension(CPLGenerateTempFilename(nullptr), "tif"), nXSize,
219 1 : nYSize, 2, GDT_Float64, m_aosGTiffCreationOptions.List());
220 1 : if (m_poGeolocTmpDataset == nullptr)
221 : {
222 0 : return false;
223 : }
224 1 : m_poGeolocTmpDataset->MarkSuppressOnClose();
225 1 : VSIUnlink(m_poGeolocTmpDataset->GetDescription());
226 :
227 1 : auto poXBand = m_poGeolocTmpDataset->GetRasterBand(1);
228 1 : auto poYBand = m_poGeolocTmpDataset->GetRasterBand(2);
229 :
230 : // Case of regular grid.
231 : // The XBAND contains the x coordinates for all lines.
232 : // The YBAND contains the y coordinates for all columns.
233 :
234 : double *padfTempX =
235 1 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
236 : double *padfTempY =
237 1 : static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
238 1 : if (padfTempX == nullptr || padfTempY == nullptr)
239 : {
240 0 : CPLFree(padfTempX);
241 0 : CPLFree(padfTempY);
242 0 : return false;
243 : }
244 :
245 : CPLErr eErr =
246 1 : GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
247 : padfTempX, nXSize, 1, GDT_Float64, 0, 0);
248 :
249 81 : for (int j = 0; j < nYSize; j++)
250 : {
251 80 : if (poXBand->RasterIO(GF_Write, 0, j, nXSize, 1, padfTempX, nXSize,
252 80 : 1, GDT_Float64, 0, 0, nullptr) != CE_None)
253 : {
254 0 : eErr = CE_Failure;
255 0 : break;
256 : }
257 : }
258 :
259 1 : if (eErr == CE_None)
260 : {
261 1 : eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
262 : 1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
263 :
264 361 : for (int i = 0; i < nXSize; i++)
265 : {
266 360 : if (poYBand->RasterIO(GF_Write, i, 0, 1, nYSize, padfTempY, 1,
267 : nYSize, GDT_Float64, 0, 0,
268 360 : nullptr) != CE_None)
269 : {
270 0 : eErr = CE_Failure;
271 0 : break;
272 : }
273 : }
274 : }
275 :
276 1 : CPLFree(padfTempX);
277 1 : CPLFree(padfTempY);
278 :
279 1 : if (eErr != CE_None)
280 0 : return false;
281 :
282 1 : geolocXAccessor.SetBand(poXBand);
283 1 : geolocYAccessor.SetBand(poYBand);
284 : }
285 : else
286 : {
287 1 : geolocXAccessor.SetBand(
288 1 : GDALRasterBand::FromHandle(m_psTransform->hBand_X));
289 1 : geolocYAccessor.SetBand(
290 1 : GDALRasterBand::FromHandle(m_psTransform->hBand_Y));
291 : }
292 :
293 2 : GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(m_psTransform);
294 2 : return true;
295 : }
296 :
297 : /*! @endcond */
|
