/src/gdal/alg/gdalgeoloc_carray_accessor.h

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

/*! @cond Doxygen_Suppress */

/************************************************************************/
/*                        GDALGeoLocCArrayAccessors                     */
/************************************************************************/

class GDALGeoLocCArrayAccessors
{
    typedef class GDALGeoLocCArrayAccessors AccessorType;

    GDALGeoLocTransformInfo *m_psTransform;
    double *m_padfGeoLocX = nullptr;
    double *m_padfGeoLocY = nullptr;
    float *m_pafBackMapX = nullptr;
    float *m_pafBackMapY = nullptr;
    float *m_wgtsBackMap = nullptr;

    bool LoadGeoloc(bool bIsRegularGrid);

  public:
    template <class Type> struct CArrayAccessor
    {
        Type *m_array;
        size_t m_nXSize;

        CArrayAccessor(Type *array, size_t nXSize)
            : m_array(array), m_nXSize(nXSize)
        {
        }

        inline Type Get(int nX, int nY, bool *pbSuccess = nullptr)
        {
            if (pbSuccess)
                *pbSuccess = true;
            return m_array[nY * m_nXSize + nX];
        }

        inline bool Set(int nX, int nY, Type val)
        {
            m_array[nY * m_nXSize + nX] = val;
            return true;
        }
    };

    CArrayAccessor<double> geolocXAccessor;
    CArrayAccessor<double> geolocYAccessor;
    CArrayAccessor<float> backMapXAccessor;
    CArrayAccessor<float> backMapYAccessor;
    CArrayAccessor<float> backMapWeightAccessor;

    explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo *psTransform)
        : m_psTransform(psTransform), geolocXAccessor(nullptr, 0),
          geolocYAccessor(nullptr, 0), backMapXAccessor(nullptr, 0),
          backMapYAccessor(nullptr, 0), backMapWeightAccessor(nullptr, 0)
    {
    }

    ~GDALGeoLocCArrayAccessors()
    {
        VSIFree(m_pafBackMapX);
        VSIFree(m_pafBackMapY);
        VSIFree(m_padfGeoLocX);
        VSIFree(m_padfGeoLocY);
        VSIFree(m_wgtsBackMap);
    }

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

    bool Load(bool bIsRegularGrid, bool bUseQuadtree);

    bool AllocateBackMap();

    GDALDataset *GetBackmapDataset();

    static void FlushBackmapCaches()
    {
    }

    static void ReleaseBackmapDataset(GDALDataset *poDS)
    {
        delete poDS;
    }

    void FreeWghtsBackMap();
};

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

bool GDALGeoLocCArrayAccessors::AllocateBackMap()
{
    m_pafBackMapX = static_cast<float *>(
        VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
                            m_psTransform->nBackMapHeight, sizeof(float)));
    m_pafBackMapY = static_cast<float *>(
        VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
                            m_psTransform->nBackMapHeight, sizeof(float)));

    m_wgtsBackMap = static_cast<float *>(
        VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
                            m_psTransform->nBackMapHeight, sizeof(float)));

    if (m_pafBackMapX == nullptr || m_pafBackMapY == nullptr ||
        m_wgtsBackMap == nullptr)
    {
        return false;
    }

    const size_t nBMXYCount =
        static_cast<size_t>(m_psTransform->nBackMapWidth) *
        m_psTransform->nBackMapHeight;
    for (size_t i = 0; i < nBMXYCount; i++)
    {
        m_pafBackMapX[i] = 0;
        m_pafBackMapY[i] = 0;
        m_wgtsBackMap[i] = 0.0;
    }

    backMapXAccessor.m_array = m_pafBackMapX;
    backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;

    backMapYAccessor.m_array = m_pafBackMapY;
    backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;

    backMapWeightAccessor.m_array = m_wgtsBackMap;
    backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;

    return true;
}

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

void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
{
    VSIFree(m_wgtsBackMap);
    m_wgtsBackMap = nullptr;
    backMapWeightAccessor.m_array = nullptr;
    backMapWeightAccessor.m_nXSize = 0;
}

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

