1 | /*M/////////////////////////////////////////////////////////////////////////////////////// |
2 | // |
3 | // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. |
4 | // |
5 | // By downloading, copying, installing or using the software you agree to this license. |
6 | // If you do not agree to this license, do not download, install, |
7 | // copy or use the software. |
8 | // |
9 | // |
10 | // License Agreement |
11 | // For Open Source Computer Vision Library |
12 | // |
13 | // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. |
14 | // Copyright (C) 2008-2011, Willow Garage Inc., all rights reserved. |
15 | // Third party copyrights are property of their respective owners. |
16 | // |
17 | // @Authors |
18 | // Nghia Ho, nghiaho12@yahoo.com |
19 | // |
20 | // Redistribution and use in source and binary forms, with or without modification, |
21 | // are permitted provided that the following conditions are met: |
22 | // |
23 | // * Redistribution's of source code must retain the above copyright notice, |
24 | // this list of conditions and the following disclaimer. |
25 | // |
26 | // * Redistribution's in binary form must reproduce the above copyright notice, |
27 | // this list of conditions and the following disclaimer in the documentation |
28 | // and/or other materials provided with the distribution. |
29 | // |
30 | // * The name of OpenCV Foundation may not be used to endorse or promote products |
31 | // derived from this software without specific prior written permission. |
32 | // |
33 | // This software is provided by the copyright holders and contributors "as is" and |
34 | // any express or implied warranties, including, but not limited to, the implied |
35 | // warranties of merchantability and fitness for a particular purpose are disclaimed. |
36 | // In no event shall the OpenCV Foundation or contributors be liable for any direct, |
37 | // indirect, incidental, special, exemplary, or consequential damages |
38 | // (including, but not limited to, procurement of substitute goods or services; |
39 | // loss of use, data, or profits; or business interruption) however caused |
40 | // and on any theory of liability, whether in contract, strict liability, |
41 | // or tort (including negligence or otherwise) arising in any way out of |
42 | // the use of this software, even if advised of the possibility of such damage. |
43 | // |
44 | //M*/ |
45 | #include "precomp.hpp" |
46 | |
47 | namespace cv |
48 | { |
49 | |
50 | static inline bool _isOnPositiveSide(const Point2f& line_vec, const Point2f& line_pt, const Point2f& pt) |
51 | { |
52 | //we are interested by the cross product between the line vector (line_vec) and the line-to-pt vector (pt-line_pt) |
53 | //the sign of the only non-null component of the result determining which side of the line 'pt' is on |
54 | //the "positive" side meaning depends on the context usage of the current function and how line_vec and line_pt were filled |
55 | return (line_vec.y*(line_pt.x-pt.x) >= line_vec.x*(line_pt.y-pt.y)); |
56 | } |
57 | |
58 | static int _rotatedRectangleIntersection( const RotatedRect& rect1, const RotatedRect& rect2, std::vector<Point2f> &intersection ) |
59 | { |
60 | CV_INSTRUMENT_REGION(); |
61 | |
62 | Point2f vec1[4], vec2[4]; |
63 | Point2f pts1[4], pts2[4]; |
64 | |
65 | rect1.points(pts: pts1); |
66 | rect2.points(pts: pts2); |
67 | |
68 | // L2 metric |
69 | float samePointEps = 1e-6f * (float)std::max(a: rect1.size.area(), b: rect2.size.area()); |
70 | |
71 | int ret = INTERSECT_FULL; |
72 | |
73 | // Specical case of rect1 == rect2 |
74 | { |
75 | bool same = true; |
76 | |
77 | for( int i = 0; i < 4; i++ ) |
78 | { |
79 | if( fabs(x: pts1[i].x - pts2[i].x) > samePointEps || (fabs(x: pts1[i].y - pts2[i].y) > samePointEps) ) |
80 | { |
81 | same = false; |
82 | break; |
83 | } |
84 | } |
85 | |
86 | if(same) |
87 | { |
88 | intersection.resize(new_size: 4); |
89 | |
90 | for( int i = 0; i < 4; i++ ) |
91 | { |
92 | intersection[i] = pts1[i]; |
93 | } |
94 | |
95 | return INTERSECT_FULL; |
96 | } |
97 | } |
98 | |
99 | // Line vector |
100 | // A line from p1 to p2 is: p1 + (p2-p1)*t, t=[0,1] |
101 | for( int i = 0; i < 4; i++ ) |
102 | { |
103 | vec1[i].x = pts1[(i+1)%4].x - pts1[i].x; |
104 | vec1[i].y = pts1[(i+1)%4].y - pts1[i].y; |
105 | |
106 | vec2[i].x = pts2[(i+1)%4].x - pts2[i].x; |
107 | vec2[i].y = pts2[(i+1)%4].y - pts2[i].y; |
108 | } |
109 | |
110 | //we adapt the epsilon to the smallest dimension of the rects |
111 | for( int i = 0; i < 4; i++ ) |
112 | { |
113 | samePointEps = std::min(a: samePointEps, b: std::sqrt(x: vec1[i].x*vec1[i].x+vec1[i].y*vec1[i].y)); |
114 | samePointEps = std::min(a: samePointEps, b: std::sqrt(x: vec2[i].x*vec2[i].x+vec2[i].y*vec2[i].y)); |
115 | } |
116 | samePointEps = std::max(a: 1e-16f, b: samePointEps); |
117 | |
118 | // Line test - test all line combos for intersection |
119 | for( int i = 0; i < 4; i++ ) |
120 | { |
121 | for( int j = 0; j < 4; j++ ) |
122 | { |
123 | // Solve for 2x2 Ax=b |
124 | const float x21 = pts2[j].x - pts1[i].x; |
125 | const float y21 = pts2[j].y - pts1[i].y; |
126 | |
127 | float vx1 = vec1[i].x; |
128 | float vy1 = vec1[i].y; |
129 | |
130 | float vx2 = vec2[j].x; |
131 | float vy2 = vec2[j].y; |
132 | |
133 | const float det = vx2*vy1 - vx1*vy2; |
134 | if (std::abs(x: det) < 1e-12)//we consider accuracy around 1e-6, i.e. 1e-12 when squared |
135 | continue; |
136 | const float detInvScaled = 1.f/det; |
137 | |
138 | const float t1 = (vx2*y21 - vy2*x21)*detInvScaled; |
139 | const float t2 = (vx1*y21 - vy1*x21)*detInvScaled; |
140 | |
141 | // This takes care of parallel lines |
142 | if( cvIsInf(value: t1) || cvIsInf(value: t2) || cvIsNaN(value: t1) || cvIsNaN(value: t2) ) |
143 | { |
144 | continue; |
145 | } |
146 | |
147 | if( t1 >= 0.0f && t1 <= 1.0f && t2 >= 0.0f && t2 <= 1.0f ) |
148 | { |
149 | const float xi = pts1[i].x + vec1[i].x*t1; |
150 | const float yi = pts1[i].y + vec1[i].y*t1; |
151 | |
152 | intersection.push_back(x: Point2f(xi,yi)); |
153 | } |
154 | } |
155 | } |
156 | |
157 | if( !intersection.empty() ) |
158 | { |
159 | ret = INTERSECT_PARTIAL; |
160 | } |
161 | |
162 | // Check for vertices from rect1 inside recct2 |
163 | for( int i = 0; i < 4; i++ ) |
164 | { |
165 | // We do a sign test to see which side the point lies. |
166 | // If the point all lie on the same sign for all 4 sides of the rect, |
167 | // then there's an intersection |
168 | int posSign = 0; |
169 | int negSign = 0; |
170 | |
171 | const Point2f& pt = pts1[i]; |
172 | |
173 | for( int j = 0; j < 4; j++ ) |
174 | { |
175 | // line equation: Ax + By + C = 0 where |
176 | // A = -vec2[j].y ; B = vec2[j].x ; C = -(A * pts2[j].x + B * pts2[j].y) |
177 | // check which side of the line this point is at |
178 | // A*x + B*y + C <> 0 |
179 | // + computation reordered for better numerical stability |
180 | |
181 | const bool isPositive = _isOnPositiveSide(line_vec: vec2[j], line_pt: pts2[j], pt); |
182 | |
183 | if( isPositive ) |
184 | { |
185 | posSign++; |
186 | } |
187 | else |
188 | { |
189 | negSign++; |
190 | } |
191 | } |
192 | |
193 | if( posSign == 4 || negSign == 4 ) |
194 | { |
195 | intersection.push_back(x: pts1[i]); |
196 | } |
197 | } |
198 | |
199 | // Reverse the check - check for vertices from rect2 inside recct1 |
200 | for( int i = 0; i < 4; i++ ) |
201 | { |
202 | // We do a sign test to see which side the point lies. |
203 | // If the point all lie on the same sign for all 4 sides of the rect, |
204 | // then there's an intersection |
205 | int posSign = 0; |
206 | int negSign = 0; |
207 | |
208 | const Point2f& pt = pts2[i]; |
209 | |
210 | for( int j = 0; j < 4; j++ ) |
211 | { |
212 | // line equation: Ax + By + C = 0 where |
213 | // A = -vec1[j].y ; B = vec1[j].x ; C = -(A * pts1[j].x + B * pts1[j].y) |
214 | // check which side of the line this point is at |
215 | // A*x + B*y + C <> 0 |
216 | // + computation reordered for better numerical stability |
217 | |
218 | const bool isPositive = _isOnPositiveSide(line_vec: vec1[j], line_pt: pts1[j], pt); |
219 | |
220 | if( isPositive ) |
221 | { |
222 | posSign++; |
223 | } |
224 | else |
225 | { |
226 | negSign++; |
227 | } |
228 | } |
229 | |
230 | if( posSign == 4 || negSign == 4 ) |
231 | { |
232 | intersection.push_back(x: pts2[i]); |
233 | } |
234 | } |
235 | |
236 | int N = (int)intersection.size(); |
237 | if (N == 0) |
238 | { |
239 | return INTERSECT_NONE; |
240 | } |
241 | |
242 | // Get rid of duplicated points |
243 | const int Nstride = N; |
244 | cv::AutoBuffer<float, 100> distPt(N * N); |
245 | cv::AutoBuffer<int> ptDistRemap(N); |
246 | for (int i = 0; i < N; ++i) |
247 | { |
248 | const Point2f pt0 = intersection[i]; |
249 | ptDistRemap[i] = i; |
250 | for (int j = i + 1; j < N; ) |
251 | { |
252 | const Point2f pt1 = intersection[j]; |
253 | const float d2 = normL2Sqr<float>(pt: pt1 - pt0); |
254 | if(d2 <= samePointEps) |
255 | { |
256 | if (j < N - 1) |
257 | intersection[j] = intersection[N - 1]; |
258 | N--; |
259 | continue; |
260 | } |
261 | distPt[i*Nstride + j] = d2; |
262 | ++j; |
263 | } |
264 | } |
265 | while (N > 8) // we still have duplicate points after samePointEps threshold (eliminate closest points) |
266 | { |
267 | int minI = 0; |
268 | int minJ = 1; |
269 | float minD = distPt[1]; |
270 | for (int i = 0; i < N - 1; ++i) |
271 | { |
272 | const float* pDist = distPt.data() + Nstride * ptDistRemap[i]; |
273 | for (int j = i + 1; j < N; ++j) |
274 | { |
275 | const float d = pDist[ptDistRemap[j]]; |
276 | if (d < minD) |
277 | { |
278 | minD = d; |
279 | minI = i; |
280 | minJ = j; |
281 | } |
282 | } |
283 | } |
284 | CV_Assert(fabs(normL2Sqr<float>(intersection[minI] - intersection[minJ]) - minD) < 1e-6); // ptDistRemap is not corrupted |
285 | // drop minJ point |
286 | if (minJ < N - 1) |
287 | { |
288 | intersection[minJ] = intersection[N - 1]; |
289 | ptDistRemap[minJ] = ptDistRemap[N - 1]; |
290 | } |
291 | N--; |
292 | } |
293 | |
294 | // order points |
295 | for (int i = 0; i < N - 1; ++i) |
296 | { |
297 | Point2f diffI = intersection[i + 1] - intersection[i]; |
298 | for (int j = i + 2; j < N; ++j) |
299 | { |
300 | Point2f diffJ = intersection[j] - intersection[i]; |
301 | if (diffI.cross(pt: diffJ) < 0) |
302 | { |
303 | std::swap(a&: intersection[i + 1], b&: intersection[j]); |
304 | diffI = diffJ; |
305 | } |
306 | } |
307 | } |
308 | |
309 | intersection.resize(new_size: N); |
310 | |
311 | return ret; |
312 | } |
313 | |
314 | int rotatedRectangleIntersection( const RotatedRect& rect1, const RotatedRect& rect2, OutputArray intersectingRegion ) |
315 | { |
316 | CV_INSTRUMENT_REGION(); |
317 | |
318 | if (rect1.size.empty() || rect2.size.empty()) |
319 | { |
320 | intersectingRegion.release(); |
321 | return INTERSECT_NONE; |
322 | } |
323 | |
324 | // Shift rectangles closer to origin (0, 0) to improve the calculation of the intesection region |
325 | // To do that, the average center of the rectangles is moved to the origin |
326 | const Point2f averageCenter = (rect1.center + rect2.center) / 2.0f; |
327 | |
328 | RotatedRect shiftedRect1(rect1); |
329 | RotatedRect shiftedRect2(rect2); |
330 | |
331 | // Move rectangles closer to origin |
332 | shiftedRect1.center -= averageCenter; |
333 | shiftedRect2.center -= averageCenter; |
334 | |
335 | std::vector <Point2f> intersection; intersection.reserve(n: 24); |
336 | |
337 | const int ret = _rotatedRectangleIntersection(rect1: shiftedRect1, rect2: shiftedRect2, intersection); |
338 | |
339 | // If return is not None, the intersection Points are shifted back to the original position |
340 | // and copied to the interesectingRegion |
341 | if (ret != INTERSECT_NONE) |
342 | { |
343 | for (size_t i = 0; i < intersection.size(); ++i) |
344 | { |
345 | intersection[i] += averageCenter; |
346 | } |
347 | |
348 | Mat(intersection).copyTo(m: intersectingRegion); |
349 | } |
350 | else |
351 | { |
352 | intersectingRegion.release(); |
353 | } |
354 | |
355 | return ret; |
356 | } |
357 | |
358 | } // end namespace |
359 | |