I'm trying to write an efficient bilinear (2D)-interpolation, after reading some recipes, as a fortran-mex for Matlab that is used extensively throughout a long algorithm of solar image processing, and therefore is one of my main bottlenecks. I'm not a Fortran expert and this is one of my self-made exercise to understand and enforce efficient programming with Matlab. I managed to write a working code, partially vectorized, with one DO loop, but my poor skills makes me sure that there must be a faster and more elegant way of doing it. I was hoping, if faster, that one could avoid a DO LOOP to iterate over each given points where i need the interpolation.
Any advice that can speed up the code is most welcome.
I'm using gcc version 4.3.6.
I cannot use a later one as it the only one fully compatible with my version of Matlab. Any later version prevented the use of some embedded matlab-mex functions necessary for interfacing / displaying error messages in the matlab prompts.
My mex-setup can compile Fortran 90 and 95 (so far...)
I already have a C++ mex equivalent used so far and was hoping that a well-written Fortran code would be faster.
I'm new to this forum, I understood from the sticky notes in this forum that this was the right place. I apologize if I misread the notes and please re-direct me to any better place for this request.
The equations for bilinear interpolation that i'm using are :
f1 = wx1*Image(y1,x1) + wx2*Image(y1,x2) f2 = wx1*Image(y2,x1) + wx2*Image(y2,x2) IntIm= wy1*f1 + wy2*f2
x1 and x2 are the x-1 and x+1, same for y, x and y being the point at which the interpolation is needed. These coordinates are input as vectors, i.e, vector of x, 1D, and vector of y, 1D (see below for my attempt at implementing it)
wx* and wy* are the weights, which goes to NaN whenever x1 = x2 and y1 = y2
The Attempt at a Solution
SUBROUTINE L2DINTERPOL3(IntIm,Image,x,y,NPts,M,N) INTEGER*8 :: NPts,M,N,jj INTEGER*8 :: x1(NPts),y1(NPts),x2(NPts),y2(NPts) DOUBLE PRECISION :: IntIm(NPts) DOUBLE PRECISION :: Image(M,N) DOUBLE PRECISION :: xjj,yjj,f1,f2 DOUBLE PRECISION :: x(NPts),y(NPts) DOUBLE PRECISION :: wx1(NPts),wx2(NPts),wy1(NPts),wy2(NPts) c Take the nearest integers around the input. x1 = FLOOR(x) x2 = CEILING(x) y1 = FLOOR(y) y2 = CEILING(y) wx1 = (x2-x)/(x2-x1) wx2 = (x-x1)/(x2-x1) wy1 = (y2-y)/(y2-y1) wy2 = (y-y1)/(y2-y1) c Whenever the input are already integers, weights go to infinity, c so set each pair to 1 and 0. WHERE(x1 .EQ. x2) wx1 = 1. wx2 = 0. END WHERE WHERE (y1 .EQ. y2) wy1 = 1. wy2 = 0. END WHERE WHERE((x1 .EQ. x2) .AND. (y1 .EQ. y2)) wx1 = 1. wx2 = 0. wy1 = 1. wy2 = 0. END WHERE c Main calculation, loop over each element of the input. DO 10 jj=1,NPts f1 = wx1(jj)*Image(y1(jj),x1(jj)) + wx2(jj)*Image(y1(jj),x2(jj)) f2 = wx1(jj)*Image(y2(jj),x1(jj)) + wx2(jj)*Image(y2(jj),x2(jj)) IntIm(jj) = wy1(jj)*f1 + wy2(jj)*f2 10 CONTINUE RETURN END