interpolating or fitting spline issues

Question

I'm having some trouble with interpolating and/or fitting splines. There are two aspects to my issue.

I am selecting points on an image manually, and want to interpolate the points between. The points I am clicking on are fixed, and any interpolation MUST go through these points. There is not a fixed number of points. One image may have 5, another may have 20 (although I keep to a minimum of 5). The points are not evenly distributed alone the line that needs interpolating. In fact, most of the time there are large gaps to be filled. All points are recorded as x and y coordinate values using:

[x,y] = ginput; 

Each (full) spline needs to be of length szG (normally 100). So I calculate the x and y interpolants, and the yy is reversed as its to be overlaid on the image and plot/image indexes are reversed:

xx = min(x):(max(x)-min(x)-1)/(szG-1):max(x);
yy = max(y):-(max(y)-min(y)-1)/(szG-1):min(y);

I then calculate the curve with:

newX = interp1(y,x,yy, 'cubic');
newY = interp1(x,y,xx, 'cubic');

However, when I plot these, the interpolated line doesn't go through the initial points. It's actually quite far away a lot of the time.

The second problem is that often two selected points have the same x or y value. This means that I get an error when it comes to interpolation as the values need to be distinct.

How can you ensure the interpolated line goes through selected data points? Also, how can the distinct value issue be resolved?

Answer 1

One problem that you have is that you aren't reversing the image coordinates properly. The origin of your axis needs to be at the top left, rather than the bottom left and so simply negating the y axis doesn't do the job for you. You need to make sure that the origin is at the top left, and so you need to take your y coordinates and subtract by the total number of rows to facilitate this. However, a trick to reversing the coordinates without you having to do it manually is to do it on the figure itself. This way, you don't have to worry about doing any calculations to reversing the y axis. Simply do this once the figure is open:

axis ij;

To revert back to the original system where the origin is at the bottom left, simply do:

axis xy;

However, if you are using imshow to display your images, then there's absolutely no need to reverse the y coordinates because this is already in effect when you show the image (a.k.a axis ij is called under the hood).

If you do want to reverse the coordinates manually yourself, simply do this assuming your image is stored in im:

y = size(im, 1) - y;

In addition, you're not using interp1 correctly. You specify control points with x and y, and you use xx as the new input coordinates. The output will give you new y coordinates that interpolate along this line, with xx being the query points to be interpolated. Therefore, you can just use the output of interp1 directly for your result and you only need to call it once.

Here's a reproducible example. Assuming you have the image processing toolbox:

im = imread('cameraman.tif');
imshow(im);
[x,y] = ginput(5);

Here, I'm choosing 5 points. Note that I don't need to reverse the y axis because imshow already does that for you. I've chosen these 5 coordinates:

>> x

x =

   23.0000
   71.0000
  143.0000
  200.0000
  238.0000

>> y

y =

    75
    43
    27
    47
    77

Now let's make the spline out of these points:

szG = 100;
xx = min(x):(max(x)-min(x)-1)/(szG-1):max(x);
yy = interp1(x, y, xx, 'cubic');

If I can make a minor point, you can do the same thing for your xx coordinates by using linspace by:

xx = linspace(min(x), max(x), szG);

Now let's show these points where the keypoints are in blue with large circles and red is the interpolated line:

hold on;
plot(x, y, 'b.', 'MarkerSize', 16);
plot(xx, yy, 'r');

This is what I get:

I get the points to go through the control points, and it's probably because you're not using interp1 properly. Now, to fix the issue where you have duplicate points. This isn't useful at all for interp1, so one thing I suggest is to filter out any 2D points that are not unique and you can do this by using the unique function. Specifically, place the x and y coordinates into a 2D array and filter out the points based on the rows. This way, you'll eliminate any non-unique x and y pairs for use in interp1. Something like this:

out = unique([x(:) y(:)], 'rows', 'stable');
x = out(:,1);
y = out(:,2);

unique also has a side-effect where the points are arranged in sorted order, such that the first column is used as a sorting key. You probably don't want this and you want to maintain the same order of points as you had when you clicked on the figure. Therefore, you need to use the 'stable' flag to ensure this doesn't happen. After this operation, x and y will now contain distinct coordinates for use in interp1.

---

Hope this helps. Good luck!