{"id":1742,"date":"2016-04-12T17:36:33","date_gmt":"2016-04-12T21:36:33","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=1742"},"modified":"2019-11-01T16:36:26","modified_gmt":"2019-11-01T20:36:26","slug":"intersecting-curves-that-dont-intersect","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2016\/04\/12\/intersecting-curves-that-dont-intersect\/","title":{"rendered":"Intersecting curves that don’t intersect"},"content":{"rendered":"

A post in MATLAB Answers<\/a> earlier this year reminded me that working on a discrete grid can really mess up apparently obvious notions about geometry.<\/p>

User Hg<\/a> offered an image containing two intersecting curves. The problem is that the intersecting curves didn't intersect!<\/p>

Here's the image:<\/p>

\"\" <\/p>

The red curve and the blue curve, which obviously cross each other, do not have any pixels in common. Hg's question: \"What can I do to estimate the intersection point?\"<\/p>

Before I get into answering the original question, let me take a brief side trip to see if I can reconstruct the original input matrix (or at least something close). The image posted on MATLAB Answers had the pixels obviously blown up, with different colors applied to distinguish between the curves. Can I get back to just a matrix with 0s, 1s, and 2s corresponding to pixels belonging to the background and the two curves?<\/p>

It helps that the image is a PNG file, which is losslessly compressed. I took a close look at the pixels using imtool<\/tt>, and I was able to determine that pixels belong to one curve were colored using [237 28 36], and pixels belonging to the second curve were colored using [0 162 232]. Also, each pixel seemed to be magnified by about a factor of 16. Let's run with those numbers.<\/p>

url = 'https:\/\/www.mathworks.com\/matlabcentral\/answers\/uploaded_files\/44249\/intersect2.png'<\/span>;\r\nrgb = imread(url);\r\nred = rgb(:,:,1);\r\ngreen = rgb(:,:,2);\r\nblue = rgb(:,:,3);\r\n<\/pre>
curve1 = (red == 237) & (green == 28) & (blue == 36);\r\ncurve2 = (red == 0) & (green == 162) & (blue == 232);\r\nL = zeros(size(curve1));\r\nL(curve1) = 1;\r\nL(curve2) = 2;\r\nimshow(L,[],'InitialMagnification'<\/span>,'fit'<\/span>)\r\n<\/pre>\"\" 
That looks good. Let's try shrinking it down by a factor of 16.<\/p>
L16 = imresize(L,1\/16,'nearest'<\/span>);\r\nimshow(L16,[],'InitialMagnification'<\/span>,'fit'<\/span>);\r\n<\/pre>\"\" 
Well, that's not perfect. But it's close enough to work with for answering Hg's question.<\/p>
Image Analyst<\/a> offered a potential solution:<\/p>
There is no one pixel where the overlap occurs. If you'll accept any of those 4 pixel locations as an overlap, then perhaps if you dilated, ANDed, then called bwulterode. Something like (untested)<\/i><\/p>
intImage = imdilate(bw1, true(3)) & imdilate(bw2, true(3));\r\nintPoints = bwulterode(intImage);<\/pre>
Image Analyst's idea is to use morphological dilation to thicken each curve individually, then use a logical AND operator (&<\/tt>) to determine the set of pixels that belong to both thickened curves, and then compute the ultimate erosion using bwerode<\/tt> to shrink the intersection area down. Here it is in several steps.<\/p>
curve1 = (L16 == 1);\r\ncurve2 = (L16 == 2);\r\ncurve1_thickened = imdilate(curve1,ones(3,3));\r\nimshow(curve1_thickened,'InitialMagnification'<\/span>,'fit'<\/span>)\r\ntitle('Curve 1 (thickened)'<\/span>)\r\n<\/pre>\"\" 
curve2_thickened = imdilate(curve2,ones(3,3));\r\nimshow(curve2_thickened,'InitialMagnification'<\/span>,'fit'<\/span>)\r\ntitle('Curve 2 (thickened)'<\/span>)\r\n<\/pre>\"\" 
curve_intersection = curve1_thickened & curve2_thickened;\r\nimshow(curve_intersection,'InitialMagnification'<\/span>,'fit'<\/span>)\r\ntitle('Intersection of thickened curves'<\/span>)\r\n<\/pre>\"\" 
ultimate_erosion = bwulterode(curve_intersection);\r\nimshow(ultimate_erosion,'InitialMagnification'<\/span>,'fit'<\/span>)\r\ntitle('Ultimate erosion of intersection'<\/span>)\r\n<\/pre>\"\" 
Let's overlay the output of the ultimate erosion on top of the original image.<\/p>
imshow(imfuse(L16,ultimate_erosion,'falsecolor'<\/span>),'InitialMagnification'<\/span>,'fit'<\/span>)\r\n<\/pre>\"\" 
You can see that the ultimate erosion is not a single point. An alternative is to compute the centroids of the intersection \"blobs.\"<\/p>
s = regionprops(ultimate_erosion,'Centroid'<\/span>)\r\n<\/pre>
\r\ns = \r\n\r\n    Centroid: [9.5000 6.5000]\r\n\r\n<\/pre>
There's just one intersection blob, and its centroid is at (9.5,6.5). Let's plot that.<\/p>
imshow(L16,[],'InitialMagnification'<\/span>,'fit'<\/span>)\r\nhold on<\/span>\r\nplot(s.Centroid(1),s.Centroid(2),'go'<\/span>,'MarkerSize'<\/span>,15,...<\/span>\r\n    'MarkerFaceColor'<\/span>,'g'<\/span>)\r\nhold off<\/span>\r\n<\/pre>\"\" 
How would you approach this problem? Leave your ideas and thoughts in the comments.<\/p>