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. My 1st attempt work but is slow. See code below. 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 using the C++ as reference for the speedup that Fortran can or cannot offer.
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 :: 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
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
[EDIT] As i only have regular grid, x2-x1 = 1 or 0, i've made another version that does not require division. This got rid of the where constructs. In addition, it seems that replace the DO loop with a "FORALL" loop improved by about 15% on 2 cpus. My understanding is that FORALL uses multithread when present, and each computation at each jj is independent from the other at jj+k or jj-k which makes the FORALL statement sensible. Doesn't it ?
SUBROUTINE L2DINTERPOL(IntIm,Image,x,y,NPts,M,N)
implicit none
mwSize, PARAMETER :: dp = kind(0.d0) ! Double precision
mwSize :: NPts, M,N ! Input
REAL(dp),DIMENSION(Npts) :: x,y ! Input
REAL(dp),DIMENSION(M,N) :: Image ! Input
mwSize :: jj
mwSize, DIMENSION(NPts) :: x1,y1,x2,y2
REAL(dp),DIMENSION(Npts) :: wx1,wx2,wy1,wy2
REAL(dp),DIMENSION(Npts) :: IntIm ! Output
x1 = FLOOR(x)
x2 = CEILING(x)
y1 = FLOOR(y)
y2 = CEILING(y)
wx1 = x2-x
wy1 = y2-y
wx2 = 1 - wx1
wy2 = 1 - wy1
FORALL( jj=1 : NPts)
IntIm(jj) = wy1(jj)*(wx1(jj)*Image(y1(jj),x1(jj))+wx2(jj)*
$ Image(y1(jj),x2(jj))) + wy2(jj)*(wx1(jj)*
$ Image(y2(jj),x1(jj))+wx2(jj)*Image(y2(jj),x2(jj)))
END FORALL
RETURN
END
NPts
) do you have? What are the dimensions (M
,N
) of the image? So that we can make sure that we benchmark your actual problem. I wrote a trilinear interpolation (in 3d) here: github.com/certik/hfsolver/blob/master/src/interp3d.f90, so you can compare the speed against it. How are the points inx
,y
distributed? Is it a regular mesh? If so, then you can exploit this feature to possibly speed things up. $\endgroup$