1 // This is core/vil1/vil1_convolve_1d_y_hxx_
2 #ifndef vil1_convolve_1d_y_hxx_
3 #define vil1_convolve_1d_y_hxx_
4
5 #ifndef fsm_dont_croak // this file should only be #included from vil1_convolve_1d.hxx
6 croak
7 #endif
8
9 #include <cstdlib>
10 #include "vil1_convolve.h"
11 #ifdef _MSC_VER
12 # include <vcl_msvc_warnings.h>
13 #endif
14
15 template <class I1, class I2, class AC, class O>
vil1_convolve_1d_y(vil1_convolve_signal_1d<I1 const> const & kernel,vil1_convolve_signal_2d<I2 const> const & input,AC *,vil1_convolve_signal_2d<O> const & output,vil1_convolve_boundary_option b,vil1_convolve_boundary_option e)16 void vil1_convolve_1d_y(vil1_convolve_signal_1d<I1 const> const &kernel,
17 vil1_convolve_signal_2d<I2 const> const &input,
18 AC *,
19 vil1_convolve_signal_2d<O> const &output,
20 vil1_convolve_boundary_option b,
21 vil1_convolve_boundary_option e)
22 {
23 // compute ranges of i, x, y here.
24 int i0 = kernel.begin_-kernel.origin_;
25 int i1 = kernel.end_ -kernel.origin_;
26
27 int x0 = output.beginx_-output.originx_;
28 int x1 = output.endx_ -output.originx_;
29
30 int y0 = output.beginy_-output.originy_;
31 int y1 = output.endy_ -output.originy_;
32
33 // compute total weight of the kernel.
34 // FIXME assumes non-negative kernel.
35 AC total_weight = 0;
36 for (int i=i0; i<i1; ++i)
37 total_weight += AC(value1d(kernel, i));
38
39
40 // this is not very efficient at the moment, but my main
41 // concern for now is that it works correctly.
42 for (int y=y0; y<y1; ++y) {
43 for (int x=x0; x<x1; ++x) {
44 AC ac = 0; // accumulated "kernel * input" terms.
45 AC wt = 0; // accumulated "kernel" terms.
46 bool zero = false;
47
48 for (int i=i0; i<i1 && !zero; ++i) {
49 // value of kernel at i :
50 AC kval = AC(value1d(kernel, i));
51
52 int yy = y-i;
53 if (yy < y0) switch (b) {
54 case vil1_convolve_no_extend:
55 zero = true; /*FIXME*/
56 break;
57 case vil1_convolve_zero_extend:
58 wt += kval;
59 break;
60 case vil1_convolve_constant_extend:
61 ac += kval * AC(value2d(input, x, y0));
62 wt += kval;
63 break;
64 case vil1_convolve_periodic_extend:
65 ac += kval * AC(value2d(input, x, yy+(y1-y0)));
66 wt += kval;
67 break;
68 case vil1_convolve_reflect_extend:
69 ac += kval * AC(value2d(input, x, 2*y0-yy));
70 wt += kval;
71 break;
72 case vil1_convolve_trim:
73 break;
74 default:
75 std::abort();
76 break;
77 }
78
79 else if (yy >= y1) switch (e) {
80 case vil1_convolve_no_extend:
81 zero = true; /*FIXME*/
82 break;
83 case vil1_convolve_zero_extend:
84 wt += kval;
85 break;
86 case vil1_convolve_constant_extend:
87 ac += kval * AC(value2d(input, x, y1-1));
88 wt += kval;
89 break;
90 case vil1_convolve_periodic_extend:
91 ac += kval * AC(value2d(input, x, yy-(y1-y0)));
92 wt += kval;
93 break;
94 case vil1_convolve_reflect_extend:
95 ac += kval * AC(value2d(input, x, 2*(y1-1)-yy));
96 wt += kval;
97 break;
98 case vil1_convolve_trim:
99 break;
100 default:
101 std::abort();
102 break;
103 }
104
105 else {
106 ac += kval * AC(value2d(input, x, yy));
107 wt += kval;
108 }
109 }
110
111 // compute and store final value.
112 if (zero)
113 value2d(output, x, y) = AC(0);
114 else if (wt)
115 value2d(output, x, y) = O(ac * total_weight / wt);
116 }
117 }
118 }
119
120 #endif // vil1_convolve_1d_y_hxx_
121