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