/src/gdal/alg/gdalgeoloc_dataset_accessor.h

Source (jump to first uncovered line)
/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  Dataset storage of geolocation array and backmap
 * Author:   Even Rouault, <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2022, Planet Labs
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "gdalcachedpixelaccessor.h"

/*! @cond Doxygen_Suppress */

/************************************************************************/
/*                        GDALGeoLocDatasetAccessors                    */
/************************************************************************/

class GDALGeoLocDatasetAccessors
{
    typedef class GDALGeoLocDatasetAccessors AccessorType;

    GDALGeoLocTransformInfo *m_psTransform;

    CPLStringList m_aosGTiffCreationOptions{};

    GDALDataset *m_poGeolocTmpDataset = nullptr;
    GDALDataset *m_poBackmapTmpDataset = nullptr;
    GDALDataset *m_poBackmapWeightsTmpDataset = nullptr;

    GDALGeoLocDatasetAccessors(const GDALGeoLocDatasetAccessors &) = delete;
    GDALGeoLocDatasetAccessors &
    operator=(const GDALGeoLocDatasetAccessors &) = delete;

    bool LoadGeoloc(bool bIsRegularGrid);

  public:
    static constexpr int TILE_SIZE = 256;
    static constexpr int TILE_COUNT = 64;

    GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocXAccessor;
    GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocYAccessor;
    GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapXAccessor;
    GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapYAccessor;
    GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapWeightAccessor;

    explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo *psTransform)
        : m_psTransform(psTransform), geolocXAccessor(nullptr),
          geolocYAccessor(nullptr), backMapXAccessor(nullptr),
          backMapYAccessor(nullptr), backMapWeightAccessor(nullptr)
    {
        m_aosGTiffCreationOptions.SetNameValue("TILED", "YES");
        m_aosGTiffCreationOptions.SetNameValue("INTERLEAVE", "BAND");
        m_aosGTiffCreationOptions.SetNameValue("BLOCKXSIZE",
                                               CPLSPrintf("%d", TILE_SIZE));
        m_aosGTiffCreationOptions.SetNameValue("BLOCKYSIZE",
                                               CPLSPrintf("%d", TILE_SIZE));
    }

    ~GDALGeoLocDatasetAccessors();

    bool Load(bool bIsRegularGrid, bool bUseQuadtree);

    bool AllocateBackMap();

    GDALDataset *GetBackmapDataset();
    void FlushBackmapCaches();

    static void ReleaseBackmapDataset(GDALDataset *)
    {
    }

    void FreeWghtsBackMap();
};

/************************************************************************/
/*                    ~GDALGeoLocDatasetAccessors()                     */
/************************************************************************/

GDALGeoLocDatasetAccessors::~GDALGeoLocDatasetAccessors()
{
    geolocXAccessor.ResetModifiedFlag();
    geolocYAccessor.ResetModifiedFlag();
    backMapXAccessor.ResetModifiedFlag();
    backMapYAccessor.ResetModifiedFlag();

    FreeWghtsBackMap();

    delete m_poGeolocTmpDataset;
    delete m_poBackmapTmpDataset;
}

/************************************************************************/
/*                         AllocateBackMap()                            */
/************************************************************************/

bool GDALGeoLocDatasetAccessors::AllocateBackMap()
{
    auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
    if (poDriver == nullptr)
        return false;

    // CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
    // so store them in a long-lived std::string
    const std::string osBackmapTmpFilename = CPLResetExtensionSafe(
        CPLGenerateTempFilenameSafe(nullptr).c_str(), "tif");
    m_poBackmapTmpDataset = poDriver->Create(
        osBackmapTmpFilename.c_str(), m_psTransform->nBackMapWidth,
        m_psTransform->nBackMapHeight, 2, GDT_Float32,
        m_aosGTiffCreationOptions.List());
    if (m_poBackmapTmpDataset == nullptr)
    {
        return false;
    }
    m_poBackmapTmpDataset->MarkSuppressOnClose();
    VSIUnlink(m_poBackmapTmpDataset->GetDescription());
    auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
    auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);

    backMapXAccessor.SetBand(poBandX);
    backMapYAccessor.SetBand(poBandY);

    // CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
    // so store them in a long-lived std::string
    const std::string osBackmapWeightsTmpFilename = CPLResetExtensionSafe(
        CPLGenerateTempFilenameSafe(nullptr).c_str(), "tif");
    m_poBackmapWeightsTmpDataset = poDriver->Create(
        osBackmapWeightsTmpFilename.c_str(), m_psTransform->nBackMapWidth,
        m_psTransform->nBackMapHeight, 1, GDT_Float32,
        m_aosGTiffCreationOptions.List());
    if (m_poBackmapWeightsTmpDataset == nullptr)
    {
        return false;
    }
    m_poBackmapWeightsTmpDataset->MarkSuppressOnClose();
    VSIUnlink(m_poBackmapWeightsTmpDataset->GetDescription());
    backMapWeightAccessor.SetBand(
        m_poBackmapWeightsTmpDataset->GetRasterBand(1));

    return true;
}

