1 // This is core/vil/algo/vil_correlate_1d.h
2 #ifndef vil_correlate_1d_h_
3 #define vil_correlate_1d_h_
4 //:
5 // \file
6 // \brief 1D Convolution with cunning boundary options
7 // \author Tim Cootes (based on work by fsm)
8
9 #include <algorithm>
10 #include <cstring>
11 #include <iostream>
12 #ifdef _MSC_VER
13 # include <vcl_msvc_warnings.h>
14 #endif
15 #include <cassert>
16 #include <vil/vil_image_view.h>
17 #include <vil/vil_image_resource.h>
18 #include <vil/vil_property.h>
19 #include <vil/algo/vil_convolve_1d.h>
20
21 //: Correlate kernel[x] (x in [k_lo,k_hi]) with srcT
22 // Assumes dest and src same size (nx)
23 template <class srcT, class destT, class kernelT, class accumT>
vil_correlate_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)24 inline void vil_correlate_1d(const srcT* src0, unsigned nx, std::ptrdiff_t s_step,
25 destT* dest0, std::ptrdiff_t d_step,
26 const kernelT* kernel,
27 std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
28 accumT ac,
29 vil_convolve_boundary_option start_option,
30 vil_convolve_boundary_option end_option)
31 {
32 // Deal with start (fill elements 0..1-k_lo of dest)
33 vil_convolve_edge_1d(src0,nx,s_step,dest0,d_step,kernel,-k_hi,-k_lo,-1,ac,start_option);
34
35 const kernelT* k_begin = kernel+k_lo;
36 const kernelT* k_end = kernel+k_hi+1;
37 const srcT* src = src0;
38
39 destT* end_dest = dest0 + d_step*(int(nx)-k_hi);
40 for (destT* dest = dest0-d_step*k_lo;dest!=end_dest;dest+=d_step,src+=s_step)
41 {
42 accumT sum = 0;
43 const srcT* s= src;
44 for (const kernelT *k = k_begin;k!=k_end; ++k,s+=s_step) sum+= (accumT)((*k)*(*s));
45 *dest = destT(sum);
46 }
47
48 // Deal with end (reflect data and kernel!)
49 vil_convolve_edge_1d(src0+(nx-1)*s_step,nx,-s_step,
50 dest0+(nx-1)*d_step,-d_step,
51 kernel,k_lo,k_hi,1,ac,end_option);
52 }
53
54 //: correlate kernel[i] (i in [k_lo,k_hi]) with srcT in i-direction
55 // On exit dest_im(i,j) = sum src(i+x,j)*kernel(x) (x=k_lo..k_hi)
56 // \note This function does not reverse the kernel. If you want the
57 // kernel reversed, use vil_convolve_1d instead.
58 // \param kernel should point to tap 0.
59 // \param dest_im will be resized to size of src_im.
60 // \relatesalso vil_image_view
61 template <class srcT, class destT, class kernelT, class accumT>
vil_correlate_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)62 inline void vil_correlate_1d(const vil_image_view<srcT>& src_im,
63 vil_image_view<destT>& dest_im,
64 const kernelT* kernel,
65 std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
66 accumT ac,
67 vil_convolve_boundary_option start_option,
68 vil_convolve_boundary_option end_option)
69 {
70 unsigned ni = src_im.ni();
71 unsigned nj = src_im.nj();
72 std::ptrdiff_t s_istep = src_im.istep(), s_jstep = src_im.jstep();
73
74 dest_im.set_size(ni,nj,src_im.nplanes());
75 std::ptrdiff_t d_istep = dest_im.istep(),d_jstep = dest_im.jstep();
76
77 for (unsigned int p=0;p<src_im.nplanes();++p)
78 {
79 // Select first row of p-th plane
80 const srcT* src_row = src_im.top_left_ptr()+p*src_im.planestep();
81 destT* dest_row = dest_im.top_left_ptr()+p*dest_im.planestep();
82
83 // Apply convolution to each row in turn
84 // First check if either istep is 1 for speed optimisation.
85 if (s_istep == 1)
86 {
87 if (d_istep == 1)
88 for (unsigned int j=0;j<nj;++j,src_row+=s_jstep,dest_row+=d_jstep)
89 vil_correlate_1d(src_row,ni,1, dest_row,1,
90 kernel,k_lo,k_hi,ac,start_option,end_option);
91 else
92 for (unsigned int j=0;j<nj;++j,src_row+=s_jstep,dest_row+=d_jstep)
93 vil_correlate_1d(src_row,ni,1, dest_row,d_istep,
94 kernel,k_lo,k_hi,ac,start_option,end_option);
95 }
96 else
97 {
98 if (d_istep == 1)
99 for (unsigned int j=0;j<nj;++j,src_row+=s_jstep,dest_row+=d_jstep)
100 vil_correlate_1d(src_row,ni,s_istep, dest_row,1,
101 kernel,k_lo,k_hi,ac,start_option,end_option);
102 else
103 for (unsigned int j=0;j<nj;++j,src_row+=s_jstep,dest_row+=d_jstep)
104 vil_correlate_1d(src_row,ni,s_istep, dest_row,d_istep,
105 kernel,k_lo,k_hi,ac,start_option,end_option);
106 }
107 }
108 }
109
110 template <class destT, class kernelT, class accumT>
111 vil_image_resource_sptr vil_correlate_1d(
112 const vil_image_resource_sptr& src_im,
113 const destT dt,
114 const kernelT* kernel, std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
115 const accumT ac,
116 vil_convolve_boundary_option start_option,
117 vil_convolve_boundary_option end_option);
118
119 //: A resource adaptor that behaves like a correlated version of its input
120 template <class kernelT, class accumT, class destT>
121 class vil_correlate_1d_resource : public vil_image_resource
122 {
123 public:
124 //: Construct a correlate filter.
125 // You can't create one of these directly, use vil_correlate_1d instead
vil_correlate_1d_resource(const vil_image_resource_sptr & src,const kernelT * kernel,std::ptrdiff_t k_lo,std::ptrdiff_t k_hi,vil_convolve_boundary_option start_option,vil_convolve_boundary_option end_option)126 vil_correlate_1d_resource(const vil_image_resource_sptr& src,
127 const kernelT* kernel, std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
128 vil_convolve_boundary_option start_option,
129 vil_convolve_boundary_option end_option) :
130 src_(src), kernel_(kernel), klo_(k_lo), khi_(k_hi),
131 start_option_(start_option), end_option_(end_option)
132 {
133 // Can't do period extension yet.
134 assert (start_option != vil_convolve_periodic_extend ||
135 end_option != vil_convolve_periodic_extend);
136 }
137
138 friend vil_image_resource_sptr vil_correlate_1d <> (
139 const vil_image_resource_sptr& src_im, const destT dt, const kernelT* kernel,
140 std::ptrdiff_t k_lo, std::ptrdiff_t k_hi, const accumT ac,
141 vil_convolve_boundary_option start_option,
142 vil_convolve_boundary_option end_option);
143
144 public:
get_copy_view(unsigned i0,unsigned ni,unsigned j0,unsigned nj)145 vil_image_view_base_sptr get_copy_view(unsigned i0, unsigned ni,
146 unsigned j0, unsigned nj) const override
147 {
148 if (i0 + ni > src_->ni() || j0 + nj > src_->nj()) return nullptr;
149 const unsigned lsrc = (unsigned)std::max(0,int(i0+klo_)); // lhs of input window
150 const unsigned hsrc = std::min(src_->ni(),(unsigned int)(i0+ni-klo_+khi_)); // 1+rhs of input window.
151 const unsigned lboundary = std::min((unsigned) -klo_, i0); // width of lhs boundary area.
152 assert (hsrc > lsrc);
153 vil_image_view_base_sptr vs = src_->get_view(lsrc, hsrc-lsrc, j0, nj);
154 vil_image_view<destT> dest(vs->ni(), vs->nj(), vs->nplanes());
155 switch (vs->pixel_format())
156 {
157 #define macro( F , T ) \
158 case F : \
159 vil_correlate_1d(static_cast<vil_image_view<T >&>(*vs),dest, \
160 kernel_, klo_, khi_, accumT(), start_option_, end_option_); \
161 return new vil_image_view<destT>(vil_crop(dest, lboundary, ni, 0, nj));
162
163 macro(VIL_PIXEL_FORMAT_BYTE , vxl_byte )
164 macro(VIL_PIXEL_FORMAT_SBYTE , vxl_sbyte )
165 macro(VIL_PIXEL_FORMAT_UINT_32 , vxl_uint_32 )
166 macro(VIL_PIXEL_FORMAT_UINT_16 , vxl_uint_16 )
167 macro(VIL_PIXEL_FORMAT_INT_32 , vxl_int_32 )
168 macro(VIL_PIXEL_FORMAT_INT_16 , vxl_int_16 )
169 macro(VIL_PIXEL_FORMAT_BOOL , bool )
170 macro(VIL_PIXEL_FORMAT_FLOAT , float )
171 macro(VIL_PIXEL_FORMAT_DOUBLE , double )
172 // complex<float> should work - but causes all manner of compiler template errors.
173 // maybe need a better compiler, maybe there is a code fix - IMS
174 #undef macro
175 default:
176 return nullptr;
177 }
178 }
179
nplanes()180 unsigned nplanes() const override { return src_->nplanes(); }
ni()181 unsigned ni() const override { return src_->ni(); }
nj()182 unsigned nj() const override { return src_->nj(); }
183
pixel_format()184 enum vil_pixel_format pixel_format() const override
185 { return vil_pixel_format_of(accumT()); }
186
187
188 //: Put the data in this view back into the image source.
put_view(const vil_image_view_base &,unsigned,unsigned)189 bool put_view(const vil_image_view_base& /*im*/, unsigned /*i0*/, unsigned /*j0*/) override
190 {
191 std::cerr << "WARNING: vil_correlate_1d_resource::put_back\n"
192 << "\tYou can't push data back into a correlate filter.\n";
193 return false;
194 }
195
196 //: Extra property information
197 bool get_property(char const* tag, void* property_value = nullptr) const override
198 {
199 if (0==std::strcmp(tag, vil_property_read_only))
200 return property_value ? (*static_cast<bool*>(property_value)) = true : true;
201
202 return src_->get_property(tag, property_value);
203 }
204
205 protected:
206 vil_image_resource_sptr src_;
207 const kernelT* kernel_;
208 std::ptrdiff_t klo_, khi_;
209 vil_convolve_boundary_option start_option_, end_option_;
210 };
211
212 //: Create an image_resource object which correlate kernel[x] x in [k_lo,k_hi] with srcT
213 // \note This function does not reverse the kernel. If you want the
214 // kernel reversed, use vil_convolve_1d instead.
215 // \param kernel should point to tap 0.
216 // \relatesalso vil_image_resource
217 template <class destT, class kernelT, class accumT>
vil_correlate_1d(const vil_image_resource_sptr & src_im,const destT,const kernelT * kernel,std::ptrdiff_t k_lo,std::ptrdiff_t k_hi,const accumT,vil_convolve_boundary_option start_option,vil_convolve_boundary_option end_option)218 vil_image_resource_sptr vil_correlate_1d(
219 const vil_image_resource_sptr& src_im,
220 const destT /*dt*/,
221 const kernelT* kernel, std::ptrdiff_t k_lo, std::ptrdiff_t k_hi,
222 const accumT,
223 vil_convolve_boundary_option start_option,
224 vil_convolve_boundary_option end_option)
225 {
226 return new vil_correlate_1d_resource<kernelT, accumT, destT>(src_im, kernel, k_lo, k_hi, start_option, end_option);
227 }
228
229 #endif // vil_correlate_1d_h_
230