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