vil_gauss_reduce.cxx
Go to the documentation of this file.
1 // This is core/vil/algo/vil_gauss_reduce.cxx
2 #include <cmath>
3 #include "vil_gauss_reduce.h"
4 //:
5 // \file
6 // \brief Functions to smooth and sub-sample an image in one direction
7 // \author Tim Cootes
8 
9 #ifdef _MSC_VER
10 # include <vcl_msvc_warnings.h>
11 #endif
12 #include <cassert>
13 #include <vxl_config.h> // for vxl_byte
14 #include <vnl/vnl_erf.h>
15 #include <vnl/vnl_math.h>
16 
17 //: Smooth and subsample single plane src_im in x to produce dest_im
18 // Applies 1-5-8-5-1 filter in x, then samples
19 // every other pixel. Fills [0,(ni+1)/2-1][0,nj-1] elements of dest
20 
21 template <>
22 void vil_gauss_reduce_1plane(const vxl_byte* src_im,
23  unsigned src_ni, unsigned src_nj,
24  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
25  vxl_byte* dest_im,
26  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
27 {
28  vxl_byte* d_row = dest_im;
29  const vxl_byte* s_row = src_im;
30  std::ptrdiff_t sxs2 = s_x_step*2;
31  unsigned ni2 = (src_ni-3)/2;
32  for (unsigned y=0;y<src_nj;++y)
33  {
34  // Set first element of row
35  *d_row = static_cast<vxl_byte>(
36  vnl_math::rnd( 0.071f * s_row[sxs2]
37  + 0.357f * s_row[s_x_step]
38  + 0.572f * s_row[0]));
39 
40  vxl_byte * d = d_row + d_x_step;
41  const vxl_byte* s = s_row + sxs2;
42  for (unsigned x=0;x<ni2;++x)
43  {
44  *d = static_cast<vxl_byte>(
45  vnl_math::rnd( 0.05*s[-sxs2] + 0.05*s[sxs2]
46  + 0.25*s[-s_x_step]+ 0.25*s[s_x_step]
47  + 0.4*s[0]));
48 
49  d += d_x_step;
50  s += sxs2;
51  }
52  // Set last elements of row
53  *d = static_cast<vxl_byte>(
54  vnl_math::rnd( 0.071f * s[-sxs2]
55  + 0.357f * s[-s_x_step]
56  + 0.572f * s[0]));
57 
58  d_row += d_y_step;
59  s_row += s_y_step;
60  }
61 }
62 
63 //: Smooth and subsample single plane src_im in x to produce dest_im
64 // Applies 1-5-8-5-1 filter in x, then samples
65 // every other pixel. Fills [0,(ni+1)/2-1][0,nj-1] elements of dest
66 template <>
67 void vil_gauss_reduce_1plane(const float* src_im,
68  unsigned src_ni, unsigned src_nj,
69  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
70  float* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
71 {
72  float* d_row = dest_im;
73  const float* s_row = src_im;
74  std::ptrdiff_t sxs2 = s_x_step*2;
75  unsigned ni2 = (src_ni-3)/2;
76  for (unsigned y=0;y<src_nj;++y)
77  {
78  // Set first element of row
79  *d_row = 0.071f * s_row[sxs2]
80  + 0.357f * s_row[s_x_step]
81  + 0.572f * s_row[0];
82  float * d = d_row + d_x_step;
83  const float* s = s_row + sxs2;
84  for (unsigned x=0;x<ni2;++x)
85  {
86  *d = 0.05f*(s[-sxs2] + s[sxs2])
87  + 0.25f*(s[-s_x_step]+ s[s_x_step])
88  + 0.40f*s[0];
89 
90  d += d_x_step;
91  s += sxs2;
92  }
93  // Set last elements of row
94  *d = 0.071f * s[-sxs2]
95  + 0.357f * s[-s_x_step]
96  + 0.572f * s[0];
97 
98  d_row += d_y_step;
99  s_row += s_y_step;
100  }
101 }
102 
103 
104 //: Smooth and subsample single plane src_im in x to produce dest_im
105 // Applies 1-5-8-5-1 filter in x, then samples
106 // every other pixel. Fills [0,(ni+1)/2-1][0,nj-1] elements of dest
107 template <>
108 void vil_gauss_reduce_1plane(const double* src_im,
109  unsigned src_ni, unsigned src_nj,
110  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
111  double* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
112 {
113  double* d_row = dest_im;
114  const double* s_row = src_im;
115  std::ptrdiff_t sxs2 = s_x_step*2;
116  unsigned ni2 = (src_ni-3)/2;
117  for (unsigned y=0;y<src_nj;++y)
118  {
119  // Set first element of row
120  *d_row = 0.071 * s_row[sxs2]
121  + 0.357 * s_row[s_x_step]
122  + 0.572 * s_row[0];
123  double * d = d_row + d_x_step;
124  const double* s = s_row + sxs2;
125  for (unsigned x=0;x<ni2;++x)
126  {
127  *d = 0.05f*(s[-sxs2] + s[sxs2])
128  + 0.25f*(s[-s_x_step]+ s[s_x_step])
129  + 0.40f*s[0];
130 
131  d += d_x_step;
132  s += sxs2;
133  }
134  // Set last elements of row
135  *d = 0.071f * s[-sxs2]
136  + 0.357f * s[-s_x_step]
137  + 0.572f * s[0];
138 
139  d_row += d_y_step;
140  s_row += s_y_step;
141  }
142 }
143 
144 //: Smooth and subsample single plane src_im in x to produce dest_im
145 // Applies 1-5-8-5-1 filter in x, then samples
146 // every other pixel. Fills [0,(ni+1)/2-1][0,nj-1] elements of dest
147 template <>
148 void vil_gauss_reduce_1plane(const int* src_im,
149  unsigned src_ni, unsigned src_nj,
150  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
151  int* dest_im,
152  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
153 {
154  int* d_row = dest_im;
155  const int* s_row = src_im;
156  std::ptrdiff_t sxs2 = s_x_step*2;
157  unsigned ni2 = (src_ni-3)/2;
158  for (unsigned y=0;y<src_nj;++y)
159  {
160  // Set first element of row
161  *d_row = vnl_math::rnd( 0.071f * s_row[sxs2]
162  + 0.357f * s_row[s_x_step]
163  + 0.572f * s_row[0]);
164 
165  int * d = d_row + d_x_step;
166  const int* s = s_row + sxs2;
167  for (unsigned x=0;x<ni2;++x)
168  {
169  *d = vnl_math::rnd( 0.05*s[-sxs2] + 0.25*s[-s_x_step] +
170  0.05*s[ sxs2] + 0.25*s[ s_x_step] +
171  0.4 *s[0]);
172 
173  d += d_x_step;
174  s += sxs2;
175  }
176  // Set last elements of row
177  *d = vnl_math::rnd( 0.071f * s[-sxs2]
178  + 0.357f * s[-s_x_step]
179  + 0.572f * s[0]);
180 
181  d_row += d_y_step;
182  s_row += s_y_step;
183  }
184 }
185 
186 //: Smooth and subsample single plane src_im in x to produce dest_im
187 // Applies 1-5-8-5-1 filter in x, then samples
188 // every other pixel. Fills [0,(ni+1)/2-1][0,nj-1] elements of dest
189 template <>
190 void vil_gauss_reduce_1plane(const vxl_int_16* src_im,
191  unsigned src_ni, unsigned src_nj,
192  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
193  vxl_int_16* dest_im,
194  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
195 {
196  vxl_int_16* d_row = dest_im;
197  const vxl_int_16* s_row = src_im;
198  std::ptrdiff_t sxs2 = s_x_step*2;
199  unsigned ni2 = (src_ni-3)/2;
200  for (unsigned y=0;y<src_nj;++y)
201  {
202  // Set first element of row
203  *d_row = static_cast<vxl_int_16>(
204  vnl_math::rnd( 0.071f * s_row[sxs2]
205  + 0.357f * s_row[s_x_step]
206  + 0.572f * s_row[0]));
207 
208  vxl_int_16 * d = d_row + d_x_step;
209  const vxl_int_16* s = s_row + sxs2;
210  for (unsigned x=0;x<ni2;++x)
211  {
212  // The 0.5 offset in the following ensures rounding
213  *d = vxl_int_16(0.5 + 0.05*s[-sxs2] + 0.25*s[-s_x_step]
214  + 0.05*s[ sxs2] + 0.25*s[ s_x_step]
215  + 0.4 *s[0]);
216 
217  d += d_x_step;
218  s += sxs2;
219  }
220  // Set last elements of row
221  *d = static_cast<vxl_int_16>(
222  vnl_math::rnd( 0.071f * s[-sxs2]
223  + 0.357f * s[-s_x_step]
224  + 0.572f * s[0]));
225 
226  d_row += d_y_step;
227  s_row += s_y_step;
228  }
229 }
230 
231 //: Smooth and subsample single plane src_im in x, result is 2/3rd size
232 // Applies alternate 1-3-1, 1-1 filter in x, then samples
233 // every other pixel. Fills [0,(2*ni+1)/3-1][0,nj-1] elements of dest
234 //
235 // Note, 131 filter only an approximation
236 template <>
237 void vil_gauss_reduce_2_3_1plane(const float* src_im,
238  unsigned src_ni, unsigned src_nj,
239  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
240  float* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
241 {
242  float* d_row = dest_im;
243  const float* s_row = src_im;
244  std::ptrdiff_t sxs2 = s_x_step*2,sxs3 = s_x_step*3;
245  unsigned d_ni = (2*src_ni+1)/3;
246  unsigned d_ni2 = d_ni/2;
247  for (unsigned y=0;y<src_nj;++y)
248  {
249  // Set first elements of row
250  d_row[0] = 0.75f*s_row[0] + 0.25f*s_row[s_x_step];
251  d_row[d_x_step] = 0.5f*s_row[s_x_step] + 0.5f*s_row[sxs2];
252  float * d = d_row + 2*d_x_step;
253  const float* s = s_row + sxs3;
254  for (unsigned x=1;x<d_ni2;++x)
255  {
256  *d = 0.2f*(s[-s_x_step] + s[s_x_step])+0.6f*s[0];
257  d += d_x_step;
258  *d = 0.5f*(s[s_x_step] + s[sxs2]);
259  d += d_x_step;
260  s += sxs3;
261  }
262  // Set last elements of row
263  if (src_ni%3==1) *d=0.75f*s[-s_x_step] + 0.25f*s[0];
264  else
265  if (src_ni%3==2) *d=0.2f*(s[-s_x_step] + s[s_x_step])+0.6f*s[0];
266 
267  d_row += d_y_step;
268  s_row += s_y_step;
269  }
270 }
271 
272 //: Smooth and subsample single plane src_im in x, result is 2/3rd size
273 // Applies alternate 1-3-1, 1-1 filter in x, then samples
274 // every other pixel. Fills [0,(2*ni+1)/3-1][0,nj-1] elements of dest
275 //
276 // Note, 131 filter only an approximation
277 template <>
278 void vil_gauss_reduce_2_3_1plane(const double* src_im,
279  unsigned src_ni, unsigned src_nj,
280  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
281  double* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
282 {
283  double* d_row = dest_im;
284  const double* s_row = src_im;
285  std::ptrdiff_t sxs2 = s_x_step*2,sxs3 = s_x_step*3;
286  unsigned d_ni = (2*src_ni+1)/3;
287  unsigned d_ni2 = d_ni/2;
288  for (unsigned y=0;y<src_nj;++y)
289  {
290  // Set first elements of row
291  d_row[0] = 0.75*s_row[0] + 0.25*s_row[s_x_step];
292  d_row[d_x_step] = 0.5*s_row[s_x_step] + 0.5*s_row[sxs2];
293  double * d = d_row + 2*d_x_step;
294  const double* s = s_row + sxs3;
295  for (unsigned x=1;x<d_ni2;++x)
296  {
297  *d = 0.2*(s[-s_x_step] + s[s_x_step])+0.6*s[0];
298  d += d_x_step;
299  *d = 0.5*(s[s_x_step] + s[sxs2]);
300  d += d_x_step;
301  s += sxs3;
302  }
303  // Set last elements of row
304  if (src_ni%3==1) *d=0.75*s[-s_x_step] + 0.25*s[0];
305  else
306  if (src_ni%3==2) *d=0.2*(s[-s_x_step] + s[s_x_step])+0.6*s[0];
307 
308  d_row += d_y_step;
309  s_row += s_y_step;
310  }
311 }
312 
313 //: Smooth and subsample single plane src_im in x, result is 2/3rd size
314 // Applies alternate 1-3-1, 1-1 filter in x, then samples
315 // every other pixel. Fills [0,(2*ni+1)/3-1][0,nj-1] elements of dest
316 //
317 // Note, 131 filter only an approximation
318 template <>
319 void vil_gauss_reduce_2_3_1plane(const vxl_byte* src_im,
320  unsigned src_ni, unsigned src_nj,
321  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
322  vxl_byte* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
323 {
324  vxl_byte* d_row = dest_im;
325  const vxl_byte* s_row = src_im;
326  std::ptrdiff_t sxs2 = s_x_step*2,sxs3 = s_x_step*3;
327  unsigned d_ni = (2*src_ni+1)/3;
328  unsigned d_ni2 = d_ni/2;
329  for (unsigned y=0;y<src_nj;++y)
330  {
331  // Set first elements of row
332  // The 0.5 offset in the following ensures rounding
333  d_row[0] = vxl_byte(0.5f + 0.75f*s_row[0] + 0.25f*s_row[s_x_step]);
334  d_row[d_x_step] = vxl_byte(0.5f + 0.5f*s_row[s_x_step] + 0.5f*s_row[sxs2]);
335  vxl_byte * d = d_row + 2*d_x_step;
336  const vxl_byte* s = s_row + sxs3;
337  for (unsigned x=1;x<d_ni2;++x)
338  {
339  *d = vxl_byte(0.5f + 0.2f*(s[-s_x_step] + s[s_x_step]) + 0.6f*s[0]);
340  d += d_x_step;
341  *d = vxl_byte(0.5f + 0.5f*(s[s_x_step] + s[sxs2]));
342  d += d_x_step;
343  s += sxs3;
344  }
345  // Set last elements of row
346  if (src_ni%3==1)
347  *d = vxl_byte(0.5f + 0.75f*s[-s_x_step] + 0.25f*s[0]);
348  else if (src_ni%3==2)
349  *d = vxl_byte(0.5f + 0.2f*(s[-s_x_step] + s[s_x_step]) + 0.6f*s[0]);
350 
351  d_row += d_y_step;
352  s_row += s_y_step;
353  }
354 }
355 
356 
357 //: Smooth and subsample single plane src_im in x, result is 2/3rd size
358 // Applies alternate 1-3-1, 1-1 filter in x, then samples
359 // every other pixel. Fills [0,(2*ni+1)/3-1][0,nj-1] elements of dest
360 //
361 // Note, 131 filter only an approximation
362 template <>
363 void vil_gauss_reduce_2_3_1plane(const int* src_im,
364  unsigned src_ni, unsigned src_nj,
365  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
366  int* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
367 {
368  int* d_row = dest_im;
369  const int* s_row = src_im;
370  std::ptrdiff_t sxs2 = s_x_step*2,sxs3 = s_x_step*3;
371  unsigned d_ni = (2*src_ni+1)/3;
372  unsigned d_ni2 = d_ni/2;
373  for (unsigned y=0;y<src_nj;++y)
374  {
375  // Set first elements of row
376  // The 0.5 offset in the following ensures rounding
377  d_row[0] = int(0.5f + 0.75f*s_row[0] + 0.25f*s_row[s_x_step]);
378  d_row[d_x_step] = int(0.5f + 0.5f*s_row[s_x_step] + 0.5f*s_row[sxs2]);
379  int * d = d_row + 2*d_x_step;
380  const int* s = s_row + sxs3;
381  for (unsigned x=1;x<d_ni2;++x)
382  {
383  *d = int(0.5f + 0.2f*(s[-s_x_step] + s[s_x_step]) + 0.6f*s[0]);
384  d += d_x_step;
385  *d = int(0.5f + 0.5f*(s[s_x_step] + s[sxs2]));
386  d += d_x_step;
387  s += sxs3;
388  }
389  // Set last elements of row
390  if (src_ni%3==1)
391  *d = int(0.5f + 0.75f*s[-s_x_step] + 0.25f*s[0]);
392  else if (src_ni%3==2)
393  *d = int(0.5f + 0.2f*(s[-s_x_step] + s[s_x_step]) + 0.6f*s[0]);
394 
395  d_row += d_y_step;
396  s_row += s_y_step;
397  }
398 }
399 
400 //: Smooth and subsample single plane src_im in x, result is 2/3rd size
401 // Applies alternate 1-3-1, 1-1 filter in x, then samples
402 // every other pixel. Fills [0,(2*ni+1)/3-1][0,nj-1] elements of dest
403 //
404 // Note, 131 filter only an approximation
405 template <>
406 void vil_gauss_reduce_2_3_1plane(const vxl_int_16* src_im,
407  unsigned src_ni, unsigned src_nj,
408  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
409  vxl_int_16* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
410 {
411  vxl_int_16* d_row = dest_im;
412  const vxl_int_16* s_row = src_im;
413  std::ptrdiff_t sxs2 = s_x_step*2,sxs3 = s_x_step*3;
414  unsigned d_ni = (2*src_ni+1)/3;
415  unsigned d_ni2 = d_ni/2;
416  for (unsigned y=0;y<src_nj;++y)
417  {
418  // Set first elements of row
419  // The 0.5 offset in the following ensures rounding
420  d_row[0] = vxl_int_16(0.5f + 0.75f*s_row[0] + 0.25f*s_row[s_x_step]);
421  d_row[d_x_step] = vxl_int_16(0.5f + 0.5f*s_row[s_x_step] + 0.5f*s_row[sxs2]);
422  vxl_int_16 * d = d_row + 2*d_x_step;
423  const vxl_int_16* s = s_row + sxs3;
424  for (unsigned x=1;x<d_ni2;++x)
425  {
426  *d = vxl_int_16(0.5f + 0.2f*(s[-s_x_step] + s[s_x_step]) + 0.6f*s[0]);
427  d += d_x_step;
428  *d = vxl_int_16(0.5f + 0.5f*(s[s_x_step] + s[sxs2]));
429  d += d_x_step;
430  s += sxs3;
431  }
432  // Set last elements of row
433  if (src_ni%3==1)
434  *d = vxl_int_16(0.5f + 0.75f*s[-s_x_step] + 0.25f*s[0]);
435  else if (src_ni%3==2)
436  *d = vxl_int_16(0.5f + 0.2f*(s[-s_x_step] + s[s_x_step]) + 0.6f*s[0]);
437 
438  d_row += d_y_step;
439  s_row += s_y_step;
440  }
441 }
442 
443 //: Smooth and subsample single plane src_im in x to produce dest_im using 121 filter in x and y
444 // Smooths with a 3x3 filter and subsamples
445 template <>
446 void vil_gauss_reduce_121_1plane(const vxl_byte* src_im,
447  unsigned src_ni, unsigned src_nj,
448  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
449  vxl_byte* dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
450 {
451  std::ptrdiff_t sxs2 = s_x_step*2;
452  std::ptrdiff_t sys2 = s_y_step*2;
453  vxl_byte* d_row = dest_im+d_y_step;
454  const vxl_byte* s_row1 = src_im + s_y_step;
455  const vxl_byte* s_row2 = s_row1 + s_y_step;
456  const vxl_byte* s_row3 = s_row2 + s_y_step;
457  unsigned ni2 = (src_ni-2)/2;
458  unsigned nj2 = (src_nj-2)/2;
459  for (unsigned y=0;y<nj2;++y)
460  {
461  // Set first element of row
462  *d_row = *s_row2;
463  vxl_byte * d = d_row + d_x_step;
464  const vxl_byte* s1 = s_row1 + sxs2;
465  const vxl_byte* s2 = s_row2 + sxs2;
466  const vxl_byte* s3 = s_row3 + sxs2;
467  for (unsigned x=0;x<ni2;++x)
468  {
469  // The following is a little inefficient - could group terms to reduce arithmetic
470  // Add 0.5 so that truncating effectively rounds
471  *d = vxl_byte( 0.0625f * s1[-s_x_step] + 0.125f * s1[0] + 0.0625f * s1[s_x_step]
472  + 0.1250f * s2[-s_x_step] + 0.250f * s2[0] + 0.1250f * s2[s_x_step]
473  + 0.0625f * s3[-s_x_step] + 0.125f * s3[0] + 0.0625f * s3[s_x_step] +0.5);
474 
475  d += d_x_step;
476  s1 += sxs2;
477  s2 += sxs2;
478  s3 += sxs2;
479  }
480  // Set last elements of row
481  if (src_ni&1)
482  *d = *s2;
483 
484  d_row += d_y_step;
485  s_row1 += sys2;
486  s_row2 += sys2;
487  s_row3 += sys2;
488  }
489 
490  // Need to set first and last rows as well
491 
492  // Dest image should be (src_ni+1)/2 x (src_nj+1)/2
493  const vxl_byte* s0 = src_im;
494  unsigned ni=(src_ni+1)/2;
495  for (unsigned i=0;i<ni;++i)
496  {
497  dest_im[i]= *s0;
498  s0+=sxs2;
499  }
500 
501  if (src_nj&1)
502  {
503  unsigned yhi = (src_nj-1)/2;
504  vxl_byte* dest_last_row = dest_im + yhi*d_y_step;
505  const vxl_byte* s_last = src_im + yhi*sys2;
506  for (unsigned i=0;i<ni;++i)
507  {
508  dest_last_row[i]= *s_last;
509  s_last+=sxs2;
510  }
511  }
512 }
513 
514 
515 //: Smooth and subsample single plane src_im in x to produce dest_im using 121 filter in x and y
516 // Smooths with a 3x3 filter and subsamples
517 template <>
518 void vil_gauss_reduce_121_1plane(const float* src_im,
519  unsigned src_ni, unsigned src_nj,
520  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
521  float* dest_im,
522  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
523 {
524  std::ptrdiff_t sxs2 = s_x_step*2;
525  std::ptrdiff_t sys2 = s_y_step*2;
526  float* d_row = dest_im+d_y_step;
527  const float* s_row1 = src_im + s_y_step;
528  const float* s_row2 = s_row1 + s_y_step;
529  const float* s_row3 = s_row2 + s_y_step;
530  unsigned ni2 = (src_ni-2)/2;
531  unsigned nj2 = (src_nj-2)/2;
532  for (unsigned y=0;y<nj2;++y)
533  {
534  // Set first element of row
535  *d_row = *s_row2;
536  float * d = d_row + d_x_step;
537  const float* s1 = s_row1 + sxs2;
538  const float* s2 = s_row2 + sxs2;
539  const float* s3 = s_row3 + sxs2;
540  for (unsigned x=0;x<ni2;++x)
541  {
542  // The following is a little inefficient - could group terms to reduce arithmetic
543  *d = 0.0625f * s1[-s_x_step] + 0.125f * s1[0] + 0.0625f * s1[s_x_step]
544  + 0.1250f * s2[-s_x_step] + 0.250f * s2[0] + 0.1250f * s2[s_x_step]
545  + 0.0625f * s3[-s_x_step] + 0.125f * s3[0] + 0.0625f * s3[s_x_step];
546 
547  d += d_x_step;
548  s1 += sxs2;
549  s2 += sxs2;
550  s3 += sxs2;
551  }
552  // Set last elements of row
553  if (src_ni&1)
554  *d = *s2;
555 
556  d_row += d_y_step;
557  s_row1 += sys2;
558  s_row2 += sys2;
559  s_row3 += sys2;
560  }
561 
562  // Need to set first and last rows as well
563 
564  // Dest image should be (src_ni+1)/2 x (src_nj+1)/2
565  const float* s0 = src_im;
566  unsigned ni=(src_ni+1)/2;
567  for (unsigned i=0;i<ni;++i)
568  {
569  dest_im[i]= *s0;
570  s0+=sxs2;
571  }
572 
573  if (src_nj&1)
574  {
575  unsigned yhi = (src_nj-1)/2;
576  float* dest_last_row = dest_im + yhi*d_y_step;
577  const float* s_last = src_im + yhi*sys2;
578  for (unsigned i=0;i<ni;++i)
579  {
580  dest_last_row[i]= *s_last;
581  s_last+=sxs2;
582  }
583  }
584 }
585 
586 
587 //: Smooth and subsample single plane src_im in x to produce dest_im using 121 filter in x and y
588 // Smooths with a 3x3 filter and subsamples
589 template <>
590 void vil_gauss_reduce_121_1plane(const double* src_im,
591  unsigned src_ni, unsigned src_nj,
592  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
593  double* dest_im,
594  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
595 {
596  std::ptrdiff_t sxs2 = s_x_step*2;
597  std::ptrdiff_t sys2 = s_y_step*2;
598  double* d_row = dest_im+d_y_step;
599  const double* s_row1 = src_im + s_y_step;
600  const double* s_row2 = s_row1 + s_y_step;
601  const double* s_row3 = s_row2 + s_y_step;
602  unsigned ni2 = (src_ni-2)/2;
603  unsigned nj2 = (src_nj-2)/2;
604  for (unsigned y=0;y<nj2;++y)
605  {
606  // Set first element of row
607  *d_row = *s_row2;
608  double * d = d_row + d_x_step;
609  const double* s1 = s_row1 + sxs2;
610  const double* s2 = s_row2 + sxs2;
611  const double* s3 = s_row3 + sxs2;
612  for (unsigned x=0;x<ni2;++x)
613  {
614  // The following is a little inefficient - could group terms to reduce arithmetic
615  *d = 0.0625 * s1[-s_x_step] + 0.125 * s1[0] + 0.0625 * s1[s_x_step]
616  + 0.1250 * s2[-s_x_step] + 0.250 * s2[0] + 0.1250 * s2[s_x_step]
617  + 0.0625 * s3[-s_x_step] + 0.125 * s3[0] + 0.0625 * s3[s_x_step];
618 
619  d += d_x_step;
620  s1 += sxs2;
621  s2 += sxs2;
622  s3 += sxs2;
623  }
624  // Set last elements of row
625  if (src_ni&1)
626  *d = *s2;
627 
628  d_row += d_y_step;
629  s_row1 += sys2;
630  s_row2 += sys2;
631  s_row3 += sys2;
632  }
633 
634  // Need to set first and last rows as well
635 
636  // Dest image should be (src_ni+1)/2 x (src_nj+1)/2
637  const double* s0 = src_im;
638  unsigned ni=(src_ni+1)/2;
639  for (unsigned i=0;i<ni;++i)
640  {
641  dest_im[i]= *s0;
642  s0+=sxs2;
643  }
644 
645  if (src_nj&1)
646  {
647  unsigned yhi = (src_nj-1)/2;
648  double* dest_last_row = dest_im + yhi*d_y_step;
649  const double* s_last = src_im + yhi*sys2;
650  for (unsigned i=0;i<ni;++i)
651  {
652  dest_last_row[i]= *s_last;
653  s_last+=sxs2;
654  }
655  }
656 }
657 
658 
659 //: Smooth and subsample single plane src_im in x to produce dest_im using 121 filter in x and y
660 // Smooths with a 3x3 filter and subsamples
661 template <>
662 void vil_gauss_reduce_121_1plane(const int* src_im,
663  unsigned src_ni, unsigned src_nj,
664  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
665  int* dest_im,
666  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
667 {
668  std::ptrdiff_t sxs2 = s_x_step*2;
669  std::ptrdiff_t sys2 = s_y_step*2;
670  int* d_row = dest_im+d_y_step;
671  const int* s_row1 = src_im + s_y_step;
672  const int* s_row2 = s_row1 + s_y_step;
673  const int* s_row3 = s_row2 + s_y_step;
674  unsigned ni2 = (src_ni-2)/2;
675  unsigned nj2 = (src_nj-2)/2;
676  for (unsigned y=0;y<nj2;++y)
677  {
678  // Set first element of row
679  *d_row = *s_row2;
680  int * d = d_row + d_x_step;
681  const int* s1 = s_row1 + sxs2;
682  const int* s2 = s_row2 + sxs2;
683  const int* s3 = s_row3 + sxs2;
684  for (unsigned x=0;x<ni2;++x)
685  {
686  // The following is a little inefficient - could group terms to reduce arithmetic
687  // Add 0.5 so that truncating effectively rounds
688  *d = int( 0.0625f * s1[-s_x_step] + 0.125f * s1[0] + 0.0625f * s1[s_x_step]
689  + 0.1250f * s2[-s_x_step] + 0.250f * s2[0] + 0.1250f * s2[s_x_step]
690  + 0.0625f * s3[-s_x_step] + 0.125f * s3[0] + 0.0625f * s3[s_x_step] +0.5);
691 
692  d += d_x_step;
693  s1 += sxs2;
694  s2 += sxs2;
695  s3 += sxs2;
696  }
697  // Set last elements of row
698  if (src_ni&1)
699  *d = *s2;
700 
701  d_row += d_y_step;
702  s_row1 += sys2;
703  s_row2 += sys2;
704  s_row3 += sys2;
705  }
706 
707  // Need to set first and last rows as well
708 
709  // Dest image should be (src_ni+1)/2 x (src_nj+1)/2
710  const int* s0 = src_im;
711  unsigned ni=(src_ni+1)/2;
712  for (unsigned i=0;i<ni;++i)
713  {
714  dest_im[i]= *s0;
715  s0+=sxs2;
716  }
717 
718  if (src_nj&1)
719  {
720  unsigned yhi = (src_nj-1)/2;
721  int* dest_last_row = dest_im + yhi*d_y_step;
722  const int* s_last = src_im + yhi*sys2;
723  for (unsigned i=0;i<ni;++i)
724  {
725  dest_last_row[i]= *s_last;
726  s_last+=sxs2;
727  }
728  }
729 }
730 
731 //: Smooth and subsample single plane src_im in x to produce dest_im using 121 filter in x and y
732 // Smooths with a 3x3 filter and subsamples
733 template <>
734 void vil_gauss_reduce_121_1plane(const vxl_int_16* src_im,
735  unsigned src_ni, unsigned src_nj,
736  std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step,
737  vxl_int_16* dest_im,
738  std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
739 {
740  std::ptrdiff_t sxs2 = s_x_step*2;
741  std::ptrdiff_t sys2 = s_y_step*2;
742  vxl_int_16* d_row = dest_im+d_y_step;
743  const vxl_int_16* s_row1 = src_im + s_y_step;
744  const vxl_int_16* s_row2 = s_row1 + s_y_step;
745  const vxl_int_16* s_row3 = s_row2 + s_y_step;
746  unsigned ni2 = (src_ni-2)/2;
747  unsigned nj2 = (src_nj-2)/2;
748  for (unsigned y=0;y<nj2;++y)
749  {
750  // Set first element of row
751  *d_row = *s_row2;
752  vxl_int_16 * d = d_row + d_x_step;
753  const vxl_int_16* s1 = s_row1 + sxs2;
754  const vxl_int_16* s2 = s_row2 + sxs2;
755  const vxl_int_16* s3 = s_row3 + sxs2;
756  for (unsigned x=0;x<ni2;++x)
757  {
758  // The following is a little inefficient - could group terms to reduce arithmetic
759  // Add 0.5 so that truncating effectively rounds
760  *d = vxl_int_16( 0.0625f * s1[-s_x_step] + 0.125f * s1[0] + 0.0625f * s1[s_x_step]
761  + 0.1250f * s2[-s_x_step] + 0.250f * s2[0] + 0.1250f * s2[s_x_step]
762  + 0.0625f * s3[-s_x_step] + 0.125f * s3[0] + 0.0625f * s3[s_x_step] +0.5);
763 
764  d += d_x_step;
765  s1 += sxs2;
766  s2 += sxs2;
767  s3 += sxs2;
768  }
769  // Set last elements of row
770  if (src_ni&1)
771  *d = *s2;
772 
773  d_row += d_y_step;
774  s_row1 += sys2;
775  s_row2 += sys2;
776  s_row3 += sys2;
777  }
778 
779  // Need to set first and last rows as well
780 
781  // Dest image should be (src_ni+1)/2 x (src_nj+1)/2
782  const vxl_int_16* s0 = src_im;
783  unsigned ni=(src_ni+1)/2;
784  for (unsigned i=0;i<ni;++i)
785  {
786  dest_im[i]= *s0;
787  s0+=sxs2;
788  }
789 
790  if (src_nj&1)
791  {
792  unsigned yhi = (src_nj-1)/2;
793  vxl_int_16* dest_last_row = dest_im + yhi*d_y_step;
794  const vxl_int_16* s_last = src_im + yhi*sys2;
795  for (unsigned i=0;i<ni;++i)
796  {
797  dest_last_row[i]= *s_last;
798  s_last+=sxs2;
799  }
800  }
801 }
802 
803 
804 vil_gauss_reduce_params::vil_gauss_reduce_params(double scaleStep)
805 {
806  assert(scaleStep> 1.0 && scaleStep<=2.0);
807  scale_step_ = scaleStep;
808  // This arrangement gives close to a 1-5-8-5-1 filter for scalestep of 2.0;
809  // and 0-0-1-0-0 for a scale step close to 1.0;
810  double z = 1/std::sqrt(2.0*(scaleStep-1.0));
811  filt0_ = vnl_erf(0.5 * z) - vnl_erf(-0.5 * z);
812  filt1_ = vnl_erf(1.5 * z) - vnl_erf(0.5 * z);
813  filt2_ = vnl_erf(2.5 * z) - vnl_erf(1.5 * z);
814 
815  double five_tap_total = 2*(filt2_ + filt1_) + filt0_;
816 #if 0
817  double four_tap_total = filt2_ + 2*(filt1_) + filt0_;
818  double three_tap_total = filt2_ + filt1_ + filt0_;
819 #endif
820 
821  // Calculate 3 tap half Gaussian filter assuming constant edge extension
822  filt_edge0_ = (filt0_ + filt1_ + filt2_) / five_tap_total;
823  filt_edge1_ = filt1_ / five_tap_total;
824  filt_edge2_ = filt2_ / five_tap_total;
825 #if 0
826  filt_edge0_ = 1.0;
827  filt_edge1_ = 0.0;
828  filt_edge2_ = 0.0;
829 #endif
830  // Calculate 4 tap skewed Gaussian filter assuming constant edge extension
831  filt_pen_edge_n1_ = (filt1_+filt2_) / five_tap_total;
832  filt_pen_edge0_ = filt0_ / five_tap_total;
833  filt_pen_edge1_ = filt1_ / five_tap_total;
834  filt_pen_edge2_ = filt2_ / five_tap_total;
835 
836  // Calculate 5 tap Gaussian filter
837  filt0_ = filt0_ / five_tap_total;
838  filt1_ = filt1_ / five_tap_total;
839  filt2_ = filt2_ / five_tap_total;
840 
841  assert(filt_edge0_ > filt_edge1_);
842  assert(filt_edge1_ > filt_edge2_);
843 }
vnl_erf.h
vil_gauss_reduce_params::filt_edge1_
double filt_edge1_
Definition: vil_gauss_reduce.h:48
vil_gauss_reduce_params::filt_pen_edge2_
double filt_pen_edge2_
Definition: vil_gauss_reduce.h:49
vnl_math::rnd
int rnd(float x)
vnl_math.h
vil_gauss_reduce_params::filt_pen_edge1_
double filt_pen_edge1_
Definition: vil_gauss_reduce.h:49
vil_gauss_reduce_1plane
void vil_gauss_reduce_1plane(const vxl_byte *src_im, unsigned src_ni, unsigned src_nj, std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step, vxl_byte *dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
Smooth and subsample single plane src_im in x to produce dest_im.
Definition: vil_gauss_reduce.cxx:22
vil_gauss_reduce_params::filt1_
double filt1_
Definition: vil_gauss_reduce.h:47
vil_gauss_reduce_params::filt_pen_edge0_
double filt_pen_edge0_
Definition: vil_gauss_reduce.h:49
vil_gauss_reduce_params::filt0_
double filt0_
Definition: vil_gauss_reduce.h:47
vil_gauss_reduce_params::filt_edge0_
double filt_edge0_
Definition: vil_gauss_reduce.h:48
vil_gauss_reduce_params::filt_pen_edge_n1_
double filt_pen_edge_n1_
Definition: vil_gauss_reduce.h:49
vnl_erf
double vnl_erf(double x)
vil_gauss_reduce_params::filt_edge2_
double filt_edge2_
Definition: vil_gauss_reduce.h:48
vil_gauss_reduce_2_3_1plane
void vil_gauss_reduce_2_3_1plane(const float *src_im, unsigned src_ni, unsigned src_nj, std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step, float *dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
Smooth and subsample single plane src_im in x, result is 2/3rd size.
Definition: vil_gauss_reduce.cxx:237
vil_gauss_reduce_params::filt2_
double filt2_
Definition: vil_gauss_reduce.h:47
vil_gauss_reduce.h
Functions to smooth and sub-sample image in one direction.
vil_gauss_reduce_params::vil_gauss_reduce_params
vil_gauss_reduce_params(double scale_step)
Definition: vil_gauss_reduce.cxx:804
vil_gauss_reduce_121_1plane
void vil_gauss_reduce_121_1plane(const vxl_byte *src_im, unsigned src_ni, unsigned src_nj, std::ptrdiff_t s_x_step, std::ptrdiff_t s_y_step, vxl_byte *dest_im, std::ptrdiff_t d_x_step, std::ptrdiff_t d_y_step)
Smooth and subsample single plane src_im in x to produce dest_im using 121 filter in x and y.
Definition: vil_gauss_reduce.cxx:446
vil_gauss_reduce_params::scale_step_
double scale_step_
Definition: vil_gauss_reduce.h:46
Generated by   doxygen 1.8.15