Thread Subject: transpose of interpn() operations

The operation

  VI = INTERPN(X1,X2,X3,...,V,Y1,Y2,Y3,...)

is a linear operation on the data V, meaning that there
exists a matrix A such that

VI(:)=A*V(:)


In many applications, however, we need to be able to
implement multiplication by the transpose

VII=A'*VI

Since A is typically gigantic, you do not of course want to
compute A or A.' directly.

Does anyone know if this transpose operation is implemented
somewhere in a more or less optimized way?

I know that one can do this with accumarray(SUBS,VAL),
however, you then need to derive the SUBS data from the
(Y1,Y2,Y3,...) coordinates which is a time-consuming and
memory-consuming computation when implemented directly in
MATLAB.

I have the vague recollection that this was available
somewhere in the library of image registration functions,
but can't remember where. In any case, I would prefer to be
able to do it without the image processing toolbox.

"Matt " <mjacobson.removethis@xorantech.com> wrote in
message <g6vt59$k9o$1@fred.mathworks.com>...

>
> Does anyone know if this transpose operation is implemented
> somewhere in a more or less optimized way?
>
> I know that one can do this with accumarray(SUBS,VAL),
> however, you then need to derive the SUBS data from the
> (Y1,Y2,Y3,...) coordinates which is a time-consuming and
> memory-consuming computation when implemented directly in
> MATLAB.
>

Are you interested specifically in tensorial linear
interpolaion (default) or a general interpolation, including
spline?

Bruno

"Matt " <mjacobson.removethis@xorantech.com> wrote in message
<g6vt59$k9o$1@fred.mathworks.com>...
>
> The operation
>
> VI = INTERPN(X1,X2,X3,...,V,Y1,Y2,Y3,...)
>
> is a linear operation on the data V, meaning that there
> exists a matrix A such that
>
> VI(:)=A*V(:)
>
>
> In many applications, however, we need to be able to
> implement multiplication by the transpose
>
> VII=A'*VI
>
> Since A is typically gigantic, you do not of course want to
> compute A or A.' directly.

Why not? Learn what a sparse matrix is,
and you can write that matrix directly,
then performing any such operations
with ease.

John

> Are you interested specifically in tensorial linear
> interpolaion (default) or a general interpolation,
including
> spline?
 
I'd settle for tensorial linear. I wasn't too confident of
even this being available.

But I'd be even happier, of course, if higher order
tensorial splines were available, too.

"Matt " <mjacobson.removethis@xorantech.com> wrote in
message <g76etp$1d5$1@fred.mathworks.com>...
> > Are you interested specifically in tensorial linear
> > interpolaion (default) or a general interpolation,
> including
> > spline?
>
> I'd settle for tensorial linear. I wasn't too confident of
> even this being available.

I don't know about availability, but you could of course
program it with not much effort. Using sparse matrix as John
suggest is also a convenient way to accomplish the task.

>
> But I'd be even happier, of course, if higher order
> tensorial splines were available, too.

However, sparse matrix trick wouldn't work here, the
interpolation matrix is dense (elements are however larger
on the diagonal band). It is possible to "transpose" the
spline algorithm to compute its adjoint. Again, I have no
idea if such tool exists.

If you need any help to program the later, just let me know.

Bruno

> Why not? Learn what a sparse matrix is,
> and you can write that matrix directly,
> then performing any such operations
> with ease.

Believe me, I'm no stranger to sparse matrices, but it
won't do the trick.

For the applications I work with, a typical size for V
might be 400x400x300.

Even if I stuck with tensorial linear interpolation, and
even if sparse matrices supported single precision floats,
the size of A would easily run to about 2GB and strain
MATLAB's memory limits. On top of that, I would have to
convert V to type double to perform the multiplication.

Aside from memory issues, the time required to generate A
is problematic. In typical applications, you may need to
repeat these operations with several sets of interpolation
coordinate data. This means you will have to regenerate A
multiple times.

Anyway, if I'm reaching for the moon, so be it. I'm just
surprised that there hasn't been more interest in a coded
routine for A.'. It seems to me like something that would
often be needed.
 

