1 // This is core/vil/algo/vil_convolve_1d.h
2 #ifndef vil_convolve_1d_h_
3 #define vil_convolve_1d_h_
4 //:
5 // \file
6 // \brief 1D Convolution with cunning boundary options
7 // \author Tim Cootes, Ian Scott (based on work by fsm)
8 //
9 // \note The convolution operation is defined by
10 // $(f*g)(x) = \int f(x-y) g(y) dy$
11 // i.e. the kernel g is reflected before the integration is performed.
12 // If you don't want this to happen, the behaviour you want is not
13 // called "convolution". So don't break the convolution routines in
14 // that particular way.
15
16 #include <algorithm>
17 #include <cstdlib>
18 #include <cstring>
19 #include <iostream>
20 #ifdef _MSC_VER
21 # include <vcl_msvc_warnings.h>
22 #endif
23 #include <cassert>
24 #include <vil/vil_image_view.h>
25 #include <vil/vil_image_resource.h>
26 #include <vil/vil_property.h>
27
28
29 //: Available options for boundary behavior
30 // When convolving a finite signal the boundaries may be
31 // treated in various ways which can often be expressed in terms
32 // of ways to extend the signal outside its original range.
33 enum vil_convolve_boundary_option
34 {
35 //: Do not fill destination edges at all.
36 // i.e. leave them unchanged.
37 vil_convolve_ignore_edge,
38
39 //: Do not to extend the signal, but pad with zeros.
40 // \verbatim
41 // | |
42 // K ----*-------
43 // in ... ---------------------------
44 // out ... --------------------0000000
45 // \endverbatim
46 vil_convolve_no_extend,
47
48 //: Zero-extend the input signal beyond the boundary.
49 // \verbatim
50 // | |
51 // K ----*--------
52 // in ... ---------------------------000000000000...
53 // out ... ---------------------------
54 // \endverbatim
55 vil_convolve_zero_extend,
56
57 //: Extend the signal to be constant beyond the boundary.
58 // \verbatim
59 // | |
60 // K ----*--------
61 // in ... --------------------------aaaaaaaaaaaaa...
62 // out ... ---------------------------
63 // \endverbatim
64 vil_convolve_constant_extend,
65
66 //: Extend the signal periodically beyond the boundary.
67 // \verbatim
68 // | |
69 // K ----*--------
70 // in abc...-------------------------abc...------..
71 // out ... ---------------------------
72 // \endverbatim
73 vil_convolve_periodic_extend,
74
75 //: Extend the signal by reflection about the boundary.
76 // \verbatim
77 // | |
78 // K ----*--------
79 // in ... -------------------...edcbabcde...
80 // out ... ---------------------------
81 // \endverbatim
82 vil_convolve_reflect_extend,
83
84 //: Kernel is trimmed and reweighed, to allow convolution up to boundary.
85 // This one is slightly different. The input signal is not
86 // extended in any way, but the kernel is trimmed to allow
87 // convolution to proceed up to the boundary and reweighed
88 // to keep the total area the same.
89 // \note may not work with kernels which take negative values.
90 vil_convolve_trim
91 };
92
93 //: Convolve edge with kernel[x*kstep] x in [k_lo,k_hi] (k_hi>=0)
94 // Fills only edge: dest[i], i=0..(k_hi-1)
95 template <class srcT, class destT, class kernelT, class accumT>
vil_convolve_edge_1d(const srcT * src,unsigned n,std::ptrdiff_t s_step,destT * dest,std::ptrdiff_t d_step,const kernelT * kernel,std::ptrdiff_t k_lo,std::ptrdiff_t k_hi,std::ptrdiff_t kstep,accumT,vil_convolve_boundary_option option)96 inline void vil_convolve_edge_1d(const srcT* src, unsigned n, std::ptrdiff_t s_step,
97 destT* dest, std::ptrdiff_t d_step,
98 const kernelT* kernel,
99 std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
100 std::ptrdiff_t kstep, accumT,
101 vil_convolve_boundary_option option)
102 {
103 switch (option)
104 {
105 case vil_convolve_ignore_edge:
106 return;
107 case vil_convolve_no_extend:
108 // Initialise first elements of row to zero
109 for (std::ptrdiff_t i=-k_hi;i<0;++i,dest+=d_step)
110 *dest = 0;
111 return;
112 case vil_convolve_zero_extend:
113 // Assume src[i]==0 for i<0
114 // for (std::ptrdiff_t i=-k_hi+1;i<=0;++i,dest+=d_step,src+=s_step)
115 for (std::ptrdiff_t i=0;i<k_hi;++i,dest+=d_step)
116 {
117 accumT sum = 0;
118 const srcT* s = src;
119 const kernelT* k = kernel+i*kstep;
120 for (std::ptrdiff_t j=i;j>=k_lo;--j,s+=s_step,k-=kstep)
121 sum+= (accumT)((*s)*(*k));
122 *dest=(destT)sum;
123 }
124 return;
125 case vil_convolve_constant_extend:
126 {
127 // Assume src[i]=src[0] for i<0
128 std::ptrdiff_t i_max = k_hi-1;
129 for (std::ptrdiff_t i=0;i<=i_max;++i)
130 {
131 accumT sum=0;
132 for (std::ptrdiff_t j=-k_hi;j<=-k_lo;++j)
133 {
134 if ((i+j)<0) sum+=(accumT)(src[0]*kernel[j*(-kstep)]);
135 else sum+=(accumT)(src[(i+j)*s_step]*kernel[j*(-kstep)]);
136 }
137 dest[i*d_step]=(destT)sum;
138 }
139 return;
140 }
141 case vil_convolve_reflect_extend:
142 {
143 // Assume src[i]=src[0] for i<0
144 std::ptrdiff_t i_max = k_hi-1;
145 for (std::ptrdiff_t i=0;i<=i_max;++i)
146 {
147 accumT sum=0;
148 for (std::ptrdiff_t j=-k_hi;j<=-k_lo;++j)
149 {
150 if ((i+j)<0) sum+=(accumT)(src[-(i+j)*s_step]*kernel[j*(-kstep)]);
151 else sum+=(accumT)(src[(i+j)*s_step]*kernel[j*(-kstep)]);
152 }
153 dest[i*d_step]=(destT)sum;
154 }
155 return;
156 }
157 case vil_convolve_periodic_extend:
158 {
159 // Assume src[i]=src[n+i] for i<0
160 std::ptrdiff_t i_max = k_hi-1;
161 for (int i=0;i<=i_max;++i)
162 {
163 accumT sum=0;
164 for (std::ptrdiff_t j=k_hi;j>=k_lo;--j)
165 sum+=(accumT)(src[((i-j+n)%n)*s_step]*kernel[j*kstep]);
166 dest[i*d_step]=(destT)sum;
167 }
168 return;
169 }
170 case vil_convolve_trim:
171 {
172 // Truncate and reweight kernel
173 accumT k_sum_all=0;
174 for (std::ptrdiff_t j=-k_hi;j<=-k_lo;++j) k_sum_all+=(accumT)(kernel[j*(-kstep)]);
175
176 std::ptrdiff_t i_max = k_hi-1;
177 for (std::ptrdiff_t i=0;i<=i_max;++i)
178 {
179 accumT sum=0;
180 accumT k_sum=0;
181 // Sum elements which overlap src
182 // ie i+j>=0 (so j starts at -i)
183 for (std::ptrdiff_t j=-i;j<=-k_lo;++j)
184 {
185 sum+=(accumT)(src[(i+j)*s_step]*kernel[j*(-kstep)]);
186 k_sum += (accumT)(kernel[j*(-kstep)]);
187 }
188 dest[i*d_step]=(destT)(sum*k_sum_all/k_sum);
189 }
190 return;
191 }
192 default:
193 std::cout<<"ERROR: vil_convolve_edge_1d: "
194 <<"Sorry, can't deal with supplied edge option.\n";
195 std::abort();
196 }
197 }
198
199 //: Convolve kernel[x] (x in [k_lo,k_hi]) with srcT
200 // Assumes dest and src same size (nx)
201 // Kernel must not be larger than nx;
202 template <class srcT, class destT, class kernelT, class accumT>
vil_convolve_1d(const srcT * src0,unsigned nx,std::ptrdiff_t s_step,destT * dest0,std::ptrdiff_t d_step,const kernelT * kernel,std::ptrdiff_t k_lo,std::ptrdiff_t k_hi,accumT ac,vil_convolve_boundary_option start_option,vil_convolve_boundary_option end_option)203 inline void vil_convolve_1d(const srcT* src0, unsigned nx, std::ptrdiff_t s_step,
204 destT* dest0, std::ptrdiff_t d_step,
205 const kernelT* kernel,
206 std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
207 accumT ac,
208 vil_convolve_boundary_option start_option,
209 vil_convolve_boundary_option end_option)
210 {
211 assert(k_hi - k_lo < int(nx));
212
213 // Deal with start (fill elements 0..1+k_hi of dest)
214 vil_convolve_edge_1d(src0,nx,s_step,dest0,d_step,kernel,k_lo,k_hi,1,ac,start_option);
215
216 const kernelT* k_rbegin = kernel+k_hi;
217 const kernelT* k_rend = kernel+k_lo-1;
218 assert(k_rbegin >= k_rend);
219 const srcT* src = src0;
220
221 for (destT * dest = dest0 + d_step*k_hi,
222 * const end_dest = dest0 + d_step*(int(nx)+k_lo);
223 dest!=end_dest;
224 dest+=d_step,src+=s_step)
225 {
226 accumT sum = 0;
227 const srcT* s= src;
228 for (const kernelT *k = k_rbegin;k!=k_rend;--k,s+=s_step)
229 sum+= (accumT)((*k)*(*s));
230 *dest = destT(sum);
231 }
232
233 // Deal with end (reflect data and kernel!)
234 vil_convolve_edge_1d(src0+(nx-1)*s_step,nx,-s_step,
235 dest0+(nx-1)*d_step,-d_step,
236 kernel,-k_hi,-k_lo,-1,ac,end_option);
237 }
238
239 //: Convolve kernel[i] (i in [k_lo,k_hi]) with srcT in i-direction
240 // On exit dest_im(i,j) = sum src(i-x,j)*kernel(x) (x=k_lo..k_hi)
241 // \note This function reverses the kernel. If you don't want the
242 // kernel reversed, use vil_correlate_1d instead. The kernel must
243 // not be larger than src_im.ni()
244 // \param kernel should point to tap 0.
245 // \param dest_im will be resized to size of src_im.
246 // \relatesalso vil_image_view
247 template <class srcT, class destT, class kernelT, class accumT>
vil_convolve_1d(const vil_image_view<srcT> & src_im,vil_image_view<destT> & dest_im,const kernelT * kernel,std::ptrdiff_t k_lo,std::ptrdiff_t k_hi,accumT ac,vil_convolve_boundary_option start_option,vil_convolve_boundary_option end_option)248 inline void vil_convolve_1d(const vil_image_view<srcT>& src_im,
249 vil_image_view<destT>& dest_im,
250 const kernelT* kernel,
251 std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
252 accumT ac,
253 vil_convolve_boundary_option start_option,
254 vil_convolve_boundary_option end_option)
255 {
256 unsigned n_i = src_im.ni();
257 unsigned n_j = src_im.nj();
258 assert(k_hi - k_lo +1 <= (int) n_i);
259 std::ptrdiff_t s_istep = src_im.istep(), s_jstep = src_im.jstep();
260
261 dest_im.set_size(n_i,n_j,src_im.nplanes());
262 std::ptrdiff_t d_istep = dest_im.istep(),d_jstep = dest_im.jstep();
263
264 for (unsigned int p=0;p<src_im.nplanes();++p)
265 {
266 // Select first row of p-th plane
267 const srcT* src_row = src_im.top_left_ptr()+p*src_im.planestep();
268 destT* dest_row = dest_im.top_left_ptr()+p*dest_im.planestep();
269
270 // Apply convolution to each row in turn
271 // First check if either istep is 1 for speed optimisation.
272
273 if (s_istep == 1)
274 {
275 if (d_istep == 1)
276 for (unsigned int j=0;j<n_j;++j,src_row+=s_jstep,dest_row+=d_jstep)
277 vil_convolve_1d(src_row,n_i,1, dest_row,1,
278 kernel,k_lo,k_hi,ac,start_option,end_option);
279 else
280 for (unsigned int j=0;j<n_j;++j,src_row+=s_jstep,dest_row+=d_jstep)
281 vil_convolve_1d(src_row,n_i,1, dest_row,d_istep,
282 kernel,k_lo,k_hi,ac,start_option,end_option);
283 }
284 else
285 {
286 if (d_istep == 1)
287 for (unsigned int j=0;j<n_j;++j,src_row+=s_jstep,dest_row+=d_jstep)
288 vil_convolve_1d(src_row,n_i,s_istep, dest_row,1,
289 kernel,k_lo,k_hi,ac,start_option,end_option);
290 else
291 for (unsigned int j=0;j<n_j;++j,src_row+=s_jstep,dest_row+=d_jstep)
292 vil_convolve_1d(src_row,n_i,s_istep, dest_row,d_istep,
293 kernel,k_lo,k_hi,ac,start_option,end_option);
294 }
295 }
296 }
297
298 template <class destT, class kernelT, class accumT>
299 vil_image_resource_sptr vil_convolve_1d(
300 const vil_image_resource_sptr& src_im,
301 const destT dt,
302 const kernelT* kernel, int k_lo, int k_hi,
303 const accumT ac,
304 vil_convolve_boundary_option start_option,
305 vil_convolve_boundary_option end_option);
306
307 //: A resource adaptor that behaves like a convolved version of its input
308 template <class kernelT, class accumT, class destT>
309 class vil_convolve_1d_resource : public vil_image_resource
310 {
311 public:
312 //: Construct a convolve filter.
313 // You can't create one of these directly, use vil_convolve_1d instead
vil_convolve_1d_resource(const vil_image_resource_sptr & src,const kernelT * kernel,int k_lo,int k_hi,vil_convolve_boundary_option start_option,vil_convolve_boundary_option end_option)314 vil_convolve_1d_resource(const vil_image_resource_sptr& src,
315 const kernelT* kernel, int k_lo, int k_hi,
316 vil_convolve_boundary_option start_option,
317 vil_convolve_boundary_option end_option) :
318 src_(src), kernel_(kernel), klo_(k_lo), khi_(k_hi),
319 start_option_(start_option), end_option_(end_option)
320 {
321 // Can't do period extension yet.
322 assert (start_option != vil_convolve_periodic_extend ||
323 end_option != vil_convolve_periodic_extend);
324 }
325
326 friend vil_image_resource_sptr vil_convolve_1d <> (
327 const vil_image_resource_sptr& src_im, const destT dt, const kernelT* kernel,
328 int k_lo, int k_hi, const accumT ac,
329 vil_convolve_boundary_option start_option,
330 vil_convolve_boundary_option end_option);
331
332 public:
get_copy_view(unsigned i0,unsigned n_i,unsigned j0,unsigned n_j)333 vil_image_view_base_sptr get_copy_view(unsigned i0, unsigned n_i,
334 unsigned j0, unsigned n_j) const override
335 {
336 if (i0 + n_i > src_->ni() || j0 + n_j > src_->nj()) return nullptr;
337 const unsigned lsrc = (unsigned) std::max(0,(int)i0 + klo_); // lhs of input window
338 const unsigned hsrc = std::min(src_->ni(),i0 + n_i - klo_ + khi_); // 1+rhs of input window.
339 const unsigned lboundary = std::min((unsigned) -klo_, i0); // width of lhs boundary area.
340 assert (hsrc > lsrc);
341 vil_image_view_base_sptr vs = src_->get_view(lsrc, hsrc-lsrc, j0, n_j);
342 vil_image_view<destT> dest(vs->ni(), vs->nj(), vs->nplanes());
343 switch (vs->pixel_format())
344 {
345 #define macro( F , T ) \
346 case F : \
347 vil_convolve_1d(static_cast<vil_image_view<T >&>(*vs),dest, \
348 kernel_, klo_, khi_, accumT(), start_option_, end_option_); \
349 return new vil_image_view<destT>(vil_crop(dest, lboundary, n_i, 0, n_j));
350
351 macro(VIL_PIXEL_FORMAT_BYTE , vxl_byte )
352 macro(VIL_PIXEL_FORMAT_SBYTE , vxl_sbyte )
353 macro(VIL_PIXEL_FORMAT_UINT_32 , vxl_uint_32 )
354 macro(VIL_PIXEL_FORMAT_UINT_16 , vxl_uint_16 )
355 macro(VIL_PIXEL_FORMAT_INT_32 , vxl_int_32 )
356 macro(VIL_PIXEL_FORMAT_INT_16 , vxl_int_16 )
357 macro(VIL_PIXEL_FORMAT_BOOL , bool )
358 macro(VIL_PIXEL_FORMAT_FLOAT , float )
359 macro(VIL_PIXEL_FORMAT_DOUBLE , double )
360 // complex<float> should work - but causes all manner of compiler template errors.
361 // maybe need a better compiler, maybe there is a code fix - IMS
362 #undef macro
363 default:
364 return nullptr;
365 }
366 }
367
nplanes()368 unsigned nplanes() const override { return src_->nplanes(); }
ni()369 unsigned ni() const override { return src_->ni(); }
nj()370 unsigned nj() const override { return src_->nj(); }
371
pixel_format()372 enum vil_pixel_format pixel_format() const override
373 { return vil_pixel_format_of(accumT()); }
374
375
376 //: Put the data in this view back into the image source.
put_view(const vil_image_view_base &,unsigned,unsigned)377 bool put_view(const vil_image_view_base& /*im*/, unsigned /*i0*/, unsigned /*j0*/) override
378 {
379 std::cerr << "WARNING: vil_convolve_1d_resource::put_back\n"
380 << "\tYou can't push data back into a convolve filter.\n";
381 return false;
382 }
383
384 //: Extra property information
385 bool get_property(char const* tag, void* property_value = nullptr) const override
386 {
387 if (0==std::strcmp(tag, vil_property_read_only))
388 return property_value ? (*static_cast<bool*>(property_value)) = true : true;
389
390 return src_->get_property(tag, property_value);
391 }
392
393 protected:
394 vil_image_resource_sptr src_;
395 const kernelT* kernel_;
396 int klo_, khi_;
397 vil_convolve_boundary_option start_option_, end_option_;
398 };
399
400 //: Create an image_resource object which convolve kernel[x] x in [k_lo,k_hi] with srcT
401 // \note This function reverses the kernel. If you don't want the
402 // kernel reversed, use vil_correlate_1d instead.
403 // \param kernel should point to tap 0.
404 // \relatesalso vil_image_resource
405 template <class destT, class kernelT, class accumT>
vil_convolve_1d(const vil_image_resource_sptr & src_im,const destT,const kernelT * kernel,int k_lo,int k_hi,const accumT,vil_convolve_boundary_option start_option,vil_convolve_boundary_option end_option)406 vil_image_resource_sptr vil_convolve_1d(
407 const vil_image_resource_sptr& src_im,
408 const destT /*dt*/,
409 const kernelT* kernel, int k_lo, int k_hi,
410 const accumT,
411 vil_convolve_boundary_option start_option,
412 vil_convolve_boundary_option end_option)
413 {
414 return new vil_convolve_1d_resource<kernelT, accumT, destT>(src_im, kernel, k_lo, k_hi, start_option, end_option);
415 }
416
417 #endif // vil_convolve_1d_h_
418