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 : * Permission is hereby granted, free of charge, to any person obtaining a
10 : * copy of this software and associated documentation files (the "Software"),
11 : * to deal in the Software without restriction, including without limitation
12 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13 : * and/or sell copies of the Software, and to permit persons to whom the
14 : * Software is furnished to do so, subject to the following conditions:
15 : *
16 : * The above copyright notice and this permission notice shall be included
17 : * in all copies or substantial portions of the Software.
18 : *
19 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25 : * DEALINGS IN THE SOFTWARE.
26 : ****************************************************************************/
27 :
28 : #include "gdal_simplesurf.h"
29 :
30 : /************************************************************************/
31 : /* ==================================================================== */
32 : /* GDALIntegralImage */
33 : /* ==================================================================== */
34 : /************************************************************************/
35 :
36 : GDALIntegralImage::GDALIntegralImage() = default;
37 :
38 0 : int GDALIntegralImage::GetHeight()
39 : {
40 0 : return nHeight;
41 : }
42 :
43 0 : int GDALIntegralImage::GetWidth()
44 : {
45 0 : return nWidth;
46 : }
47 :
48 0 : void GDALIntegralImage::Initialize(const double **padfImg, int nHeightIn,
49 : int nWidthIn)
50 : {
51 0 : if (pMatrix)
52 : {
53 0 : for (int i = 0; i < nHeight; i++)
54 0 : delete[] pMatrix[i];
55 0 : delete[] pMatrix;
56 : }
57 :
58 : // Memory allocation.
59 0 : pMatrix = new double *[nHeightIn];
60 0 : for (int i = 0; i < nHeightIn; i++)
61 0 : pMatrix[i] = new double[nWidthIn];
62 :
63 0 : nHeight = nHeightIn;
64 0 : nWidth = nWidthIn;
65 :
66 : // Integral image calculation.
67 0 : for (int i = 0; i < nHeight; i++)
68 0 : for (int j = 0; j < nWidth; j++)
69 : {
70 0 : const double val = padfImg[i][j];
71 0 : double a = 0.0;
72 0 : double b = 0.0;
73 0 : double c = 0.0;
74 :
75 0 : if (i - 1 >= 0 && j - 1 >= 0)
76 0 : a = pMatrix[i - 1][j - 1];
77 0 : if (j - 1 >= 0)
78 0 : b = pMatrix[i][j - 1];
79 0 : if (i - 1 >= 0)
80 0 : c = pMatrix[i - 1][j];
81 :
82 : // New value based on previous calculations.
83 0 : pMatrix[i][j] = val - a + b + c;
84 : }
85 0 : }
86 :
87 : /*
88 : * Returns value of specified cell.
89 : */
90 0 : double GDALIntegralImage::GetValue(int nRow, int nCol)
91 : {
92 0 : if (!((nRow >= 0 && nRow < nHeight) && (nCol >= 0 && nCol < nWidth)))
93 0 : return 0;
94 :
95 0 : return pMatrix[nRow][nCol];
96 : }
97 :
98 0 : double GDALIntegralImage::GetRectangleSum(int nRow, int nCol, int nWidthIn,
99 : int nHeightIn)
100 : {
101 : // Left top point of rectangle is first.
102 0 : const int w = nWidthIn - 1;
103 0 : const int h = nHeightIn - 1;
104 :
105 0 : const int row = nRow;
106 0 : const int col = nCol;
107 :
108 : // Left top point.
109 0 : const int lt_row = (row <= nHeight) ? (row - 1) : -1;
110 0 : const int lt_col = (col <= nWidth) ? (col - 1) : -1;
111 : // Right bottom point of the rectangle.
112 0 : const int rb_row = (row + h < nHeight) ? (row + h) : (nHeight - 1);
113 0 : const int rb_col = (col + w < nWidth) ? (col + w) : (nWidth - 1);
114 :
115 0 : double a = 0.0;
116 0 : double b = 0.0;
117 0 : double c = 0.0;
118 0 : double d = 0.0;
119 :
120 0 : if (lt_row >= 0 && lt_col >= 0)
121 0 : a = GetValue(lt_row, lt_col);
122 :
123 0 : if (lt_row >= 0 && rb_col >= 0)
124 0 : b = GetValue(lt_row, rb_col);
125 :
126 0 : if (rb_row >= 0 && rb_col >= 0)
127 0 : c = GetValue(rb_row, rb_col);
128 :
129 0 : if (rb_row >= 0 && lt_col >= 0)
130 0 : d = GetValue(rb_row, lt_col);
131 :
132 0 : const double res = a + c - b - d;
133 :
134 0 : return res > 0 ? res : 0;
135 : }
136 :
137 0 : double GDALIntegralImage::HaarWavelet_X(int nRow, int nCol, int nSize)
138 : {
139 0 : return GetRectangleSum(nRow, nCol + nSize / 2, nSize / 2, nSize) -
140 0 : GetRectangleSum(nRow, nCol, nSize / 2, nSize);
141 : }
142 :
143 0 : double GDALIntegralImage::HaarWavelet_Y(int nRow, int nCol, int nSize)
144 : {
145 0 : return GetRectangleSum(nRow + nSize / 2, nCol, nSize, nSize / 2) -
146 0 : GetRectangleSum(nRow, nCol, nSize, nSize / 2);
147 : }
148 :
149 0 : GDALIntegralImage::~GDALIntegralImage()
150 : {
151 : // Clean up memory.
152 0 : for (int i = 0; i < nHeight; i++)
153 0 : delete[] pMatrix[i];
154 :
155 0 : delete[] pMatrix;
156 0 : }
157 :
158 : /************************************************************************/
159 : /* ==================================================================== */
160 : /* GDALOctaveLayer */
161 : /* ==================================================================== */
162 : /************************************************************************/
163 :
164 0 : GDALOctaveLayer::GDALOctaveLayer(int nOctave, int nInterval)
165 : : octaveNum(nOctave),
166 0 : filterSize(3 * static_cast<int>(pow(2.0, nOctave)) * nInterval + 1),
167 0 : radius((filterSize - 1) / 2), scale(static_cast<int>(pow(2.0, nOctave))),
168 0 : width(0), height(0), detHessians(nullptr), signs(nullptr)
169 : {
170 0 : }
171 :
172 0 : void GDALOctaveLayer::ComputeLayer(GDALIntegralImage *poImg)
173 : {
174 0 : width = poImg->GetWidth();
175 0 : height = poImg->GetHeight();
176 :
177 : // Allocate memory for arrays.
178 0 : detHessians = new double *[height];
179 0 : signs = new int *[height];
180 :
181 0 : for (int i = 0; i < height; i++)
182 : {
183 0 : detHessians[i] = new double[width];
184 0 : signs[i] = new int[width];
185 : }
186 :
187 : // 1/3 of filter side.
188 0 : const int lobe = filterSize / 3;
189 :
190 : // Length of the longer side of the lobe in dxx and dyy filters.
191 0 : const int longPart = 2 * lobe - 1;
192 :
193 0 : const int normalization = filterSize * filterSize;
194 :
195 : // Loop over image pixels.
196 : // Filter should remain into image borders.
197 0 : for (int r = radius; r <= height - radius; r++)
198 0 : for (int c = radius; c <= width - radius; c++)
199 : {
200 : // Values of Fast Hessian filters.
201 : double dxx =
202 0 : poImg->GetRectangleSum(r - lobe + 1, c - radius, filterSize,
203 : longPart) -
204 0 : 3 * poImg->GetRectangleSum(r - lobe + 1, c - (lobe - 1) / 2,
205 0 : lobe, longPart);
206 0 : double dyy = poImg->GetRectangleSum(r - radius, c - lobe - 1,
207 : longPart, filterSize) -
208 0 : 3 * poImg->GetRectangleSum(r - lobe + 1, c - lobe + 1,
209 0 : longPart, lobe);
210 : double dxy =
211 0 : poImg->GetRectangleSum(r - lobe, c - lobe, lobe, lobe) +
212 0 : poImg->GetRectangleSum(r + 1, c + 1, lobe, lobe) -
213 0 : poImg->GetRectangleSum(r - lobe, c + 1, lobe, lobe) -
214 0 : poImg->GetRectangleSum(r + 1, c - lobe, lobe, lobe);
215 :
216 0 : dxx /= normalization;
217 0 : dyy /= normalization;
218 0 : dxy /= normalization;
219 :
220 : // Memorize Hessian values and their signs.
221 0 : detHessians[r][c] = dxx * dyy - 0.9 * 0.9 * dxy * dxy;
222 0 : signs[r][c] = (dxx + dyy >= 0) ? 1 : -1;
223 : }
224 0 : }
225 :
226 0 : GDALOctaveLayer::~GDALOctaveLayer()
227 : {
228 0 : for (int i = 0; i < height; i++)
229 : {
230 0 : delete[] detHessians[i];
231 0 : delete[] signs[i];
232 : }
233 :
234 0 : delete[] detHessians;
235 0 : delete[] signs;
236 0 : }
237 :
238 : /************************************************************************/
239 : /* ==================================================================== */
240 : /* GDALOctaveMap */
241 : /* ==================================================================== */
242 : /************************************************************************/
243 :
244 0 : GDALOctaveMap::GDALOctaveMap(int nOctaveStartIn, int nOctaveEndIn)
245 0 : : pMap(new GDALOctaveLayer **[nOctaveEndIn]), octaveStart(nOctaveStartIn),
246 0 : octaveEnd(nOctaveEndIn)
247 : {
248 0 : for (int i = 0; i < octaveEnd; i++)
249 0 : pMap[i] = new GDALOctaveLayer *[INTERVALS];
250 :
251 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
252 0 : for (int i = 1; i <= INTERVALS; i++)
253 0 : pMap[oct - 1][i - 1] = new GDALOctaveLayer(oct, i);
254 0 : }
255 :
256 0 : void GDALOctaveMap::ComputeMap(GDALIntegralImage *poImg)
257 : {
258 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
259 0 : for (int i = 1; i <= INTERVALS; i++)
260 0 : pMap[oct - 1][i - 1]->ComputeLayer(poImg);
261 0 : }
262 :
263 0 : bool GDALOctaveMap::PointIsExtremum(int row, int col, GDALOctaveLayer *bot,
264 : GDALOctaveLayer *mid, GDALOctaveLayer *top,
265 : double threshold)
266 : {
267 : // Check that point in middle layer has all neighbors.
268 0 : if (row <= top->radius || col <= top->radius ||
269 0 : row + top->radius >= top->height || col + top->radius >= top->width)
270 0 : return false;
271 :
272 0 : const double curPoint = mid->detHessians[row][col];
273 :
274 : // Hessian should be higher than threshold.
275 0 : if (curPoint < threshold)
276 0 : return false;
277 :
278 : // Hessian should be higher than Hessians of all neighbors.
279 0 : for (int i = -1; i <= 1; i++)
280 0 : for (int j = -1; j <= 1; j++)
281 : {
282 0 : const double topPoint = top->detHessians[row + i][col + j];
283 0 : const double midPoint = mid->detHessians[row + i][col + j];
284 0 : const double botPoint = bot->detHessians[row + i][col + j];
285 :
286 0 : if (topPoint >= curPoint || botPoint >= curPoint)
287 0 : return false;
288 :
289 0 : if (i != 0 || j != 0)
290 0 : if (midPoint >= curPoint)
291 0 : return false;
292 : }
293 :
294 0 : return true;
295 : }
296 :
297 0 : GDALOctaveMap::~GDALOctaveMap()
298 : {
299 : // Clean up Octave layers.
300 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
301 0 : for (int i = 0; i < INTERVALS; i++)
302 0 : delete pMap[oct - 1][i];
303 :
304 : // Clean up allocated memory.
305 0 : for (int oct = 0; oct < octaveEnd; oct++)
306 0 : delete[] pMap[oct];
307 :
308 0 : delete[] pMap;
309 0 : }
|