"Matt " <mjacobson.removethis@xorantech.com> wrote in
message <g76fpu$7lt$1@fred.mathworks.com>...

>
> For the applications I work with, a typical size for V
> might be 400x400x300.
>

And what is the typical size of your VI Matt?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in
message <g76gia$d07$1@fred.mathworks.com>...
> "Matt " <mjacobson.removethis@xorantech.com> wrote in
> message <g76fpu$7lt$1@fred.mathworks.com>...
>
> >
> > For the applications I work with, a typical size for V
> > might be 400x400x300.
> >
>
> And what is the typical size of your VI Matt?
>
> Bruno

In the cases I'm thinking of, you almost always have
size(VI) equal to or comparable to size(V).

Think, for example, of deforming a 3D image. You apply a
geometric transformation to the original coordinates of the
image samples, and then you need to interpolate/resample at
all of the new locations.

> If you need any help to program the later, just let me
know.
>
> Bruno

Thanks. I don't think it would be so bad to program, but it
would have been great if someone else had already done it
for me :)

Hi Matt, I am sorry to dig out this old message but have you found a solution for this 'Transpose' interpolation/transformation operation?

Currently, I use your Bspline Matrix Interpolation; however, I really need a transpose operation of it. (Just to test the program, the efficiency could be considered later).

Thanks a lot,
Aaronne.

"Aaronne" wrote in message <j4v7un$sbu$1@newscl01ah.mathworks.com>...
> Hi Matt, I am sorry to dig out this old message but have you found a solution for this 'Transpose' interpolation/transformation operation?
>
> Currently, I use your Bspline Matrix Interpolation; however, I really need a transpose operation of it. (Just to test the program, the efficiency could be considered later).
>

I should add that this is NOT a "transpose", but an
inverse operator. Transpose and inverse are only the
same when the operator is of a specific orthogonal
class.

Next, interpn and its cousins are only locally linear.
So while you may find an inverse, there may exist
MANY inverses for any point. Worse, since interpn
is an n --> 1 operator, there will be infinitely many
such inverse points, an n-1 manifold of solutions.

And finally, higher order spline interpolants will not
be as simply inverted.

John

"Aaronne" wrote in message <j4v7un$sbu$1@newscl01ah.mathworks.com>...
> Hi Matt, I am sorry to dig out this old message but have you found a solution for this 'Transpose' interpolation/transformation operation?
>
> Currently, I use your Bspline Matrix Interpolation; however, I really need a transpose operation of it. (Just to test the program, the efficiency could be considered later).
=========

I have a MEX routine that does it for 2D bilinear interpolation, but nothing more general. Also, it's proprietary. Sorry.

As I told you before, though, I don't think it would be too hard to modify any of the interpolation codes I've linked you to (e.g. the FEX) to do what you're looking for.

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <j4vg9e$scv$1@newscl01ah.mathworks.com>...
>
>
> I should add that this is NOT a "transpose", but an
> inverse operator.

No, we are definitely talking about a transpose. Again the statement of the problem is this: the following function of V

f(V) = INTERPN(X1,X2,X3,...,V,Y1,Y2,Y3,...)

is equivalent to a matrix-vector multiplication

z1=A*V(:)

We would like a routine, preferably a highly optimized mex, that computes

z2=A.'*z1

We are NOT trying to recover V(:).

"Matt J" wrote in message <j4vm4f$o46$1@newscl01ah.mathworks.com>...
>
> No, we are definitely talking about a transpose. Again the statement of the problem is this: the following function of V
>
> f(V) = INTERPN(X1,X2,X3,...,V,Y1,Y2,Y3,...)
>
> is equivalent to a matrix-vector multiplication
>
> z1=A*V(:)
>
> We would like a routine, preferably a highly optimized mex, that computes
>
> z2=A.'*z1
>
> We are NOT trying to recover V(:).
=====================

And, as per Aaronne's point, we would like to be able to do this for arbitrary interpolation kernels, e.g. cardinal B-splines of order 3, and not simply the kernels offered by INTERPN.