Find 8 extreme pixels to determine skewness

Hello,

I'm currently taking a course in Computer Vision and I'm trying to solve a problem, which thus far seems impossible.

I've taken a photo of a map, but the photo is slightly skewed. I want to draw timezones on the map, but they end up in wrong places because of this shear, thus I need to correct the position of timelines for this to work.

The original image has had small blobs removed with 'bwareaopen'. In the image 'Points I want to find' I've circled the exact pixels I want to find in red.

I'm able to draw a bounding box around the entire map using the answer on this page but this leads to the following result: Problematic image.

I want to achieve something like this. This box can be calculated with the X- and Y-positions of the pixels I've circled. If I can find the corners of the box I want immediately, that's fine as well.

Can anyone point me in the right direction to find the exact pixels I've marked?

Regards, Bas

Answer 1

[Image Analyst]

Use bwconvhull().

[Bas Dalenoord]

Aah, thanks for that! This seems like the right direction, but I have a small problem remaining:

This is the convhull-output: image.

I need to draw lines (of the timezones) in between the top and the bottom. If the image is perfectly straight in the photo, these lines are spread evenly apart. If the image is skewed however, like in my convhull-image, the spreading of these lines is not even.

In the convhull-images, the timezones from the left are wider than the timezones on the right... To determine with what factor the wideness of each timezone has to be calculated I need to know which side is shorter and which side is longer.

I hope everything is clear, can you help me any further?

[Image Analyst]

Time zones should not be straight. I'm not sure how you're plotting them. Do you have the mapping toolbox? Maybe it has some tools for that.

If you need to warp your image, you can use imwarp(). See this page: http://blogs.mathworks.com/steve/2006/08/04/spatial-transformations-defining-and-applying-custom-transforms/

[Bas Dalenoord]

I really appreciate you trying to help me out, thanks so far!

I should have stated that this is an initial simple implementation, I am allowed to just draw straight lines as timezones (I will add the special cases later). If you look at the top of this image you can see that all timezones start with a straight line, evenly spaced apart. I need to be able to draw those lines and continue them all the way to the bottom.

If you take a look at this image, you see that the timezones to the right are actually wider than they were in the original image.

And for clarification: I want to just be able to draw lines like I did in this image. This is the synthetic image (which I printed an photographed in the original question) loaded into matlab, which I detected by taking only the topmost, rightmost, bottommost and leftmost pixels and drawing a box around that. That approach doesnt work with skewed image though.

If only I could find the actual corners of my convhull-image (the corners of the red rectangle in this image).

[Image Analyst]

Who says that those pins in the corners were necessarily pushed into the map in the 100% accurate place? The bounding quadrilateral could vary depending on the angle the pins are stuck in the paper.

[Bas Dalenoord]

The corners are magnets, their angle is thus fixed. There is some problem with shadows, but I can leave that out of the equation. It is not mission-critical to get it spot on, but I've now got the UTC-line half way across Europe, which is kinda wrong. It's at least got to be around England.

I'm pretty sure it must be possible, but am having a hard time even thinking it out at the moment. I can visualize what I want (which is in first case this), but I cant seem to think of a way to achieve this, even when I don't think in Matlab but for instance in C++.

[Bas Dalenoord]

I've been trying a bit over the weekend but haven't found a solution. I think the issue has been over-complexed at this point, so let me take a step back to a point from where I think I can solve the rest of the problem myself.

---

I've got this convex hull image.

I want to be able to get the length of all sides of the convex hull. How can I do this?

A problem with the convex hull is the rounded corners, for they might not be evenly round on each corner which distorts the length value of each side. How can I make the corners square?

I want to have the length of the sides of the rectangle in this image finally, I can use these values in further calculations.

Could you help me any further with this issue?

[Image Analyst]

Since the pins are not necessarily pushed in there accurately, and you don't have a perfect angle and aren't using a telecentric lens, then how about if you just use the bounding box? A bounding quadrilateral is not unique - there could be several that would work, but a bounding box is unique.

[Bas Dalenoord]

The bounding box does not give me information about the angle, which the individual sides of the convex hull do. This information is vital in calculating what the distance between the timezones should be, as there is a linear relation between the angle of the picture and the distance between the timezones.