GDALDataset *GDALGeoLocCArrayAccessors::GetBackmapDataset()
{
    auto poMEMDS = MEMDataset::Create("", m_psTransform->nBackMapWidth,
                                      m_psTransform->nBackMapHeight, 0,
                                      GDT_Float32, nullptr);

    for (int i = 1; i <= 2; i++)
    {
        void *ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
        GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
            poMEMDS, i, static_cast<GByte *>(ptr), GDT_Float32, 0, 0, false);
        poMEMDS->AddMEMBand(hMEMBand);
        poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
    }
    return poMEMDS;
}

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

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

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

bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)

{
    const int nXSize = m_psTransform->nGeoLocXSize;
    const int nYSize = m_psTransform->nGeoLocYSize;

    m_padfGeoLocY = static_cast<double *>(
        VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
    m_padfGeoLocX = static_cast<double *>(
        VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));

    if (m_padfGeoLocX == nullptr || m_padfGeoLocY == nullptr)
    {
        return false;
    }

    if (bIsRegularGrid)
    {
        // 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 (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
        {
            memcpy(m_padfGeoLocX + j * nXSize, padfTempX,
                   nXSize * sizeof(double));
        }

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

            for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
            {
                for (size_t i = 0; i < static_cast<size_t>(nXSize); i++)
                {
                    m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
                }
            }
        }

        CPLFree(padfTempX);
        CPLFree(padfTempY);

        if (eErr != CE_None)
            return false;
    }
    else
    {
        if (GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, nYSize,
                         m_padfGeoLocX, nXSize, nYSize, GDT_Float64, 0,
                         0) != CE_None ||
            GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nXSize, nYSize,
                         m_padfGeoLocY, nXSize, nYSize, GDT_Float64, 0,
                         0) != CE_None)
            return false;
    }

    geolocXAccessor.m_array = m_padfGeoLocX;
    geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;

    geolocYAccessor.m_array = m_padfGeoLocY;
    geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;

    GDALGeoLoc<GDALGeoLocCArrayAccessors>::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: C-Array 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		/! @cond Doxygen_Suppress /
