00001
00002 #ifndef vil_suppress_non_max_edges_txx_
00003 #define vil_suppress_non_max_edges_txx_
00004
00005
00006
00007
00008
00009 #include "vil_suppress_non_max_edges.h"
00010 #include <vil/vil_bilin_interp.h>
00011 #include <vil/vil_fill.h>
00012 #include <vcl_cmath.h>
00013 #include <vcl_cassert.h>
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 template<class srcT, class destT>
00024 void vil_suppress_non_max_edges(const vil_image_view<srcT>& grad_i,
00025 const vil_image_view<srcT>& grad_j,
00026 double grad_mag_threshold,
00027 vil_image_view<destT>& grad_mag)
00028 {
00029 assert(grad_i.nplanes()==grad_j.nplanes());
00030 assert(grad_i.nplanes()==1);
00031 unsigned ni = grad_i.ni(), nj = grad_i.nj();
00032 assert(ni>2 && nj>2);
00033 assert(grad_j.ni()==ni && grad_j.nj()==nj);
00034 grad_mag.set_size(ni,nj,1);
00035
00036
00037 vil_fill_col(grad_mag,0,destT(0));
00038 vil_fill_col(grad_mag,1,destT(0));
00039 vil_fill_col(grad_mag,ni-1,destT(0));
00040 vil_fill_col(grad_mag,ni-2,destT(0));
00041 vil_fill_row(grad_mag,0,destT(0));
00042 vil_fill_row(grad_mag,1,destT(0));
00043 vil_fill_row(grad_mag,nj-1,destT(0));
00044 vil_fill_row(grad_mag,nj-2,destT(0));
00045
00046 const vcl_ptrdiff_t gi_istep = grad_i.istep(), gi_jstep = grad_i.jstep();
00047 const vcl_ptrdiff_t gj_istep = grad_j.istep(), gj_jstep = grad_j.jstep();
00048 const vcl_ptrdiff_t gm_istep = grad_mag.istep(), gm_jstep = grad_mag.jstep();
00049
00050 const srcT * gi_data = &grad_i(0,0);
00051 const srcT * gj_data = &grad_j(0,0);
00052 const srcT * gi_row = &grad_i(2,2);
00053 const srcT * gj_row = &grad_j(2,2);
00054 destT * gm_row = &grad_mag(2,2);
00055 unsigned ihi=ni-3;
00056 unsigned jhi=nj-3;
00057
00058 for (unsigned j=2; j<=jhi; ++j, gi_row+=gi_jstep, gj_row+=gj_jstep,
00059 gm_row+=gm_jstep)
00060 {
00061 const srcT* pgi = gi_row;
00062 const srcT* pgj = gj_row;
00063 destT *pgm = gm_row;
00064 for (unsigned i=2; i<=ihi; ++i, pgi+=gi_istep, pgj+=gj_istep,
00065 pgm+=gm_istep)
00066 {
00067 double gmag=vcl_sqrt(double(pgi[0]*pgi[0] + pgj[0]*pgj[0]));
00068 if (gmag<grad_mag_threshold) *pgm=0;
00069 else
00070 {
00071 double dx=pgi[0]/gmag;
00072 double dy=pgj[0]/gmag;
00073
00074 double gx1=vil_bilin_interp_unsafe(i+dx,j+dy,gi_data,gi_istep,gi_jstep);
00075 double gy1=vil_bilin_interp_unsafe(i+dx,j+dy,gj_data,gj_istep,gj_jstep);
00076 if (dx*gx1+dy*gy1>gmag) *pgm=0;
00077 else
00078 {
00079
00080 double gx2=vil_bilin_interp_unsafe(i-dx,j-dy,gi_data,gi_istep,gi_jstep);
00081 double gy2=vil_bilin_interp_unsafe(i-dx,j-dy,gj_data,gj_istep,gj_jstep);
00082 if (dx*gx2+dy*gy2>gmag) *pgm=0;
00083 else
00084 *pgm = destT(gmag);
00085 }
00086 }
00087 }
00088 }
00089 }
00090
00091 namespace {
00092
00093
00094
00095 double interpolate_parabola(double y_1, double y_0, double y_2,
00096 double& y_peak)
00097 {
00098 y_peak = y_0;
00099 double diff1 = y_2 - y_1;
00100 double diff2 = 2 * y_0 - y_1 - y_2;
00101
00102 if (diff2 == 0.0)
00103 return 0.0;
00104 y_peak += diff1 * diff1 / (8 * diff2);
00105 return diff1 / (2 * diff2);
00106 }
00107 }
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 template<class srcT, class destT>
00139 void vil_suppress_non_max_edges_subpixel(const vil_image_view<srcT>& grad_i,
00140 const vil_image_view<srcT>& grad_j,
00141 double grad_mag_threshold,
00142 vil_image_view<destT>& grad_moo)
00143 {
00144 assert(grad_i.nplanes()==grad_j.nplanes());
00145 assert(grad_i.nplanes()==1);
00146 unsigned ni = grad_i.ni(), nj = grad_i.nj();
00147 assert(ni>2 && nj>2);
00148 assert(grad_j.ni()==ni && grad_j.nj()==nj);
00149 grad_moo.set_size(ni,nj,3);
00150
00151
00152 vil_fill_col(grad_moo,0,destT(0));
00153 vil_fill_col(grad_moo,1,destT(0));
00154 vil_fill_col(grad_moo,ni-1,destT(0));
00155 vil_fill_col(grad_moo,ni-2,destT(0));
00156 vil_fill_row(grad_moo,0,destT(0));
00157 vil_fill_row(grad_moo,1,destT(0));
00158 vil_fill_row(grad_moo,nj-1,destT(0));
00159 vil_fill_row(grad_moo,nj-2,destT(0));
00160
00161 const vcl_ptrdiff_t gi_istep = grad_i.istep(), gi_jstep = grad_i.jstep();
00162 const vcl_ptrdiff_t gj_istep = grad_j.istep(), gj_jstep = grad_j.jstep();
00163 const vcl_ptrdiff_t gm_istep = grad_moo.istep(), gm_jstep = grad_moo.jstep();
00164 const vcl_ptrdiff_t gm_pstep = grad_moo.planestep();
00165
00166 const srcT * gi_data = &grad_i(0,0);
00167 const srcT * gj_data = &grad_j(0,0);
00168 const srcT * gi_row = &grad_i(2,2);
00169 const srcT * gj_row = &grad_j(2,2);
00170 destT * gm_row = &grad_moo(2,2);
00171 unsigned ihi=ni-3;
00172 unsigned jhi=nj-3;
00173
00174 for (unsigned j=2; j<=jhi; ++j, gi_row+=gi_jstep, gj_row+=gj_jstep,
00175 gm_row+=gm_jstep)
00176 {
00177 const srcT* pgi = gi_row;
00178 const srcT* pgj = gj_row;
00179 destT *pgm = gm_row;
00180 for (unsigned i=2; i<=ihi; ++i, pgi+=gi_istep, pgj+=gj_istep,
00181 pgm+=gm_istep)
00182 {
00183 double gmag=vcl_sqrt(double(pgi[0]*pgi[0] + pgj[0]*pgj[0]));
00184 if (gmag<grad_mag_threshold){
00185 *pgm=0;
00186 *(pgm+gm_pstep)=0;
00187 *(pgm+2*gm_pstep)=0;
00188 }
00189 else
00190 {
00191 double dx=pgi[0]/gmag;
00192 double dy=pgj[0]/gmag;
00193
00194 double gx1=vil_bilin_interp_unsafe(i+dx,j+dy,gi_data,gi_istep,gi_jstep);
00195 double gy1=vil_bilin_interp_unsafe(i+dx,j+dy,gj_data,gj_istep,gj_jstep);
00196 double g1mag = dx*gx1+dy*gy1;
00197 if (g1mag>gmag){
00198 *pgm=0;
00199 *(pgm+gm_pstep)=0;
00200 *(pgm+2*gm_pstep)=0;
00201 }
00202 else
00203 {
00204
00205 double gx2=vil_bilin_interp_unsafe(i-dx,j-dy,gi_data,gi_istep,gi_jstep);
00206 double gy2=vil_bilin_interp_unsafe(i-dx,j-dy,gj_data,gj_istep,gj_jstep);
00207 double g2mag = dx*gx2+dy*gy2;
00208 if (g2mag>gmag){
00209 *pgm=0;
00210 *(pgm+gm_pstep)=0;
00211 *(pgm+2*gm_pstep)=0;
00212 }
00213 else
00214 {
00215
00216 double peak;
00217 double offset = interpolate_parabola(g2mag, gmag, g1mag, peak);
00218 *pgm = destT(peak);
00219 *(pgm+gm_pstep) = destT(vcl_atan2(dy,dx));
00220 *(pgm+2*gm_pstep) = destT(offset);
00221 }
00222 }
00223 }
00224 }
00225 }
00226 }
00227
00228 #undef VIL_SUPPRESS_NON_MAX_EDGES_INSTANTIATE
00229 #define VIL_SUPPRESS_NON_MAX_EDGES_INSTANTIATE(srcT, destT) \
00230 template void vil_suppress_non_max_edges(const vil_image_view<srcT >& grad_i,\
00231 const vil_image_view<srcT >& grad_j,\
00232 double grad_mag_threshold,\
00233 vil_image_view<destT >& grad_mag);\
00234 template void vil_suppress_non_max_edges_subpixel(const vil_image_view<srcT >& grad_i,\
00235 const vil_image_view<srcT >& grad_j,\
00236 double grad_mag_threshold,\
00237 vil_image_view<destT >& grad_moo)
00238
00239 #endif // vil_suppress_non_max_edges_txx_