I cannot see a way to do this with the bounding box only, which is actually what I tried before I added the corner markers. The bounding box of the map is quite easy to find, but this gives me no information about the angle as I said previously.

---

I suppose it is possible in Matlab to get the length of various sides of an object, even if these sides are at an angle (i.e. length of the sides of my somewhat diamond-shaped hull)?

Only influence in this information is the rounded corners which aren't evenly rounded around each corner, thus I want to compensate for that by making these corners square so the length of the side is more representative. I suppose that's possible as well, or am I wrong in this?

[Image Analyst]

See minboundrect() or minboundquad() in John D'Errico's suite: http://www.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects

[Bas Dalenoord]

Ah, I'd found that one before but was scared by the following warning:

Caveat - if you have only an image, don't expect these tools to work directly. (The image processing toolbox can do much here anyway - look there first.) In order to use these tools, I require a set of extracted points.

I haven't got a set of extracted points, for if I had i could calculate my own bounding quadrilateral. I'll look into the sources to see if I can find something useful though.

[Image Analyst]

Well you'd have to use find() to get the coordinates of the binary image so that his functions could then work. But be careful, find() returns (rows, columns), i.e. (y,x), not (x,y).

[Bas Dalenoord]

Hm. `minboundrect` gives me the bounding rectangle just like I was able to do previously. This does not give me info about the skewness.

`minboundquad` tells me the following:

Index exceeds matrix dimensions

I'm starting with this image, using the following code for bounding rectangle:

 [y, x] = find(imgWithoutSmallBlobs);
 [rx, ry, area] = minboundrect(x,y);

This gives me just the bounding box

The same call but with 'minboundquad' gives me the following error:

Index exceeds matrix dimensions.

Error in minboundquad (line 133)
k = (edgeangles(edgelist(:,4)) - edgeangles(edgelist(:,1)) <= pi);

'y' and 'x' are both '315271x1 double arrays', minboundquad takes quite a while to run as well, which is quite problematic.

Am I still doing something wrong? I've flipped the output of the find-command, so that shouldnt the problem...

[Bas Dalenoord]

I've been struggling with this issue for quite some time now, haven't really found a solution with convex hulls thusfar.

I ditched what I had before this afternoon and came up with the following brute-force solution:

• Call regionprops on the image to calculate circularity and the bounding box per blob
• Loop over all blobs, removing those that aren't circular enough
• Take the remaining blobs and their respective bounding boxes. The various corners of these boxes are the corners of my ROI.

This works quite well for specific images where the map is not rotated too much, it might not work for situations where the acquisition of the image was less optimal. I'll try and improve this a bit, for it is an actual solution to my problem.

The code I've written is the following:

% Regionprops
measurements = regionprops(imgWithoutSmallBlobs, 'Area', 'Perimeter', 'PixelIdxList', 'BoundingBox');

% Calculate circularity, only circular blobs will remain
allAreas = [measurements.Area];
allPerims = [measurements.Perimeter];
circularities = allPerims .^ 2 ./ (4*pi*allAreas);

% Remove non-circular blobs
numberOfBlobs = size(circularities,2);
removed = [];
for blobIdx = 1:numberOfBlobs
    circularity = circularities(blobIdx);
    measurement = measurements(blobIdx);

      if((circularity > 2) == 1)
          imgWithoutSmallBlobs( vertcat( measurement.PixelIdxList ) ) = 0; 
          removed = vertcat(removed,blobIdx);
      end
  end
  measurements(removed) = [];

numberOfBlobs = size(measurements,1);
if numberOfBlobs ~= 4
    numberOfBlobs
    error('Invalid number of blobs');
end

% Calculate the corners of the overall bounding box
topLeft = measurements(1).BoundingBox;
topLeftX = topLeft(1);
topLeftY = topLeft(2);

topRight = measurements(3).BoundingBox;
topRightX = topRight(1) + topRight(3);
topRightY = topRight(2);

bottomRight = measurements(4).BoundingBox;
bottomRightX = bottomRight(1) + bottomRight(3);
bottomRightY = bottomRight(2) + bottomRight(4);

bottomLeft = measurements(2).BoundingBox;
bottomLeftX = bottomLeft(1);
bottomLeftY = bottomLeft(2) + bottomLeft(4);

Dirty, but at least usable for this specific image... I'll try and improve it for different images later