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