/************************************************************************/
/*                         FreeWghtsBackMap()                           */
/************************************************************************/

void GDALGeoLocDatasetAccessors::FreeWghtsBackMap()
{
    if (m_poBackmapWeightsTmpDataset)
    {
        backMapWeightAccessor.ResetModifiedFlag();
        delete m_poBackmapWeightsTmpDataset;
        m_poBackmapWeightsTmpDataset = nullptr;
    }
}

/************************************************************************/
/*                        GetBackmapDataset()                           */
/************************************************************************/

GDALDataset *GDALGeoLocDatasetAccessors::GetBackmapDataset()
{
    auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
    auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
    poBandX->SetNoDataValue(INVALID_BMXY);
    poBandY->SetNoDataValue(INVALID_BMXY);
    return m_poBackmapTmpDataset;
}

/************************************************************************/
/*                       FlushBackmapCaches()                           */
/************************************************************************/

void GDALGeoLocDatasetAccessors::FlushBackmapCaches()
{
    backMapXAccessor.FlushCache();
    backMapYAccessor.FlushCache();
}

/************************************************************************/
/*                             Load()                                   */
/************************************************************************/

bool GDALGeoLocDatasetAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
{
    return LoadGeoloc(bIsRegularGrid) &&
           ((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) ||
            (!bUseQuadtree &&
             GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
}

/************************************************************************/
/*                          LoadGeoloc()                                */
/************************************************************************/

bool GDALGeoLocDatasetAccessors::LoadGeoloc(bool bIsRegularGrid)

{
    if (bIsRegularGrid)
    {
        const int nXSize = m_psTransform->nGeoLocXSize;
        const int nYSize = m_psTransform->nGeoLocYSize;

        auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
        if (poDriver == nullptr)
            return false;

        // CPLResetExtension / CPLGenerateTempFilename generate short-lived
        // strings, so store them in a long-lived std::string
        const std::string osGeolocTmpFilename = CPLResetExtensionSafe(
            CPLGenerateTempFilenameSafe(nullptr).c_str(), "tif");
        m_poGeolocTmpDataset =
            poDriver->Create(osGeolocTmpFilename.c_str(), nXSize, nYSize, 2,
                             GDT_Float64, m_aosGTiffCreationOptions.List());
        if (m_poGeolocTmpDataset == nullptr)
        {
            return false;
        }
        m_poGeolocTmpDataset->MarkSuppressOnClose();
        VSIUnlink(m_poGeolocTmpDataset->GetDescription());

        auto poXBand = m_poGeolocTmpDataset->GetRasterBand(1);
        auto poYBand = m_poGeolocTmpDataset->GetRasterBand(2);

        // Case of regular grid.
        // The XBAND contains the x coordinates for all lines.
        // The YBAND contains the y coordinates for all columns.

        double *padfTempX =
            static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
        double *padfTempY =
            static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
        if (padfTempX == nullptr || padfTempY == nullptr)
        {
            CPLFree(padfTempX);
            CPLFree(padfTempY);
            return false;
        }

        CPLErr eErr =
            GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
                         padfTempX, nXSize, 1, GDT_Float64, 0, 0);

        for (int j = 0; j < nYSize; j++)
        {
            if (poXBand->RasterIO(GF_Write, 0, j, nXSize, 1, padfTempX, nXSize,
                                  1, GDT_Float64, 0, 0, nullptr) != CE_None)
            {
                eErr = CE_Failure;
                break;
            }
        }

        if (eErr == CE_None)
        {
            eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
                                1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);

            for (int i = 0; i < nXSize; i++)
            {
                if (poYBand->RasterIO(GF_Write, i, 0, 1, nYSize, padfTempY, 1,
                                      nYSize, GDT_Float64, 0, 0,
                                      nullptr) != CE_None)
                {
                    eErr = CE_Failure;
                    break;
                }
            }
        }

        CPLFree(padfTempX);
        CPLFree(padfTempY);

        if (eErr != CE_None)
            return false;

        geolocXAccessor.SetBand(poXBand);
        geolocYAccessor.SetBand(poYBand);
    }
    else
    {
        geolocXAccessor.SetBand(
            GDALRasterBand::FromHandle(m_psTransform->hBand_X));
        geolocYAccessor.SetBand(
            GDALRasterBand::FromHandle(m_psTransform->hBand_Y));
    }

    GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(m_psTransform);
    return true;
}

