1 /**
2 Module introduces Randomized Hough Transform implementation.
3
4 Example:
5 ----
6
7 // Load an image from filesystem.
8 Image image = imread("/path/to/image.png");
9
10 // Calculate edges by use of canny algorithm
11 auto edges = image
12 .sliced
13 .rgb2gray
14 .as!float
15 .slice
16 .conv(gaussian!float(1.0f, 3, 3))
17 .canny!ubyte(100);
18
19 // Setup RHT line extraction context
20 auto lines = RhtLines().epouchs(50).iterations(250).minCurve(25);
21
22 // Lazily iterate and detect lines in pre-processed image
23 foreach(line; lines(edges)) {
24 // do stuff with lines..
25 }
26 ----
27
28 For more elaborated module description visit the $(LINK2 http://dcv.dlang.io/?loc=example_rht,RHT example).
29
30 Copyright: Copyright © 2016, Dmitry Olshansky
31
32 Authors: Dmitry Olshansky
33
34 License: $(LINK3 http://www.boost.org/LICENSE_1_0.txt, Boost Software License - Version 1.0).
35 */
36
37 module dcv.features.rht;
38
39 import std.typecons, std.range.primitives;
40
41 import mir.ndslice;
42
43 /++
44 A template that bootstraps a full Randomized Hough transform implementation.
45 The basic primitives required are as follows.
46
47 Types: Curve and Key tuples that define curve parameters and accumulator key for the curve.
48 Functions: curveKey - accumulator key for a curve, onCurve - test if a point is on curve,
49 fitCurve - fit a curve to a given random access range of points.
50 +/
51 mixin template BaseRht()
52 {
53 alias This = typeof(this);
54 alias Point = Tuple!(int, "x", int, "y");
55 int _thrd = 2; // threshold for a maximum to be detected in accumulator
56 int _epouchs = 50; // number of epouchs to iterate (attempts to find a shape)
57 int _iters = 1000; // iterations in each epouch
58 int _minCurve = 25; // minimal curve length in pixels
59
60 static auto toInt(double a)
61 {
62 import std.math;
63
64 return cast(int)lrint(a);
65 }
66
67 // invalid (or empty) curve if any of parameters is NaN
68 static bool isInvalidCurve(Curve c)
69 {
70 import std.math;
71
72 foreach (v; c.tupleof)
73 if (v.isNaN)
74 return true;
75 return false;
76 }
77
78 /// Set threshold for a curve to be considered in an accumulator
79 ref threshold(int threshold)
80 {
81 _thrd = threshold;
82 return this;
83 }
84
85 ref epouchs(int epouchs)
86 {
87 _epouchs = epouchs;
88 return this;
89 }
90
91 ref iterations(int iters)
92 {
93 _iters = iters;
94 return this;
95 }
96
97 ref minCurve(int minimalCurve)
98 {
99 _minCurve = minimalCurve;
100 return this;
101 }
102
103 static auto opCall()
104 {
105 This zis;
106 return zis;
107 }
108
109 /// Run RHT using non-zero points in image as edge points.
110 auto opCall(T, SliceKind kind)(Slice!(kind, [2], T*) image)
111 {
112 import std.array : appender;
113
114 auto points = appender!(Point[])();
115 int x, y = 0;
116
117 foreach(row; image)
118 {
119 x = 0;
120 foreach(e; row)
121 {
122 if (e > 0)
123 points.put(Point(x, y));
124 ++x;
125 }
126 ++y;
127 }
128 return this.opCall(image, points.data);
129 }
130
131 /// Run RHT using prepopullated array of edge points (that may be filtered beforehand).
132 auto opCall(T, SliceKind kind, Range)(Slice!(kind, [2], T*) image, Range points) if (isInputRange!Range)
133 {
134 auto r = RhtRange!(T, kind, ElementType!Range)(this, image, points);
135 r.popFront(); // prime the detection process
136 return r;
137 }
138
139 static struct RhtRange(T, SliceKind kind, P)
140 {
141 private:
142 import std.container, std.random;
143
144 This _rht; // RHT struct with key primitives and parameters
145 Slice!(kind, [2], T*) _image; // image with edge points
146 Tuple!(Curve, int)[Key] _accum; // accumulator, parametrized on Key/Curve tuples
147 Array!P _points; // extracted edge points
148 int _epouch; // current epouch of iteration
149 Curve _current; // current detected curve
150 Xorshift rng;
151 this(Range)(This rht, Slice!(kind, [2], T*) image, Range points)
152 {
153 _rht = rht;
154 _image = image;
155 _points = make!(Array!P)(points);
156 rng = Xorshift(unpredictableSeed);
157 }
158
159 void accumulate(Curve curve)
160 {
161 auto key = _rht.curveKey(curve);
162 // if have a curve with the same key
163 // replace it with weighted average
164 if (auto p = (key in _accum))
165 {
166 auto prior = (*p)[0];
167 auto w = (*p)[1];
168 Curve weighted;
169 foreach (i, _; prior.tupleof)
170 {
171 weighted[i] = (prior[i] * w + curve[i]) / (w + 1.0);
172 }
173 auto newKey = _rht.curveKey(weighted);
174 if (key != newKey)
175 {
176 _accum.remove(key);
177 _accum[newKey] = tuple(weighted, w + 1);
178 }
179 else
180 *p = tuple(weighted, w + 1);
181 }
182 // simply add current curve to accum
183 else
184 {
185 _accum[key] = tuple(curve, 1);
186 }
187 }
188
189 Curve bestCurve()
190 {
191 Curve best;
192 int maxW = 0;
193 foreach (_, v; _accum)
194 {
195 if (v[1] > maxW)
196 {
197 maxW = v[1];
198 best = v[0];
199 }
200 }
201 return best;
202 }
203
204 public:
205 Curve front()
206 {
207 return _current;
208 }
209
210 bool empty()
211 {
212 return isInvalidCurve(_current);
213 }
214
215 void popFront()
216 {
217 import std.algorithm, std.range;
218
219 foreach (e; _epouch .. _rht._epouchs)
220 {
221 _accum.clear();
222 Curve best;
223 foreach (i; 0 .. _rht._iters)
224 {
225 if (_points.length < _rht.sampleSize)
226 break;
227 // TODO: avoid heap allocation
228 auto sample = randomSample(_points[], _rht.sampleSize, &rng).array;
229 auto curve = _rht.fitCurve(_image, sample);
230 if (!isInvalidCurve(curve))
231 accumulate(curve);
232 }
233 best = bestCurve();
234 if (isInvalidCurve(best))
235 continue;
236 auto newPoints = make!Array(_points[].filter!(x => !_rht.onCurve(best, x)));
237 if (_points.length - newPoints.length > _rht._minCurve)
238 {
239 // remove fitted curve from the set of points
240 copy(newPoints[], _points[0 .. newPoints.length]);
241 _points.length = newPoints.length;
242 _current = best;
243 _epouch = e + 1; // skip current epouch
244 return; // stop prematurely
245 }
246 }
247 // spinned through all of epouchs - no more curves to detect
248 _current = Curve.init;
249 }
250
251 /// Resulting points that are not fitted to any shape yet
252 @property auto points()
253 {
254 return _points;
255 }
256 }
257 }
258
259 /// Randomized Hough Transform for lines
260 struct RhtLines
261 {
262 alias Key = Tuple!(double, double);
263 alias Curve = Tuple!(double, "m", double, "b");
264 enum sampleSize = 2; // sample points to use in each iteration
265 double _angleTol = 1.0; // tolerance in angle approximation (deg)
266 double _interceptTol = 1.0; // tolerance in intercept approximation (pixels)
267 double _curveTol = 1.5; // tolerance of on-curve check (pixels)
268 mixin BaseRht;
269
270 // does coarsening to the multiple of tolerance
271 auto curveKey(Curve c)
272 {
273 import std.math;
274
275 auto angle = atan(c.m) / PI * 180;
276 return Key(rint(angle / _angleTol), rint(c.b / _interceptTol));
277 }
278
279 auto fitCurve(Range, SliceKind kind, Sample)(Slice!(kind, [2], Range) image, Sample sample)
280 {
281 int x1 = sample[0].x, y1 = sample[0].y;
282 int x2 = sample[1].x, y2 = sample[1].y;
283 if (x1 == x2)
284 return Curve(double.infinity, x1);
285 else
286 {
287 auto m = cast(double)(y1 - y2) / (x1 - x2);
288 auto b = y1 - m * x1;
289 return Curve(m, b);
290 }
291 }
292
293 bool onCurve(Curve curve, Point p)
294 {
295 import std.math;
296
297 int x = p.x, y = p.y;
298 double m = curve.m, b = curve.b;
299 if (m == double.infinity && fabs(x - b) < _curveTol)
300 return true;
301 else if (fabs(x * m + b - y) < _curveTol)
302 return true;
303 return false;
304 }
305 }
306
307 struct RhtCircles
308 {
309 alias Key = Tuple!(int, int, int);
310 alias Curve = Tuple!(double, "x", double, "y", double, "r");
311 enum sampleSize = 3;
312 double _centerTol = 5.0; // tolerance of center location (pixels)
313 double _radiusTol = 5.0; // tolerance of radius approfixmation (pixels)
314 double _curveTol = 8; // tolerance of on-curve check (proportional to radius)
315 mixin BaseRht;
316
317 // does coarsening to the multiple of tolerance
318 auto curveKey(Curve c)
319 {
320 return Key(toInt(c.x / _centerTol), toInt(c.y / _centerTol), toInt(c.r / _radiusTol));
321 }
322
323 auto fitCurve(Range, SliceKind kind, Sample)(Slice!(kind, [2], Range) image, Sample sample)
324 {
325 import std.math : sqrt, pow;
326
327 double x1 = sample[0].x, y1 = sample[0].y;
328 double x2 = sample[1].x, y2 = sample[1].y;
329 double x3 = sample[2].x, y3 = sample[2].y;
330 double ynorm = (y2 * y2 - y3 * y3) * (x2 - x1) / 2 + (x3 - x2) * (y2 * y2 - y1 * y1) / 2 + (
331 x1 - x3) * (y2 - y3) * (x2 - x1);
332 double y = ynorm / ((y2 - y3) * (x2 - x1) + (x3 - x2) * (y2 - y1));
333 double x = (y1 - y2) * (2 * y - y1 - y2) / (x2 - x1) / 2 + (x1 + x2) / 2;
334 double r = sqrt(pow(x - x1, 2) + pow(y - y1, 2));
335 return Curve(x, y, r);
336 }
337
338 bool onCurve(Curve curve, Point p)
339 {
340 import std.math : fabs;
341
342 double x = p.x, y = p.y;
343 double x0 = curve.x, y0 = curve.y, R = curve.r;
344 x -= x0;
345 y -= y0;
346 if (fabs(x * x + y * y - R * R) < _curveTol * _curveTol)
347 return true;
348 return false;
349 }
350 }
351
352 struct RhtEllipses
353 {
354 alias Key = Tuple!(int, int, int, int, int);
355 alias Curve = Tuple!(double, "x", double, "y", double, "a", double, "b", double, "phi");
356 enum sampleSize = 3;
357 double _centerTol = 2.0; // tolerance of center location (pixels)
358 double _axisTol = 2.0; // tolerance of axis approximation (pixels)
359 double _phiTol = 5.0; // tolerance of angle approximation (degrees)
360 double _curveTol = 0.05; // tolerance of on-curve check (proportional to axis)
361 mixin BaseRht;
362
363 // does coarsening to the multiple of tolerance
364 auto curveKey(Curve c)
365 {
366 import std.math;
367
368 return Key(toInt(c.x / _centerTol), toInt(c.y / _centerTol), toInt(c.a / _axisTol),
369 toInt(c.b / _axisTol), toInt(c.phi / _phiTol));
370 }
371
372 auto fitCurve(Range, SliceKind kind, Sample)(Slice!(kind, [2], Range) image, Sample sample)
373 {
374 Curve c;
375 //TODO:
376 return c;
377 }
378
379 bool onCurve(Curve curve, Point p)
380 {
381 import std.math;
382
383 double x = p.x, y = p.y;
384 double x0 = curve.x, y0 = curve.y, a = curve.a, b = curve.b;
385 x -= x0;
386 y -= y0;
387 double th = curve.phi * PI / 180;
388 double nx = x * cos(th) + y * sin(th);
389 double ny = -x * sin(th) + y * cos(th);
390 if (fabs(nx * nx / (a * a) + ny * ny / (b * b) - 1.0) < _curveTol)
391 return true;
392 return false;
393 }
394 }