Multi-variable interpolation, grids and 2D guessing

Hi!

I have a regular grid of (x,y) points. Each point in turn has an associated value of a range and frequency variable, that means that each (xi,yi) point has an associated (ri, fi) value. Another way of seeing it is that for each (x,y) point I have a r grid and a f grid ( r = function(x,y) and f = function(x,y)). Hence the r and f are not 'regular' let's say, as they are constructed with a specific function.

Now my question is: I want to find, given a set of (r,f) points, which is the corresponding (x,y) point that would yield those (r,f) values. Very likely the corresponding (x,y) will not exist in the original grid so they have to be 'guessed' somehow (probably with a given tolerance?). So as a result, for the set of (r,f) points I would like to find the set of guessed (x,y) points.

I don't know whether I made myself clear...hopefully yes!!

Thanks everyone in advance!!

Answer 1

[John D'Errico]

Yes. you make yourself clear.

This is a nonlinear optimization problem. Sometimes it might be called a calibration problem.

Actually, it is essentially a rootfinding problem in 2-dimensions, as long as the transformation from the (x,y) domain to the (r,f) domain is nonlinear. If that mapping was linear, then the solution is simple and linear itself too. That is, if the transformation was a linear one, then the solution is no more difficult than a solution of a linear system of equations. If any nonlineararity exists at all, then you might use a tool like fsolve. However fsolve can sometimes fail to produce a valid result, even when one exists.

Be careful here. If you are using a tool like fsolve to solve the problem, but you treat the mapping forward from (x,y) into the (r,f) domain using interp2 with linear interpolation, then fsolve will have a problem, since interp2 need not produce a differentiable function.

There are several problems in either case. In some circumstances there will be no solution at all, even if the problem can be solved using a linear mapping. Or, if the mapping is sufficiently nonlinear, then there may be several solutions for any (r,f), all equally valid.

What is the best way to solve it? There is one fairly foolproof solution that is not too difficult to implement. However, if the transformation was actually a linear mapping, then there is an even easier solution. Some people might know the phrase "affine mapping" or "affine transformation". Either correspond to what I described as linear above.

https://mathworld.wolfram.com/AffineTransformation.html

So which is it, before I spend an hour to write something up sufficiently well that you will understand what I am doing? (And then respond to your questions, which people always have anyway.) Is the mapping linear or nonlinear?

EDIT: How to solve the NONLINEAR inverse mapping.

This really is just a short explanation of how to solve a nonlinear rootfinding problem in 2-dimensions. First, I'll return to a more clear explanation of the problem.

Thus, imagine you start out with a sheet of paper, laid out in a grid, thus (x,y). Now, imagine you taking the sheet, and stretch it nonlinearly, deforming it like a piece of rubber. We can visualize it as...

rfun = @(x,y) x + (y+1).^2/5;
ffun = @(x,y) (x+2).*sin((2*x+y)/10 - 0.5);

[xg,yg] = meshgrid(-1:.1:1);

rg = rfun(xg,yg);
fg = ffun(xg,yg);

plot(xg,yg,'-k',xg',yg','-k')
xlabel X
ylabel Y
title 'Undeformed lattice'

Next, look at the deformed sheet. Be careful. The first example I cooked up had that sheet twisted and deformed in VERY strange ways. It was just too nonlinear of a transformation. The example I show below is a VERY MILD transformation, but even here you should see it is becoming problematic in one corner.

plot(rg,fg,'-r',rg',fg','-r')
xlabel R
ylabel F
title 'Rubber sheet deformation'

As you can see, in the upper right corner, the rubber sheet has actually become just slightly twisted on itself.

How should we interpret these plots? Here, I want you to visualize a set of arrows connecting the two figures. Every point in the (x,y) lattice has a corresponding point in the (r,f) lattice. Your goal of course, is to solve the inverse problem. That is, for any point in the deformed sheet, you want to know what the origin of that point would have been, in the original lattice.

I hope we now understand the problem. Now, can we solve it? How would that be done? ... (Writing more)

[John D'Errico]

By the way, please don't add answers just to make a comment. Your comment moved here:

"Thanks very much for the detailed answer. Indeed I realize now that if the problem were linear the solution may be easier. I first thought it could be treated like it, but it is not. However, I would not like to go through the full equation to invert it somehow, because r = fun(x,y) is a relatively complex function, f = fun(x,y) is simpler, but they are both implemented in a set of .m function files, and I would not like to go through them trying to invert them...and in any case as you mention if it is not linear this is not straightforward.

I'm having a look at affline if I understand it, and fsolve to do some testing...

Thanks!!"

[John D'Errico]

No. Affine is NOT an option, since the problem is nonlinear. NOT SO. Period.

You can use tools from the optimization toolbox. In fact, lsqnonlin will be a better choice than fsolve here, as I think about it. I can also show how to solve this using contourc. Let me write up an example. It will take a few minutes to do so though. I'll put it as part of my answer as an edit, since if I write it in the comments, things can get lost in the flow of later comments.

[Albert]

Ok, I understand with your example that due to the non-linearity of the problem there may be parts of the problem where a solution is not really stable, while other parts behave more as a pseudo-linear problem, that's how I interpret it.

Answer 2

[Albert]

Thanks very much for the detailed answer. Indeed I realize now that if the problem were linear the solution may be easier. I first thought it could be treated like it, but it is not. However, I would not like to go through the full equation to invert it somehow, because r = fun(x,y) is a relatively complex function, f = fun(x,y) is simpler, but they are both implemented in a set of .m function files, and I would not like to go through them trying to invert them...and in any case as you mention if it is not linear this is not straightforward.

I'm having a look at affline if I understand it, and fsolve to do some testing...

Thanks!!

[John D'Errico]

No. Affine is NOT an option, since the problem is nonlinear. NOT SO. Period.

You can use tools from the optimization toolbox. In fact, lsqnonlin will be a better choice than fsolve here, as I think about it. I can also show how to solve this using contourc. Let me write up an example. It will take a few minutes to do so though. I'll put it as part of my answer as an edit, since if I write it in the comments, things can get lost in the flow of later comments.

[Albert]

I think something in these terms works fine!

    contswath = contourc(x,y,r),[ri ri]);
    iso_r = [contswath(1,2:end);contswath(2,2:end)].';
    
    contswath = contourc(x,y,f),[fi fi]);
    iso_f = [contswath(1,2:end);contswath(2,2:end)].';
    
    [xi,yi] = polyxpoly(iso_r(:,1),iso_r(:,2),iso_f(:,1),iso_f(:,2));

[Albert]

But as I mentioned earlier I would like to use a set of ri and fi input values and retrieve the corresponding (xi,yi) in one or a few calls, not having to call contourc a lot of times...

[John D'Errico]

WILL YOU PLEASE STOP ADDING ANSWERS WILLY NILLY? Learn to use comments!

Given a r value, there are many (x,y) points that fulfill it (there is an iso-contour in fact), same thing happens with the f value. Maybe I could try to find the intersection between r and f contour lines, but this is probably not very efficient anyway...

[Albert]

Yes, ok sorry I mixed answers and comments!!

[John D'Errico]

If you already know how to use contourc, then just use it. That was going to be one of the ways I show how to solve it.

NO. You cannot solve this problem in one call, in a vectorized fashion. However there is ABSOLUTELY NOTHING that stops you from using a loop. You can even write a function that does exactly that - take a set of (r,f) pairs, then uses a loop, solving those problems one at a time, but returns the set of all solutions.

This is the virtue of MATLAB, in that if you don't find the tool that does exactly what you want, just write your own function that works in exactly the way you want it to work.

Note that the polyxpoly solution may return multiple solutions.

Finally, note that IF the mapping were a linear one, then a vectorized solution is trivial to write. But it aprently is not so. So there is no vectorized solution easily written, unless you also can compute the inverse mapping analytically.

[Albert]

Thanks John, yes, in the end I used contourc. I understand looping is always an option, but sometimes I just go straight to that solution because it's more intuitive, but many times there is a more efficient solution to this.

Thanks!!