00001
00002 #ifndef vil_gauss_filter_txx_
00003 #define vil_gauss_filter_txx_
00004
00005
00006
00007
00008
00009 #include "vil_gauss_filter.h"
00010 #include <vil/vil_transpose.h>
00011 #include <vil/algo/vil_convolve_1d.h>
00012 #include <vcl_iostream.h>
00013 #include <vcl_cassert.h>
00014
00015
00016
00017
00018 inline unsigned char vl_round(double x, unsigned char ) { return (unsigned char )(x<0?x-0.5:x+0.5);}
00019 inline signed char vl_round(double x, signed char ) { return ( signed char )(x<0?x-0.5:x+0.5);}
00020 inline unsigned short vl_round(double x, unsigned short ) { return (unsigned short)(x<0?x-0.5:x+0.5);}
00021 inline signed short vl_round(double x, signed short ) { return ( signed short)(x<0?x-0.5:x+0.5);}
00022 inline unsigned int vl_round(double x, unsigned int ) { return (unsigned int )(x<0?x-0.5:x+0.5);}
00023 inline signed int vl_round(double x, signed int ) { return ( signed int )(x<0?x-0.5:x+0.5);}
00024 inline unsigned long vl_round(double x, unsigned long ) { return (unsigned long )(x<0?x-0.5:x+0.5);}
00025 inline signed long vl_round(double x, signed long ) { return ( signed long )(x<0?x-0.5:x+0.5);}
00026 inline double vl_round(double x, double ) { return x; }
00027 inline float vl_round(double x, float ) { return (float)x; }
00028
00029
00030
00031
00032
00033
00034
00035 template <class srcT, class destT>
00036 void vil_gauss_filter_5tap(const srcT* src_im, vcl_ptrdiff_t src_istep, vcl_ptrdiff_t src_jstep,
00037 destT* dest_im, vcl_ptrdiff_t dest_istep, vcl_ptrdiff_t dest_jstep,
00038 unsigned nx, unsigned ny, const vil_gauss_filter_5tap_params& params,
00039 destT* work, vcl_ptrdiff_t work_jstep)
00040 {
00041 assert (nx > 3 && ny > 3);
00042
00043
00044
00045 for (unsigned int y=0;y<ny;y++)
00046 {
00047 destT* work_row = work + y*work_jstep;
00048 const srcT* src_col3 = src_im + y*src_jstep;
00049 const srcT* src_col2 = src_col3 - src_istep;
00050 const srcT* src_col1 = src_col3 - 2 * src_istep;
00051 const srcT* src_col4 = src_col3 + src_istep;
00052 const srcT* src_col5 = src_col3 + 2 * src_istep;
00053
00054 int x;
00055 int nx2 = nx-2;
00056 for (x=2;x<nx2;x++)
00057 work_row[x] = vl_round( params.filt2() * src_col1[x*src_istep]
00058 + params.filt1() * src_col2[x*src_istep]
00059 + params.filt0() * src_col3[x*src_istep]
00060 + params.filt1() * src_col4[x*src_istep]
00061 + params.filt2() * src_col5[x*src_istep], (destT)0);
00062
00063
00064 work_row[0] = vl_round( params.filt_edge0() * src_col3[0]
00065 + params.filt_edge1() * src_col4[0]
00066 + params.filt_edge2() * src_col5[0], (destT)0);
00067
00068 work_row[1] = vl_round( params.filt_pen_edge_n1() * src_col2[src_istep]
00069 + params.filt_pen_edge0() * src_col3[src_istep]
00070 + params.filt_pen_edge1() * src_col4[src_istep]
00071 + params.filt_pen_edge2() * src_col5[src_istep], (destT)0);
00072
00073 work_row[nx-2] = vl_round( params.filt_pen_edge2() * src_col1[(nx-2)*src_istep]
00074 + params.filt_pen_edge1() * src_col2[(nx-2)*src_istep]
00075 + params.filt_pen_edge0() * src_col3[(nx-2)*src_istep]
00076 + params.filt_pen_edge_n1() * src_col4[(nx-2)*src_istep], (destT)0);
00077
00078 work_row[nx-1] = vl_round( params.filt_edge2() * src_col1[(nx-1)*src_istep]
00079 + params.filt_edge1() * src_col2[(nx-1)*src_istep]
00080 + params.filt_edge0() * src_col3[(nx-1)*src_istep], (destT)0);
00081 }
00082
00083
00084
00085 for (unsigned int y=2;y<ny-2;y++)
00086 {
00087 destT* dest_row = dest_im + y*dest_jstep;
00088
00089 const destT* work_row3 = work + y*work_jstep;
00090 const destT* work_row2 = work_row3 - work_jstep;
00091 const destT* work_row1 = work_row3 - 2 * work_jstep;
00092 const destT* work_row4 = work_row3 + work_jstep;
00093 const destT* work_row5 = work_row3 + 2 * work_jstep;
00094
00095 for (unsigned int x=0; x<nx; x++)
00096 dest_row[x*dest_istep] = vl_round( params.filt2() * work_row1[x]
00097 + params.filt1() * work_row2[x]
00098 + params.filt0() * work_row3[x]
00099 + params.filt1() * work_row4[x]
00100 + params.filt2() * work_row5[x], (destT)0);
00101 }
00102
00103
00104
00105 const destT* work_row_bottom_1 = work;
00106 const destT* work_row_bottom_2 = work_row_bottom_1 + work_jstep;
00107 const destT* work_row_bottom_3 = work_row_bottom_1 + 2 * work_jstep;
00108 const destT* work_row_bottom_4 = work_row_bottom_1 + 3 * work_jstep;
00109
00110 const destT* work_row_top_5 = work + (ny-1) * work_jstep;
00111 const destT* work_row_top_4 = work_row_top_5 - work_jstep;
00112 const destT* work_row_top_3 = work_row_top_5 - 2 * work_jstep;
00113 const destT* work_row_top_2 = work_row_top_5 - 3 * work_jstep;
00114
00115 destT* dest_row_top = dest_im + (ny-1) * dest_jstep;
00116 destT* dest_row_next_top = dest_row_top - dest_jstep;
00117 destT* dest_row_bottom = dest_im;
00118 destT* dest_row_next_bottom = dest_row_bottom + dest_jstep;
00119
00120 for (unsigned int x=0;x<nx;x++)
00121 {
00122 dest_row_top[x*dest_istep] = vl_round( params.filt_edge0() * work_row_top_5[x]
00123 + params.filt_edge1() * work_row_top_4[x]
00124 + params.filt_edge2() * work_row_top_3[x], (destT)0);
00125
00126 dest_row_next_top[x*dest_istep] = vl_round( params.filt_pen_edge2() * work_row_top_2[x]
00127 + params.filt_pen_edge1() * work_row_top_3[x]
00128 + params.filt_pen_edge0() * work_row_top_4[x]
00129 + params.filt_pen_edge_n1()*work_row_top_5[x], (destT)0);
00130
00131 dest_row_next_bottom[x*dest_istep] = vl_round( params.filt_pen_edge2() * work_row_bottom_4[x]
00132 + params.filt_pen_edge1() * work_row_bottom_3[x]
00133 + params.filt_pen_edge0() * work_row_bottom_2[x]
00134 + params.filt_pen_edge_n1()*work_row_bottom_1[x], (destT)0);
00135
00136 dest_row_bottom[x*dest_istep] = vl_round( params.filt_edge2() * work_row_bottom_3[x]
00137 + params.filt_edge1() * work_row_bottom_2[x]
00138 + params.filt_edge0() * work_row_bottom_1[x], (destT)0);
00139 }
00140
00141 }
00142
00143 template <class srcT, class destT>
00144 void vil_gauss_filter_5tap(const vil_image_view<srcT>& src_im,
00145 vil_image_view<destT>& dest_im,
00146 const vil_gauss_filter_5tap_params& params,
00147 vil_image_view<destT>& work)
00148 {
00149 unsigned ni = src_im.ni();
00150 unsigned nj = src_im.nj();
00151 unsigned n_planes = src_im.nplanes();
00152 dest_im.set_size(ni, nj, n_planes);
00153 work.set_size(ni,nj,1);
00154 assert (work.jstep() == (vcl_ptrdiff_t)ni);
00155
00156 if (ni > 3 && nj > 3)
00157
00158 for (unsigned p=0;p<n_planes;++p)
00159 vil_gauss_filter_5tap(&src_im(0,0,p), src_im.istep(), src_im.jstep(),
00160 &dest_im(0,0,p), dest_im.istep(), dest_im.jstep(), ni,nj,
00161 params, work.top_left_ptr(), work.jstep());
00162 else
00163 {
00164 if (ni==0 || nj==0) return;
00165 if (ni==1)
00166 {
00167 for (unsigned p=0;p<n_planes;++p)
00168 for (unsigned j=0;j<nj;++j)
00169 work(0,j,p) = vl_round(src_im(0,j,p), destT());
00170 }
00171 else if (ni==2)
00172 {
00173 const double k0 = params.filt0()/(params.filt0() + params.filt1());
00174 const double k1 = params.filt1()/(params.filt0() + params.filt1());
00175 for (unsigned p=0;p<n_planes;++p)
00176 for (unsigned j=0;j<nj;++j)
00177 {
00178 work(0,j,p) = vl_round(k0*src_im(0,j,p) + k1*src_im(1,j,p), destT());
00179 work(1,j,p) = vl_round(k1*src_im(0,j,p) + k0*src_im(1,j,p), destT());
00180 }
00181 }
00182 else if (ni==3)
00183 {
00184 const double ke0 = params.filt0()/(params.filt0() + params.filt1());
00185 const double ke1 = params.filt1()/(params.filt0() + params.filt1());
00186 const double k0 = params.filt0()/(params.filt0() + 2.0*params.filt1());
00187 const double k1 = params.filt1()/(params.filt0() + 2.0*params.filt1());
00188 for (unsigned p=0;p<n_planes;++p)
00189 for (unsigned j=0;j<nj;++j)
00190 {
00191 work(0,j,p) = vl_round(ke0*src_im(0,j,p) + ke1*src_im(1,j,p) , destT());
00192 work(1,j,p) = vl_round(k1 *src_im(0,j,p) + k0 *src_im(1,j,p) + k1 *src_im(2,j,p), destT());
00193 work(2,j,p) = vl_round( ke1*src_im(1,j,p) + ke0*src_im(2,j,p), destT());
00194 }
00195 }
00196 else if (ni==4)
00197 {
00198 const double ke0 = params.filt_edge0();
00199 const double ke1 = params.filt_edge1();
00200 const double ke2 = params.filt_edge2();
00201 const double kpn1= params.filt_pen_edge_n1();
00202 const double kp0 = params.filt_pen_edge0();
00203 const double kp1 = params.filt_pen_edge1();
00204 const double kp2 = params.filt_pen_edge2();
00205 for (unsigned p=0;p<n_planes;++p)
00206 for (unsigned j=0;j<nj;++j)
00207 {
00208 work(0,j,p) = vl_round(ke0 *src_im(0,j,p) + ke1*src_im(1,j,p) + ke2*src_im(2,j,p) , destT());
00209 work(1,j,p) = vl_round(kpn1*src_im(0,j,p) + kp0*src_im(1,j,p) + kp1*src_im(2,j,p) + kp2 *src_im(3,j,p), destT());
00210 work(2,j,p) = vl_round(kp2 *src_im(0,j,p) + kp1*src_im(1,j,p) + kp0*src_im(2,j,p) + kpn1*src_im(3,j,p), destT());
00211 work(3,j,p) = vl_round( ke2*src_im(1,j,p) + ke1*src_im(2,j,p) + ke0 *src_im(3,j,p), destT());
00212 }
00213 }
00214 else if (ni>4)
00215 {
00216 double kernel[5];
00217 kernel[0] = kernel[4] = params.filt2();
00218 kernel[1] = kernel[3] = params.filt1();
00219 kernel[2] = params.filt0();
00220 vil_convolve_1d(src_im,work, kernel+2, -2, +2,
00221 double(), vil_convolve_trim, vil_convolve_trim);
00222 }
00223 if (nj==1)
00224 dest_im=work;
00225 else if (nj==2)
00226 {
00227 const double k0 = params.filt0()/(params.filt0() + params.filt1());
00228 const double k1 = params.filt1()/(params.filt0() + params.filt1());
00229 for (unsigned p=0;p<n_planes;++p)
00230 for (unsigned i=0;i<ni;++i)
00231 {
00232 dest_im(i,0,p) = vl_round(k0 * work(i,0,p) + k1 * work(i,1,p), destT());
00233 dest_im(i,1,p) = vl_round(k1 * work(i,0,p) + k0 * work(i,1,p), destT());
00234 }
00235 }
00236 else if (nj==3)
00237 {
00238 const double ke0 = params.filt0()/(params.filt0() + params.filt1());
00239 const double ke1 = params.filt1()/(params.filt0() + params.filt1());
00240 const double k0 = params.filt0()/(params.filt0() + 2.0*params.filt1());
00241 const double k1 = params.filt1()/(params.filt0() + 2.0*params.filt1());
00242 for (unsigned p=0;p<n_planes;++p)
00243 for (unsigned i=0;i<ni;++i)
00244 {
00245 dest_im(i,0,p) = vl_round(ke0*work(i,0,p) + ke1*work(i,1,p) , destT());
00246 dest_im(i,1,p) = vl_round(k1 *work(i,0,p) + k0 *work(i,1,p) + k1 *work(i,2,p), destT());
00247 dest_im(i,2,p) = vl_round( ke1*work(i,1,p) + ke0*work(i,2,p), destT());
00248 }
00249 }
00250 else if (nj==4)
00251 {
00252 const double ke0 = params.filt_edge0();
00253 const double ke1 = params.filt_edge1();
00254 const double ke2 = params.filt_edge2();
00255 const double kpn1= params.filt_pen_edge_n1();
00256 const double kp0 = params.filt_pen_edge0();
00257 const double kp1 = params.filt_pen_edge1();
00258 const double kp2 = params.filt_pen_edge2();
00259 for (unsigned p=0;p<n_planes;++p)
00260 for (unsigned i=0;i<ni;++i)
00261 {
00262 dest_im(i,0,p) = vl_round(ke0 *work(i,0,p) + ke1*work(i,1,p) + ke2*work(i,2,p) , destT());
00263 dest_im(i,1,p) = vl_round(kpn1*work(i,0,p) + kp0*work(i,1,p) + kp1*work(i,2,p) + kp2 *work(i,3,p), destT());
00264 dest_im(i,2,p) = vl_round(kp2 *work(i,0,p) + kp1*work(i,1,p) + kp0*work(i,2,p) + kpn1*work(i,3,p), destT());
00265 dest_im(i,3,p) = vl_round( ke2*work(i,1,p) + ke1*work(i,2,p) + ke0 *work(i,3,p), destT());
00266 }
00267 }
00268 else if (nj>4)
00269 {
00270 double kernel[5];
00271 kernel[0] = kernel[4] = params.filt2();
00272 kernel[1] = kernel[3] = params.filt1();
00273 kernel[2] = params.filt0();
00274 vil_image_view<destT> rdest = vil_transpose(dest_im);
00275 vil_convolve_1d(vil_transpose(work), rdest, kernel+2, -2, +2,
00276 double(), vil_convolve_trim, vil_convolve_trim);
00277 }
00278 }
00279
00280 #if 0
00281 vsl_indent_inc(vcl_cout);
00282 vcl_cout << vsl_indent() << "Work image B\n";
00283 workb_.print_all(vcl_cout);
00284 vsl_indent_dec(vcl_cout);
00285 #endif
00286 }
00287
00288 #undef VIL_GAUSS_FILTER_INSTANTIATE
00289 #define VIL_GAUSS_FILTER_INSTANTIATE(srcT, destT) \
00290 template void vil_gauss_filter_5tap(const vil_image_view<srcT >& src_im, \
00291 vil_image_view<destT >& dest_im, \
00292 const vil_gauss_filter_5tap_params& params, \
00293 vil_image_view<destT >& work)
00294
00295 #endif // vil_gauss_filter_txx_