Line data Source code
1 : /******************************************************************************
2 : * Project: GDAL
3 : * Purpose: Correlator
4 : * Author: Andrew Migal, migal.drew@gmail.com
5 : *
6 : ******************************************************************************
7 : * Copyright (c) 2012, Andrew Migal
8 : *
9 : * SPDX-License-Identifier: MIT
10 : ****************************************************************************/
11 :
12 : #include "gdal_simplesurf.h"
13 :
14 : /************************************************************************/
15 : /* ==================================================================== */
16 : /* GDALIntegralImage */
17 : /* ==================================================================== */
18 : /************************************************************************/
19 :
20 : GDALIntegralImage::GDALIntegralImage() = default;
21 :
22 0 : int GDALIntegralImage::GetHeight()
23 : {
24 0 : return nHeight;
25 : }
26 :
27 0 : int GDALIntegralImage::GetWidth()
28 : {
29 0 : return nWidth;
30 : }
31 :
32 0 : void GDALIntegralImage::Initialize(const double **padfImg, int nHeightIn,
33 : int nWidthIn)
34 : {
35 0 : if (pMatrix)
36 : {
37 0 : for (int i = 0; i < nHeight; i++)
38 0 : delete[] pMatrix[i];
39 0 : delete[] pMatrix;
40 : }
41 :
42 : // Memory allocation.
43 0 : pMatrix = new double *[nHeightIn];
44 0 : for (int i = 0; i < nHeightIn; i++)
45 0 : pMatrix[i] = new double[nWidthIn];
46 :
47 0 : nHeight = nHeightIn;
48 0 : nWidth = nWidthIn;
49 :
50 : // Integral image calculation.
51 0 : for (int i = 0; i < nHeight; i++)
52 0 : for (int j = 0; j < nWidth; j++)
53 : {
54 0 : const double val = padfImg[i][j];
55 0 : double a = 0.0;
56 0 : double b = 0.0;
57 0 : double c = 0.0;
58 :
59 0 : if (i - 1 >= 0 && j - 1 >= 0)
60 0 : a = pMatrix[i - 1][j - 1];
61 0 : if (j - 1 >= 0)
62 0 : b = pMatrix[i][j - 1];
63 0 : if (i - 1 >= 0)
64 0 : c = pMatrix[i - 1][j];
65 :
66 : // New value based on previous calculations.
67 0 : pMatrix[i][j] = val - a + b + c;
68 : }
69 0 : }
70 :
71 : /*
72 : * Returns value of specified cell.
73 : */
74 0 : double GDALIntegralImage::GetValue(int nRow, int nCol)
75 : {
76 0 : if (!((nRow >= 0 && nRow < nHeight) && (nCol >= 0 && nCol < nWidth)))
77 0 : return 0;
78 :
79 0 : return pMatrix[nRow][nCol];
80 : }
81 :
82 0 : double GDALIntegralImage::GetRectangleSum(int nRow, int nCol, int nWidthIn,
83 : int nHeightIn)
84 : {
85 : // Left top point of rectangle is first.
86 0 : const int w = nWidthIn - 1;
87 0 : const int h = nHeightIn - 1;
88 :
89 0 : const int row = nRow;
90 0 : const int col = nCol;
91 :
92 : // Left top point.
93 0 : const int lt_row = (row <= nHeight) ? (row - 1) : -1;
94 0 : const int lt_col = (col <= nWidth) ? (col - 1) : -1;
95 : // Right bottom point of the rectangle.
96 0 : const int rb_row = (row + h < nHeight) ? (row + h) : (nHeight - 1);
97 0 : const int rb_col = (col + w < nWidth) ? (col + w) : (nWidth - 1);
98 :
99 0 : double a = 0.0;
100 0 : double b = 0.0;
101 0 : double c = 0.0;
102 0 : double d = 0.0;
103 :
104 0 : if (lt_row >= 0 && lt_col >= 0)
105 0 : a = GetValue(lt_row, lt_col);
106 :
107 0 : if (lt_row >= 0 && rb_col >= 0)
108 0 : b = GetValue(lt_row, rb_col);
109 :
110 0 : if (rb_row >= 0 && rb_col >= 0)
111 0 : c = GetValue(rb_row, rb_col);
112 :
113 0 : if (rb_row >= 0 && lt_col >= 0)
114 0 : d = GetValue(rb_row, lt_col);
115 :
116 0 : const double res = a + c - b - d;
117 :
118 0 : return res > 0 ? res : 0;
119 : }
120 :
121 0 : double GDALIntegralImage::HaarWavelet_X(int nRow, int nCol, int nSize)
122 : {
123 0 : return GetRectangleSum(nRow, nCol + nSize / 2, nSize / 2, nSize) -
124 0 : GetRectangleSum(nRow, nCol, nSize / 2, nSize);
125 : }
126 :
127 0 : double GDALIntegralImage::HaarWavelet_Y(int nRow, int nCol, int nSize)
128 : {
129 0 : return GetRectangleSum(nRow + nSize / 2, nCol, nSize, nSize / 2) -
130 0 : GetRectangleSum(nRow, nCol, nSize, nSize / 2);
131 : }
132 :
133 0 : GDALIntegralImage::~GDALIntegralImage()
134 : {
135 : // Clean up memory.
136 0 : for (int i = 0; i < nHeight; i++)
137 0 : delete[] pMatrix[i];
138 :
139 0 : delete[] pMatrix;
140 0 : }
141 :
142 : /************************************************************************/
143 : /* ==================================================================== */
144 : /* GDALOctaveLayer */
145 : /* ==================================================================== */
146 : /************************************************************************/
147 :
148 0 : GDALOctaveLayer::GDALOctaveLayer(int nOctave, int nInterval)
149 : : octaveNum(nOctave),
150 0 : filterSize(3 * static_cast<int>(pow(2.0, nOctave)) * nInterval + 1),
151 0 : radius((filterSize - 1) / 2), scale(static_cast<int>(pow(2.0, nOctave))),
152 0 : width(0), height(0), detHessians(nullptr), signs(nullptr)
153 : {
154 0 : }
155 :
156 0 : void GDALOctaveLayer::ComputeLayer(GDALIntegralImage *poImg)
157 : {
158 0 : width = poImg->GetWidth();
159 0 : height = poImg->GetHeight();
160 :
161 : // Allocate memory for arrays.
162 0 : detHessians = new double *[height];
163 0 : signs = new int *[height];
164 :
165 0 : for (int i = 0; i < height; i++)
166 : {
167 0 : detHessians[i] = new double[width];
168 0 : signs[i] = new int[width];
169 : }
170 :
171 : // 1/3 of filter side.
172 0 : const int lobe = filterSize / 3;
173 :
174 : // Length of the longer side of the lobe in dxx and dyy filters.
175 0 : const int longPart = 2 * lobe - 1;
176 :
177 0 : const int normalization = filterSize * filterSize;
178 :
179 : // Loop over image pixels.
180 : // Filter should remain into image borders.
181 0 : for (int r = radius; r <= height - radius; r++)
182 0 : for (int c = radius; c <= width - radius; c++)
183 : {
184 : // Values of Fast Hessian filters.
185 : double dxx =
186 0 : poImg->GetRectangleSum(r - lobe + 1, c - radius, filterSize,
187 : longPart) -
188 0 : 3 * poImg->GetRectangleSum(r - lobe + 1, c - (lobe - 1) / 2,
189 0 : lobe, longPart);
190 0 : double dyy = poImg->GetRectangleSum(r - radius, c - lobe - 1,
191 : longPart, filterSize) -
192 0 : 3 * poImg->GetRectangleSum(r - lobe + 1, c - lobe + 1,
193 0 : longPart, lobe);
194 : double dxy =
195 0 : poImg->GetRectangleSum(r - lobe, c - lobe, lobe, lobe) +
196 0 : poImg->GetRectangleSum(r + 1, c + 1, lobe, lobe) -
197 0 : poImg->GetRectangleSum(r - lobe, c + 1, lobe, lobe) -
198 0 : poImg->GetRectangleSum(r + 1, c - lobe, lobe, lobe);
199 :
200 0 : dxx /= normalization;
201 0 : dyy /= normalization;
202 0 : dxy /= normalization;
203 :
204 : // Memorize Hessian values and their signs.
205 0 : detHessians[r][c] = dxx * dyy - 0.9 * 0.9 * dxy * dxy;
206 0 : signs[r][c] = (dxx + dyy >= 0) ? 1 : -1;
207 : }
208 0 : }
209 :
210 0 : GDALOctaveLayer::~GDALOctaveLayer()
211 : {
212 0 : for (int i = 0; i < height; i++)
213 : {
214 0 : delete[] detHessians[i];
215 0 : delete[] signs[i];
216 : }
217 :
218 0 : delete[] detHessians;
219 0 : delete[] signs;
220 0 : }
221 :
222 : /************************************************************************/
223 : /* ==================================================================== */
224 : /* GDALOctaveMap */
225 : /* ==================================================================== */
226 : /************************************************************************/
227 :
228 0 : GDALOctaveMap::GDALOctaveMap(int nOctaveStartIn, int nOctaveEndIn)
229 0 : : pMap(new GDALOctaveLayer **[nOctaveEndIn]), octaveStart(nOctaveStartIn),
230 0 : octaveEnd(nOctaveEndIn)
231 : {
232 0 : for (int i = 0; i < octaveEnd; i++)
233 0 : pMap[i] = new GDALOctaveLayer *[INTERVALS];
234 :
235 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
236 0 : for (int i = 1; i <= INTERVALS; i++)
237 0 : pMap[oct - 1][i - 1] = new GDALOctaveLayer(oct, i);
238 0 : }
239 :
240 0 : void GDALOctaveMap::ComputeMap(GDALIntegralImage *poImg)
241 : {
242 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
243 0 : for (int i = 1; i <= INTERVALS; i++)
244 0 : pMap[oct - 1][i - 1]->ComputeLayer(poImg);
245 0 : }
246 :
247 0 : bool GDALOctaveMap::PointIsExtremum(int row, int col, GDALOctaveLayer *bot,
248 : GDALOctaveLayer *mid, GDALOctaveLayer *top,
249 : double threshold)
250 : {
251 : // Check that point in middle layer has all neighbors.
252 0 : if (row <= top->radius || col <= top->radius ||
253 0 : row + top->radius >= top->height || col + top->radius >= top->width)
254 0 : return false;
255 :
256 0 : const double curPoint = mid->detHessians[row][col];
257 :
258 : // Hessian should be higher than threshold.
259 0 : if (curPoint < threshold)
260 0 : return false;
261 :
262 : // Hessian should be higher than Hessians of all neighbors.
263 0 : for (int i = -1; i <= 1; i++)
264 0 : for (int j = -1; j <= 1; j++)
265 : {
266 0 : const double topPoint = top->detHessians[row + i][col + j];
267 0 : const double midPoint = mid->detHessians[row + i][col + j];
268 0 : const double botPoint = bot->detHessians[row + i][col + j];
269 :
270 0 : if (topPoint >= curPoint || botPoint >= curPoint)
271 0 : return false;
272 :
273 0 : if (i != 0 || j != 0)
274 0 : if (midPoint >= curPoint)
275 0 : return false;
276 : }
277 :
278 0 : return true;
279 : }
280 :
281 0 : GDALOctaveMap::~GDALOctaveMap()
282 : {
283 : // Clean up Octave layers.
284 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
285 0 : for (int i = 0; i < INTERVALS; i++)
286 0 : delete pMap[oct - 1][i];
287 :
288 : // Clean up allocated memory.
289 0 : for (int oct = 0; oct < octaveEnd; oct++)
290 0 : delete[] pMap[oct];
291 :
292 0 : delete[] pMap;
293 0 : }
|