14
15		/************************************************************************/
16		/* GDALGeoLocCArrayAccessors */
17		/************************************************************************/
18
19		class GDALGeoLocCArrayAccessors
20		{
21		typedef class GDALGeoLocCArrayAccessors AccessorType;
22
23		GDALGeoLocTransformInfo *m_psTransform;
24		double *m_padfGeoLocX = nullptr;
25		double *m_padfGeoLocY = nullptr;
26		float *m_pafBackMapX = nullptr;
27		float *m_pafBackMapY = nullptr;
28		float *m_wgtsBackMap = nullptr;
29
30		bool LoadGeoloc(bool bIsRegularGrid);
31
32		public:
33		template <class Type> struct CArrayAccessor
34		{
35		Type *m_array;
36		size_t m_nXSize;
37
38		CArrayAccessor(Type *array, size_t nXSize)
39	0	: m_array(array), m_nXSize(nXSize)
40	0	{
41	0	} Unexecuted instantiation: GDALGeoLocCArrayAccessors::CArrayAccessor<double>::CArrayAccessor(double, unsigned long) Unexecuted instantiation: GDALGeoLocCArrayAccessors::CArrayAccessor<float>::CArrayAccessor(float, unsigned long)
42
43		inline Type Get(int nX, int nY, bool *pbSuccess = nullptr)
44	0	{
45	0	if (pbSuccess)
46	0	*pbSuccess = true;
47	0	return m_array[nY * m_nXSize + nX];
48	0	} Unexecuted instantiation: GDALGeoLocCArrayAccessors::CArrayAccessor<double>::Get(int, int, bool) Unexecuted instantiation: GDALGeoLocCArrayAccessors::CArrayAccessor<float>::Get(int, int, bool)
49
50		inline bool Set(int nX, int nY, Type val)
51	0	{
52	0	m_array[nY * m_nXSize + nX] = val;
53	0	return true;
54	0	}
55		};
56
57		CArrayAccessor<double> geolocXAccessor;
58		CArrayAccessor<double> geolocYAccessor;
59		CArrayAccessor<float> backMapXAccessor;
60		CArrayAccessor<float> backMapYAccessor;
61		CArrayAccessor<float> backMapWeightAccessor;
62
63		explicit GDALGeoLocCArrayAccessors(GDALGeoLocTransformInfo *psTransform)
64	0	: m_psTransform(psTransform), geolocXAccessor(nullptr, 0),
65	0	geolocYAccessor(nullptr, 0), backMapXAccessor(nullptr, 0),
66	0	backMapYAccessor(nullptr, 0), backMapWeightAccessor(nullptr, 0)
67	0	{
68	0	}
69
70		~GDALGeoLocCArrayAccessors()
71	0	{
72	0	VSIFree(m_pafBackMapX);
73	0	VSIFree(m_pafBackMapY);
74	0	VSIFree(m_padfGeoLocX);
75	0	VSIFree(m_padfGeoLocY);
76	0	VSIFree(m_wgtsBackMap);
77	0	}
78
79		GDALGeoLocCArrayAccessors(const GDALGeoLocCArrayAccessors &) = delete;
80		GDALGeoLocCArrayAccessors &
81		operator=(const GDALGeoLocCArrayAccessors &) = delete;
82
83		bool Load(bool bIsRegularGrid, bool bUseQuadtree);
84
85		bool AllocateBackMap();
86
87		GDALDataset *GetBackmapDataset();
88
89		static void FlushBackmapCaches()
90	0	{
91	0	}
92
93		static void ReleaseBackmapDataset(GDALDataset *poDS)
94	0	{
95	0	delete poDS;
96	0	}
97
98		void FreeWghtsBackMap();
99		};
100
101		/************************************************************************/
102		/* AllocateBackMap() */
103		/************************************************************************/
104
105		bool GDALGeoLocCArrayAccessors::AllocateBackMap()
106	0	{
107	0	m_pafBackMapX = static_cast<float *>(
108	0	VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
109	0	m_psTransform->nBackMapHeight, sizeof(float)));
110	0	m_pafBackMapY = static_cast<float *>(
111	0	VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
112	0	m_psTransform->nBackMapHeight, sizeof(float)));
113
114	0	m_wgtsBackMap = static_cast<float *>(
115	0	VSI_MALLOC3_VERBOSE(m_psTransform->nBackMapWidth,
116	0	m_psTransform->nBackMapHeight, sizeof(float)));
117
118	0	if (m_pafBackMapX == nullptr \|\| m_pafBackMapY == nullptr \|\|
119	0	m_wgtsBackMap == nullptr)
120	0	{
121	0	return false;
122	0	}
123
124	0	const size_t nBMXYCount =
125	0	static_cast<size_t>(m_psTransform->nBackMapWidth) *
126	0	m_psTransform->nBackMapHeight;
127	0	for (size_t i = 0; i < nBMXYCount; i++)
128	0	{
129	0	m_pafBackMapX[i] = 0;
130	0	m_pafBackMapY[i] = 0;
131	0	m_wgtsBackMap[i] = 0.0;
132	0	}
133
134	0	backMapXAccessor.m_array = m_pafBackMapX;
135	0	backMapXAccessor.m_nXSize = m_psTransform->nBackMapWidth;
136
137	0	backMapYAccessor.m_array = m_pafBackMapY;
138	0	backMapYAccessor.m_nXSize = m_psTransform->nBackMapWidth;
139
140	0	backMapWeightAccessor.m_array = m_wgtsBackMap;
141	0	backMapWeightAccessor.m_nXSize = m_psTransform->nBackMapWidth;
142
143	0	return true;
144	0	}
145
146		/************************************************************************/
147		/* FreeWghtsBackMap() */
148		/************************************************************************/
149
150		void GDALGeoLocCArrayAccessors::FreeWghtsBackMap()
151	0	{
152	0	VSIFree(m_wgtsBackMap);
153	0	m_wgtsBackMap = nullptr;
154	0	backMapWeightAccessor.m_array = nullptr;
155	0	backMapWeightAccessor.m_nXSize = 0;
156	0	}
157
158		/************************************************************************/
159		/* GetBackmapDataset() */
160		/************************************************************************/
161
162		GDALDataset *GDALGeoLocCArrayAccessors::GetBackmapDataset()
163	0	{
164	0	auto poMEMDS = MEMDataset::Create("", m_psTransform->nBackMapWidth,
165	0	m_psTransform->nBackMapHeight, 0,
166	0	GDT_Float32, nullptr);
167
168	0	for (int i = 1; i <= 2; i++)
169	0	{
170	0	void *ptr = (i == 1) ? m_pafBackMapX : m_pafBackMapY;
171	0	GDALRasterBandH hMEMBand = MEMCreateRasterBandEx(
172	0	poMEMDS, i, static_cast<GByte *>(ptr), GDT_Float32, 0, 0, false);
173	0	poMEMDS->AddMEMBand(hMEMBand);
174	0	poMEMDS->GetRasterBand(i)->SetNoDataValue(INVALID_BMXY);
175	0	}
176	0	return poMEMDS;
177	0	}
178
179		/************************************************************************/
180		/* Load() */
181		/************************************************************************/
182
183		bool GDALGeoLocCArrayAccessors::Load(bool bIsRegularGrid, bool bUseQuadtree)
184	0	{
185	0	return LoadGeoloc(bIsRegularGrid) &&
186	0	((bUseQuadtree && GDALGeoLocBuildQuadTree(m_psTransform)) \|\|
187	0	(!bUseQuadtree &&
188	0	GDALGeoLoc<AccessorType>::GenerateBackMap(m_psTransform)));
189	0	}
190
191		/************************************************************************/
192		/* LoadGeoloc() */
193		/************************************************************************/
194
195		bool GDALGeoLocCArrayAccessors::LoadGeoloc(bool bIsRegularGrid)
196
197	0	{
198	0	const int nXSize = m_psTransform->nGeoLocXSize;
199	0	const int nYSize = m_psTransform->nGeoLocYSize;
200
201	0	m_padfGeoLocY = static_cast<double *>(
202	0	VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
203	0	m_padfGeoLocX = static_cast<double *>(
204	0	VSI_MALLOC3_VERBOSE(sizeof(double), nXSize, nYSize));
205
206	0	if (m_padfGeoLocX == nullptr \|\| m_padfGeoLocY == nullptr)
207	0	{
208	0	return false;
209	0	}
210
211	0	if (bIsRegularGrid)
212	0	{
213		// Case of regular grid.
214		// The XBAND contains the x coordinates for all lines.
215		// The YBAND contains the y coordinates for all columns.
216
217	0	double *padfTempX =
218	0	static_cast<double *>(VSI_MALLOC2_VERBOSE(nXSize, sizeof(double)));
219	0	double *padfTempY =
220	0	static_cast<double *>(VSI_MALLOC2_VERBOSE(nYSize, sizeof(double)));
221	0	if (padfTempX == nullptr \|\| padfTempY == nullptr)
222	0	{
223	0	CPLFree(padfTempX);
224	0	CPLFree(padfTempY);
225	0	return false;
226	0	}
227
228	0	CPLErr eErr =
229	0	GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, 1,
230	0	padfTempX, nXSize, 1, GDT_Float64, 0, 0);
231
232	0	for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
233	0	{
234	0	memcpy(m_padfGeoLocX + j * nXSize, padfTempX,
235	0	nXSize * sizeof(double));
236	0	}
237
238	0	if (eErr == CE_None)
239	0	{
240	0	eErr = GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nYSize,
241	0	1, padfTempY, nYSize, 1, GDT_Float64, 0, 0);
242
243	0	for (size_t j = 0; j < static_cast<size_t>(nYSize); j++)
244	0	{
245	0	for (size_t i = 0; i < static_cast<size_t>(nXSize); i++)
246	0	{
247	0	m_padfGeoLocY[j * nXSize + i] = padfTempY[j];
248	0	}
249	0	}
250	0	}
251
252	0	CPLFree(padfTempX);
253	0	CPLFree(padfTempY);
254
255	0	if (eErr != CE_None)
256	0	return false;
257	0	}
258	0	else
259	0	{
260	0	if (GDALRasterIO(m_psTransform->hBand_X, GF_Read, 0, 0, nXSize, nYSize,
261	0	m_padfGeoLocX, nXSize, nYSize, GDT_Float64, 0,
262	0	0) != CE_None \|\|
263	0	GDALRasterIO(m_psTransform->hBand_Y, GF_Read, 0, 0, nXSize, nYSize,
264	0	m_padfGeoLocY, nXSize, nYSize, GDT_Float64, 0,
265	0	0) != CE_None)
266	0	return false;
267	0	}
268
269	0	geolocXAccessor.m_array = m_padfGeoLocX;
270	0	geolocXAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
271
272	0	geolocYAccessor.m_array = m_padfGeoLocY;
273	0	geolocYAccessor.m_nXSize = m_psTransform->nGeoLocXSize;
274
275	0	GDALGeoLoc<GDALGeoLocCArrayAccessors>::LoadGeolocFinish(m_psTransform);
276	0	return true;
277	0	}
278
279		/! @endcond /