1 /*********************************************************************
2 * convolve.c
3 *********************************************************************/
4
5 /* Standard includes */
6 #include <math.h>
7 #include <stdlib.h> /* malloc(), realloc() */
8
9 /* Our includes */
10 #include "base.h"
11 #include "error.h"
12 #include "convolve.h"
13 #include "klt_util.h" /* printing */
14
15 #define MAX_KERNEL_WIDTH 71
16
17
18 typedef struct {
19 int width;
20 float data[MAX_KERNEL_WIDTH];
21 } ConvolutionKernel;
22
23 /* Kernels */
24 ConvolutionKernel gauss_kernel;
25 ConvolutionKernel gaussderiv_kernel;
26 float sigma_last = -10.0;
27
28
29 /*********************************************************************
30 * _KLTToFloatImage
31 *
32 * Given a pointer to image data (probably unsigned chars), copy
33 * data to a float image.
34 */
35
_KLTToFloatImage(KLT_PixelType * img,int ncols,int nrows,_KLT_FloatImage floatimg)36 void _KLTToFloatImage(
37 KLT_PixelType *img,
38 int ncols, int nrows,
39 _KLT_FloatImage floatimg)
40 {
41 KLT_PixelType *ptrend = img + ncols*nrows;
42 float *ptrout = floatimg->data;
43
44 floatimg->ncols = ncols;
45 floatimg->nrows = nrows;
46
47 while (img < ptrend) *ptrout++ = (float) *img++;
48 }
49
50
51 /*********************************************************************
52 * _computeKernels
53 */
54
_computeKernels(float sigma,ConvolutionKernel * gauss,ConvolutionKernel * gaussderiv)55 void _computeKernels(
56 float sigma,
57 ConvolutionKernel *gauss,
58 ConvolutionKernel *gaussderiv)
59 {
60 const float factor = 0.01f; /* for truncating tail */
61 int i;
62
63 /* Compute kernels, and automatically determine widths */
64 {
65 const int hw = MAX_KERNEL_WIDTH / 2;
66 float max_gauss = 1.0f, max_gaussderiv = (float) (sigma*exp(-0.5f));
67
68 /* Compute gauss and deriv */
69 for (i = -hw ; i <= hw ; i++) {
70 gauss->data[i+hw] = (float) exp(-i*i / (2*sigma*sigma));
71 gaussderiv->data[i+hw] = -i * gauss->data[i+hw];
72 }
73
74 /* Compute widths */
75 gauss->width = MAX_KERNEL_WIDTH;
76 for (i = -hw ; fabs(gauss->data[i+hw] / max_gauss) < factor ;
77 i++, gauss->width -= 2);
78 gaussderiv->width = MAX_KERNEL_WIDTH;
79 for (i = -hw ; fabs(gaussderiv->data[i+hw] / max_gaussderiv) < factor ;
80 i++, gaussderiv->width -= 2);
81 if (gauss->width == MAX_KERNEL_WIDTH ||
82 gaussderiv->width == MAX_KERNEL_WIDTH)
83 KLTError("(_computeKernels) MAX_KERNEL_WIDTH %d is too small for "
84 "a sigma of %f", MAX_KERNEL_WIDTH, sigma);
85 }
86
87 /* Shift if width less than MAX_KERNEL_WIDTH */
88 for (i = 0 ; i < gauss->width ; i++)
89 gauss->data[i] = gauss->data[i+(MAX_KERNEL_WIDTH-gauss->width)/2];
90 for (i = 0 ; i < gaussderiv->width ; i++)
91 gaussderiv->data[i] = gaussderiv->data[i+(MAX_KERNEL_WIDTH-gaussderiv->width)/2];
92 /* Normalize gauss and deriv */
93 {
94 const int hw = gaussderiv->width / 2;
95 float den;
96
97 den = 0.0;
98 for (i = 0 ; i < gauss->width ; i++) den += gauss->data[i];
99 for (i = 0 ; i < gauss->width ; i++) gauss->data[i] /= den;
100 den = 0.0;
101 for (i = -hw ; i <= hw ; i++) den -= i*gaussderiv->data[i+hw];
102 for (i = -hw ; i <= hw ; i++) gaussderiv->data[i+hw] /= den;
103 }
104
105 sigma_last = sigma;
106 }
107
108
109 /*********************************************************************
110 * _KLTGetKernelWidths
111 *
112 */
113
_KLTGetKernelWidths(float sigma,int * gauss_width,int * gaussderiv_width)114 void _KLTGetKernelWidths(
115 float sigma,
116 int *gauss_width,
117 int *gaussderiv_width)
118 {
119 _computeKernels(sigma, &gauss_kernel, &gaussderiv_kernel);
120 *gauss_width = gauss_kernel.width;
121 *gaussderiv_width = gaussderiv_kernel.width;
122 }
123
124
125 /*********************************************************************
126 * _convolveImageHoriz
127 */
128
_convolveImageHoriz(_KLT_FloatImage imgin,ConvolutionKernel kernel,_KLT_FloatImage imgout)129 void _convolveImageHoriz(
130 _KLT_FloatImage imgin,
131 ConvolutionKernel kernel,
132 _KLT_FloatImage imgout)
133 {
134 float *ptrrow = imgin->data; /* Points to row's first pixel */
135 float *ptrout = imgout->data, /* Points to next output pixel */
136 *ppp;
137 float sum;
138 int radius = kernel.width / 2;
139 int ncols = imgin->ncols, nrows = imgin->nrows;
140 int i, j, k;
141
142 /* For each row, do ... */
143 for (j = 0 ; j < nrows ; j++) {
144
145 /* Zero leftmost columns */
146 for (i = 0 ; i < radius ; i++)
147 *ptrout++ = 0.0;
148
149 /* Convolve middle columns with kernel */
150 for ( ; i < ncols - radius ; i++) {
151 ppp = ptrrow + i - radius;
152 sum = 0.0;
153 for (k = kernel.width-1 ; k >= 0 ; k--)
154 sum += *ppp++ * kernel.data[k];
155 *ptrout++ = sum;
156 }
157
158 /* Zero rightmost columns */
159 for ( ; i < ncols ; i++)
160 *ptrout++ = 0.0;
161
162 ptrrow += ncols;
163 }
164 }
165
166
167 /*********************************************************************
168 * _convolveImageVert
169 */
170
_convolveImageVert(_KLT_FloatImage imgin,ConvolutionKernel kernel,_KLT_FloatImage imgout)171 void _convolveImageVert(
172 _KLT_FloatImage imgin,
173 ConvolutionKernel kernel,
174 _KLT_FloatImage imgout)
175 {
176 float *ptrcol = imgin->data; /* Points to row's first pixel */
177 float *ptrout = imgout->data, /* Points to next output pixel */
178 *ppp;
179 float sum;
180 int radius = kernel.width / 2;
181 int ncols = imgin->ncols, nrows = imgin->nrows;
182 int i, j, k;
183
184 /* For each column, do ... */
185 for (i = 0 ; i < ncols ; i++) {
186
187 /* Zero topmost rows */
188 for (j = 0 ; j < radius ; j++) {
189 *ptrout = 0.0;
190 ptrout += ncols;
191 }
192
193 /* Convolve middle rows with kernel */
194 for ( ; j < nrows - radius ; j++) {
195 ppp = ptrcol + ncols * (j - radius);
196 sum = 0.0;
197 for (k = kernel.width-1 ; k >= 0 ; k--) {
198 sum += *ppp * kernel.data[k];
199 ppp += ncols;
200 }
201 *ptrout = sum;
202 ptrout += ncols;
203 }
204
205 /* Zero bottommost rows */
206 for ( ; j < nrows ; j++) {
207 *ptrout = 0.0;
208 ptrout += ncols;
209 }
210
211 ptrcol++;
212 ptrout -= nrows * ncols - 1;
213 }
214 }
215
216
217 /*********************************************************************
218 * _convolveSeparate
219 */
220
_convolveSeparate(_KLT_FloatImage imgin,ConvolutionKernel horiz_kernel,ConvolutionKernel vert_kernel,_KLT_FloatImage imgout)221 void _convolveSeparate(
222 _KLT_FloatImage imgin,
223 ConvolutionKernel horiz_kernel,
224 ConvolutionKernel vert_kernel,
225 _KLT_FloatImage imgout)
226 {
227 /* Create temporary image */
228 _KLT_FloatImage tmpimg;
229 tmpimg = _KLTCreateFloatImage(imgin->ncols, imgin->nrows);
230
231 /* Do convolution */
232 _convolveImageHoriz(imgin, horiz_kernel, tmpimg);
233
234 _convolveImageVert(tmpimg, vert_kernel, imgout);
235
236 /* Free memory */
237 _KLTFreeFloatImage(tmpimg);
238 }
239
240
241 /*********************************************************************
242 * _KLTComputeGradients
243 */
244
_KLTComputeGradients(_KLT_FloatImage img,float sigma,_KLT_FloatImage gradx,_KLT_FloatImage grady)245 void _KLTComputeGradients(
246 _KLT_FloatImage img,
247 float sigma,
248 _KLT_FloatImage gradx,
249 _KLT_FloatImage grady)
250 {
251
252 /* Compute kernels, if necessary */
253 if (fabs(sigma - sigma_last) > 0.05)
254 _computeKernels(sigma, &gauss_kernel, &gaussderiv_kernel);
255
256 _convolveSeparate(img, gaussderiv_kernel, gauss_kernel, gradx);
257 _convolveSeparate(img, gauss_kernel, gaussderiv_kernel, grady);
258
259 }
260
261
262 /*********************************************************************
263 * _KLTComputeSmoothedImage
264 */
265
_KLTComputeSmoothedImage(_KLT_FloatImage img,float sigma,_KLT_FloatImage smooth)266 void _KLTComputeSmoothedImage(
267 _KLT_FloatImage img,
268 float sigma,
269 _KLT_FloatImage smooth)
270 {
271 /* Compute kernel, if necessary; gauss_deriv is not used */
272 if (fabs(sigma - sigma_last) > 0.05)
273 _computeKernels(sigma, &gauss_kernel, &gaussderiv_kernel);
274
275 _convolveSeparate(img, gauss_kernel, gauss_kernel, smooth);
276 }
277
278
279
280