1 // This is mul/vimt/algo/vimt_find_peaks.h
2 #ifndef vimt_find_peaks_h_
3 #define vimt_find_peaks_h_
4 //:
5 // \file
6 // \brief Find peaks in image
7 // \author Tim Cootes, VP (Sept03)
8 
9 #include <vimt/vimt_image_2d_of.h>
10 #include <vgl/vgl_point_2d.h>
11 
12 //: True if pixel at *im is strictly above 8 neighbours
13 template <class T>
vimt_is_peak_3x3(const T * im,std::ptrdiff_t i_step,std::ptrdiff_t j_step)14 inline bool vimt_is_peak_3x3(const T* im, std::ptrdiff_t i_step, std::ptrdiff_t j_step)
15 {
16   T v = *im;
17   return v>im[i_step] &&
18          v>im[-i_step] &&
19          v>im[j_step] &&
20          v>im[-j_step] &&
21          v>im[i_step+j_step] &&
22          v>im[i_step-j_step] &&
23          v>im[j_step-i_step] &&
24          v>im[-i_step-j_step];
25 }
26 
27 //: True if pixel at *im is strictly above its neighbours in a 2*radius+1 neighbourhood
28 template <class T>
vimt_is_peak(const T * im,int radius,std::ptrdiff_t i_step,std::ptrdiff_t j_step)29 inline bool vimt_is_peak(const T* im, int radius, std::ptrdiff_t i_step, std::ptrdiff_t j_step)
30 {
31   T v = *im;
32   for (int i=-radius; i<radius+1; i++)
33     for (int j=-radius; j<radius+1; j++)
34       if (i!=0 || j!=0)
35         if (v<=im[i_step*i+j_step*j]) return false;      // One of the
36   return true;
37 }
38 
39 
40 //: Return image co-ordinates of all points in image strictly above their 8 neighbours
41 // \param clear_list: If true (the default) then empty list before adding new examples
42 template <class T>
43 inline void vimt_find_image_peaks_3x3(std::vector<vgl_point_2d<unsigned> >& peaks,
44                                       const vil_image_view<T>& image,
45                                       unsigned plane=0, bool clear_list=true)
46 {
47   if (clear_list) peaks.resize(0);
48   unsigned ni=image.ni(),nj=image.nj();
49   std::ptrdiff_t istep = image.istep(),jstep=image.jstep();
50   const T* row = image.top_left_ptr()+plane*image.planestep()+istep+jstep;
51   for (unsigned j=1;j<nj-1;++j,row+=jstep)
52   {
53     const T* pixel = row;
54     for (unsigned i=1;i<ni-1;++i,pixel+=istep)
55       if (vimt_is_peak_3x3(pixel,istep,jstep)) peaks.emplace_back(i,j);
56   }
57 }
58 
59 //: Return image co-ordinates of all points in image strictly above their 8 neighbours
60 // \param peak_value: Value at peak
61 // \param clear_list: If true (the default) then empty list before adding new examples
62 template <class T>
63 inline void vimt_find_image_peaks_3x3(std::vector<vgl_point_2d<unsigned> >& peaks,
64                                       std::vector<T>& peak_value,
65                                       const vil_image_view<T>& image,
66                                       unsigned plane=0, bool clear_list=true)
67 {
68   if (clear_list) { peaks.resize(0); peak_value.resize(0); }
69   unsigned ni=image.ni(),nj=image.nj();
70   std::ptrdiff_t istep = image.istep(),jstep=image.jstep();
71   const T* row = image.top_left_ptr()+plane*image.planestep()+istep+jstep;
72   for (unsigned j=1;j<nj-1;++j,row+=jstep)
73   {
74     const T* pixel = row;
75     for (unsigned i=1;i<ni-1;++i,pixel+=istep)
76       if (vimt_is_peak_3x3(pixel,istep,jstep))
77       {
78         peaks.emplace_back(i,j);
79         peak_value.push_back(*pixel);
80       }
81   }
82 }
83 
84 //: Return image co-ordinates of all points in image strictly above their neighbours
85 //  In a 2*radius+1 x 2*radius+1 neighbourhood of pixels (e.g. r=2 equivalent to 5x5; default: r=1)
86 // \param peak_value: Value at peak
87 // \param clear_list: If true (the default) then empty list before adding new examples
88 template <class T>
89 inline void vimt_find_image_peaks(std::vector<vgl_point_2d<unsigned> >& peaks,
90                                   std::vector<T>& peak_value,
91                                   const vil_image_view<T>& image,
92                                   unsigned radius=1,
93                                   unsigned plane=0, bool clear_list=true)
94 {
95   if (clear_list) { peaks.resize(0); peak_value.resize(0); }
96   unsigned ni=image.ni(),nj=image.nj();
97   std::ptrdiff_t istep = image.istep(),jstep=image.jstep();
98   const T* row = image.top_left_ptr()+plane*image.planestep()+radius*istep+radius*jstep;
99   for (unsigned j=radius;j<nj-radius;++j,row+=jstep)
100   {
101     const T* pixel = row;
102     for (unsigned i=radius;i<ni-radius;++i,pixel+=istep)
103       if (vimt_is_peak(pixel,radius,istep,jstep))
104       {
105         peaks.emplace_back(i,j);
106         peak_value.push_back(*pixel);
107       }
108   }
109 }
110 
111 //: Return image co-ordinates of all points in image strictly above their neighbours
112 //  In a 2*radius+1 x 2*radius+1 neighbourhood of pixels (e.g. r=2 equivalent to 5x5; default: r=1)
113 //  Additionally, only peaks of the value higher than threshold (thresh) are returned.
114 // \param peak_value: Value at peak
115 // \param clear_list: If true (the default) then empty list before adding new examples
116 template <class T>
117 inline void vimt_find_image_peaks(std::vector<vgl_point_2d<unsigned> >& peaks,
118                                   std::vector<T>& peak_value,
119                                   const vil_image_view<T>& image,
120                                   T thresh, unsigned radius=1,
121                                   unsigned plane=0, bool clear_list=true)
122 {
123   if (clear_list) { peaks.resize(0); peak_value.resize(0); }
124   unsigned ni=image.ni(),nj=image.nj();
125   std::ptrdiff_t istep = image.istep(),jstep=image.jstep();
126   // Getting to the location of the starting point in the image (radius,radius)
127   const T* row = image.top_left_ptr()+plane*image.planestep()+radius*istep+radius*jstep;
128   for (unsigned j=radius;j<nj-radius;++j,row+=jstep)
129   {
130     const T* pixel = row;
131     for (unsigned i=radius;i<ni-radius;++i,pixel+=istep)
132     {
133       if (*pixel>thresh)
134       {
135         if (vimt_is_peak(pixel,radius,istep,jstep))
136         {
137           peaks.emplace_back(i,j);
138           peak_value.push_back(*pixel);
139         }
140       }
141     }
142   }
143 }
144 
145 //: Return world co-ordinates of all points in image strictly above their 8 neighbours
146 // \param clear_list: If true (the default) then empty list before adding new examples
147 template <class T>
148 inline void vimt_find_world_peaks_3x3(std::vector<vgl_point_2d<double> >& peaks,
149                                       const vimt_image_2d_of<T>& image,
150                                       unsigned plane=0, bool clear_list=true)
151 {
152   if (clear_list) peaks.resize(0);
153   const vil_image_view<T>& im = image.image();
154   vimt_transform_2d im2w = image.world2im().inverse();
155   unsigned ni=im.ni(),nj=im.nj();
156   std::ptrdiff_t istep = im.istep(),jstep=im.jstep();
157   const T* row = im.top_left_ptr()+plane*im.planestep()+istep+jstep;
158   for (unsigned j=1;j<nj-1;++j,row+=jstep)
159   {
160     const T* pixel = row;
161     for (unsigned i=1;i<ni-1;++i,pixel+=istep)
162       if (vimt_is_peak_3x3(pixel,istep,jstep)) peaks.push_back(im2w(i,j));
163   }
164 }
165 
166 //: Return world co-ordinates of all points in image strictly above their 8 neighbours
167 // \param peak_pos: Position of each peak
168 // \param peak_value: Value at peak
169 // \param clear_list: If true (the default) then empty list before adding new examples
170 template <class T>
171 inline void vimt_find_world_peaks_3x3(std::vector<vgl_point_2d<double> >& peak_pos,
172                                       std::vector<T>& peak_value,
173                                       const vimt_image_2d_of<T>& image,
174                                       unsigned plane=0, bool clear_list=true)
175 {
176   if (clear_list) { peak_pos.resize(0); peak_value.resize(0); }
177   const vil_image_view<T>& im = image.image();
178   vimt_transform_2d im2w = image.world2im().inverse();
179   unsigned ni=im.ni(),nj=im.nj();
180   std::ptrdiff_t istep = im.istep(),jstep=im.jstep();
181   const T* row = im.top_left_ptr()+plane*im.planestep()+istep+jstep;
182   for (unsigned j=1;j<nj-1;++j,row+=jstep)
183   {
184     const T* pixel = row;
185     for (unsigned i=1;i<ni-1;++i,pixel+=istep)
186       if (vimt_is_peak_3x3(pixel,istep,jstep))
187       {
188         peak_pos.push_back(im2w(i,j));
189         peak_value.push_back(*pixel);
190       }
191   }
192 }
193 
194 //: Return image co-ordinates of maximum value in image
195 //  (Or first one found if multiple equivalent maxima)
196 template <class T>
197 inline
198 vgl_point_2d<unsigned> vimt_find_max(const vil_image_view<T>& im, unsigned plane=0)
199 {
200   vgl_point_2d<unsigned> p(0,0);
201   T max_val = im(0,0,plane);
202   unsigned ni=im.ni(),nj=im.nj();
203   std::ptrdiff_t istep = im.istep(),jstep=im.jstep();
204   const T* row = im.top_left_ptr()+plane*im.planestep();
205   for (unsigned j=0;j<nj;++j,row+=jstep)
206   {
207     const T* pixel = row;
208     for (unsigned i=0;i<ni;++i,pixel+=istep)
209       if (*pixel>max_val)
210       {
211         max_val = *pixel;
212         p = vgl_point_2d<unsigned>(i,j);
213       }
214   }
215   return p;
216 }
217 
218 //: Return world co-ordinates of maximum value in image
219 //  (Or first one found if multiple equivalent maxima)
220 template <class T>
221 inline
222 vgl_point_2d<double> vimt_find_max(const vimt_image_2d_of<T>& image,unsigned plane=0)
223 {
224   vgl_point_2d<unsigned> im_p = vimt_find_max(image.image(),plane);
225   return image.world2im().inverse()(im_p.x(),im_p.y());
226 }
227 
228 
229 #endif // vimt_find_peaks_h_
230