- #1
octopode
- 14
- 0
Homework Statement
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.
Homework Equations
The equations for bilinear interpolation that I'm using are :
Code:
f1 = wx1*Image(y1,x1) + wx2*Image(y1,x2)
f2 = wx1*Image(y2,x1) + wx2*Image(y2,x2)
IntIm= wy1*f1 + wy2*f2x1 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
Code:
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