/*! @endcond */

Coverage Report

Created: 2025-06-13 06:29

Line	Count	Source (jump to first uncovered line)
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		* SPDX-License-Identifier: MIT
11		****************************************************************************/
12
13		#include "gdalcachedpixelaccessor.h"
14
15		/! @cond Doxygen_Suppress /
16
17		/************************************************************************/
18		/* GDALGeoLocDatasetAccessors */
19		/************************************************************************/
20
21		class GDALGeoLocDatasetAccessors
22		{
23		typedef class GDALGeoLocDatasetAccessors AccessorType;
24
25		GDALGeoLocTransformInfo *m_psTransform;
26
27		CPLStringList m_aosGTiffCreationOptions{};
28
29		GDALDataset *m_poGeolocTmpDataset = nullptr;
30		GDALDataset *m_poBackmapTmpDataset = nullptr;
31		GDALDataset *m_poBackmapWeightsTmpDataset = nullptr;
32
33		GDALGeoLocDatasetAccessors(const GDALGeoLocDatasetAccessors &) = delete;
34		GDALGeoLocDatasetAccessors &
35		operator=(const GDALGeoLocDatasetAccessors &) = delete;
36
37		bool LoadGeoloc(bool bIsRegularGrid);
38
39		public:
40		static constexpr int TILE_SIZE = 256;
41		static constexpr int TILE_COUNT = 64;
42
43		GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocXAccessor;
44		GDALCachedPixelAccessor<double, TILE_SIZE, TILE_COUNT> geolocYAccessor;
45		GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapXAccessor;
46		GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapYAccessor;
47		GDALCachedPixelAccessor<float, TILE_SIZE, TILE_COUNT> backMapWeightAccessor;
48
49		explicit GDALGeoLocDatasetAccessors(GDALGeoLocTransformInfo *psTransform)
50	0	: m_psTransform(psTransform), geolocXAccessor(nullptr),
51	0	geolocYAccessor(nullptr), backMapXAccessor(nullptr),
52	0	backMapYAccessor(nullptr), backMapWeightAccessor(nullptr)
53	0	{
54	0	m_aosGTiffCreationOptions.SetNameValue("TILED", "YES");
55	0	m_aosGTiffCreationOptions.SetNameValue("INTERLEAVE", "BAND");
56	0	m_aosGTiffCreationOptions.SetNameValue("BLOCKXSIZE",
57	0	CPLSPrintf("%d", TILE_SIZE));
58	0	m_aosGTiffCreationOptions.SetNameValue("BLOCKYSIZE",
59	0	CPLSPrintf("%d", TILE_SIZE));
60	0	}
61
62		~GDALGeoLocDatasetAccessors();
63
64		bool Load(bool bIsRegularGrid, bool bUseQuadtree);
65
66		bool AllocateBackMap();
67
68		GDALDataset *GetBackmapDataset();
69		void FlushBackmapCaches();
70
71		static void ReleaseBackmapDataset(GDALDataset *)
72	0	{
73	0	}
74
75		void FreeWghtsBackMap();
76		};
77
78		/************************************************************************/
79		/* ~GDALGeoLocDatasetAccessors() */
80		/************************************************************************/
81
82		GDALGeoLocDatasetAccessors::~GDALGeoLocDatasetAccessors()
83	0	{
84	0	geolocXAccessor.ResetModifiedFlag();
85	0	geolocYAccessor.ResetModifiedFlag();
86	0	backMapXAccessor.ResetModifiedFlag();
87	0	backMapYAccessor.ResetModifiedFlag();
88
89	0	FreeWghtsBackMap();
90
91	0	delete m_poGeolocTmpDataset;
92	0	delete m_poBackmapTmpDataset;
93	0	}
94
95		/************************************************************************/
96		/* AllocateBackMap() */
97		/************************************************************************/
98
99		bool GDALGeoLocDatasetAccessors::AllocateBackMap()
100	0	{
101	0	auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
102	0	if (poDriver == nullptr)
103	0	return false;
104
105		// CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
106		// so store them in a long-lived std::string
107	0	const std::string osBackmapTmpFilename = CPLResetExtensionSafe(
108	0	CPLGenerateTempFilenameSafe(nullptr).c_str(), "tif");
109	0	m_poBackmapTmpDataset = poDriver->Create(
110	0	osBackmapTmpFilename.c_str(), m_psTransform->nBackMapWidth,
111	0	m_psTransform->nBackMapHeight, 2, GDT_Float32,
112	0	m_aosGTiffCreationOptions.List());
113	0	if (m_poBackmapTmpDataset == nullptr)
114	0	{
115	0	return false;
116	0	}
117	0	m_poBackmapTmpDataset->MarkSuppressOnClose();
118	0	VSIUnlink(m_poBackmapTmpDataset->GetDescription());
119	0	auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
120	0	auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
121
122	0	backMapXAccessor.SetBand(poBandX);
123	0	backMapYAccessor.SetBand(poBandY);
124
125		// CPLResetExtension / CPLGenerateTempFilename generate short-lived strings,
126		// so store them in a long-lived std::string
127	0	const std::string osBackmapWeightsTmpFilename = CPLResetExtensionSafe(
128	0	CPLGenerateTempFilenameSafe(nullptr).c_str(), "tif");
129	0	m_poBackmapWeightsTmpDataset = poDriver->Create(
130	0	osBackmapWeightsTmpFilename.c_str(), m_psTransform->nBackMapWidth,
131	0	m_psTransform->nBackMapHeight, 1, GDT_Float32,
132	0	m_aosGTiffCreationOptions.List());
133	0	if (m_poBackmapWeightsTmpDataset == nullptr)
134	0	{
135	0	return false;
136	0	}
137	0	m_poBackmapWeightsTmpDataset->MarkSuppressOnClose();
138	0	VSIUnlink(m_poBackmapWeightsTmpDataset->GetDescription());
139	0	backMapWeightAccessor.SetBand(
140	0	m_poBackmapWeightsTmpDataset->GetRasterBand(1));
141
142	0	return true;
143	0	}
144
145		/************************************************************************/
146		/* FreeWghtsBackMap() */
147		/************************************************************************/
148
149		void GDALGeoLocDatasetAccessors::FreeWghtsBackMap()
150	0	{
151	0	if (m_poBackmapWeightsTmpDataset)
152	0	{
153	0	backMapWeightAccessor.ResetModifiedFlag();
154	0	delete m_poBackmapWeightsTmpDataset;
155	0	m_poBackmapWeightsTmpDataset = nullptr;
156	0	}
157	0	}
158
159		/************************************************************************/
160		/* GetBackmapDataset() */
161		/************************************************************************/
162
163		GDALDataset *GDALGeoLocDatasetAccessors::GetBackmapDataset()
164	0	{
165	0	auto poBandX = m_poBackmapTmpDataset->GetRasterBand(1);
166	0	auto poBandY = m_poBackmapTmpDataset->GetRasterBand(2);
167	0	poBandX->SetNoDataValue(INVALID_BMXY);
168	0	poBandY->SetNoDataValue(INVALID_BMXY);
169	0	return m_poBackmapTmpDataset;
170	0	}
171
172		/************************************************************************/
173		/* FlushBackmapCaches() */
174		/************************************************************************/
175
176		void GDALGeoLocDatasetAccessors::FlushBackmapCaches()
177	0	{
178	0	backMapXAccessor.FlushCache();
179	0	backMapYAccessor.FlushCache();
180	0	}
181
182		/************************************************************************/
183		/* Load() */
184		/************************************************************************/
185
186		bool GDALGeoLocDatasetAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
187	0	{
188	0	return LoadGeoloc(bIsRegularGrid) &&
189	0	((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) \|\|
190	0	(!bUseQuadtree &&
191	0	GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
192	0	}
193
194		/************************************************************************/
195		/* LoadGeoloc() */
196		/************************************************************************/
197
198		bool GDALGeoLocDatasetAccessors::LoadGeoloc(bool bIsRegularGrid)
199
200	0	{
201	0	if (bIsRegularGrid)
202	0	{
203	0	const int nXSize = m_psTransform->nGeoLocXSize;
204	0	const int nYSize = m_psTransform->nGeoLocYSize;
205
206	0	auto poDriver = GDALDriver::FromHandle(GDALGetDriverByName("GTiff"));
207	0	if (poDriver == nullptr)
208	0	return false;
209
210		// CPLResetExtension / CPLGenerateTempFilename generate short-lived
211		// strings, so store them in a long-lived std::string
212	0	const std::string osGeolocTmpFilename = CPLResetExtensionSafe(
213	0	CPLGenerateTempFilenameSafe(nullptr).c_str(), "tif");
214	0	m_poGeolocTmpDataset =
215	0	poDriver->Create(osGeolocTmpFilename.c_str(), nXSize, nYSize, 2,
216	0	GDT_Float64, m_aosGTiffCreationOptions.List());
217	0	if (m_poGeolocTmpDataset == nullptr)
218	0	{
219	0	return false;
220	0	}
221	0	m_poGeolocTmpDataset->MarkSuppressOnClose();
222	0	VSIUnlink(m_poGeolocTmpDataset->GetDescription());
223
224	0	auto poXBand = m_poGeolocTmpDataset->GetRasterBand(1);
225	0	auto poYBand = m_poGeolocTmpDataset->GetRasterBand(2);
226
227		// Case of regular grid.
228		// The XBAND contains the x coordinates for all lines.
229		// The YBAND contains the y coordinates for all columns.
230
231	0	double *padfTempX =
232	0	static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
233	0	double *padfTempY =
234	0	static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
235	0	if (padfTempX == nullptr \|\| padfTempY == nullptr)
236	0	{
237	0	CPLFree(padfTempX);
238	0	CPLFree(padfTempY);
239	0	return false;
240	0	}
241
242	0	CPLErr eErr =
243	0	GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
244	0	padfTempX, nXSize, 1, GDT_Float64, 0, 0);
245
246	0	for (int j = 0; j < nYSize; j++)
247	0	{
248	0	if (poXBand->RasterIO(GF_Write, 0, j, nXSize, 1, padfTempX, nXSize,
249	0	1, GDT_Float64, 0, 0, nullptr) != CE_None)
250	0	{
251	0	eErr = CE_Failure;
252	0	break;
253	0	}
254	0	}
255
256	0	if (eErr == CE_None)
257	0	{
258	0	eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
259	0	1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
260
261	0	for (int i = 0; i < nXSize; i++)
262	0	{
263	0	if (poYBand->RasterIO(GF_Write, i, 0, 1, nYSize, padfTempY, 1,
264	0	nYSize, GDT_Float64, 0, 0,
265	0	nullptr) != CE_None)
266	0	{
267	0	eErr = CE_Failure;
268	0	break;
269	0	}
270	0	}
271	0	}
272
273	0	CPLFree(padfTempX);
274	0	CPLFree(padfTempY);
275
276	0	if (eErr != CE_None)
277	0	return false;
278
279	0	geolocXAccessor.SetBand(poXBand);
280	0	geolocYAccessor.SetBand(poYBand);
281	0	}
282	0	else
283	0	{
284	0	geolocXAccessor.SetBand(
285	0	GDALRasterBand::FromHandle(m_psTransform->hBand_X));
286	0	geolocYAccessor.SetBand(
287	0	GDALRasterBand::FromHandle(m_psTransform->hBand_Y));
288	0	}
289
290	0	GDALGeoLoc<GDALGeoLocDatasetAccessors>::LoadGeolocFinish(m_psTransform);
291	0	return true;
292	0	}
293
294		/! @endcond /