Efficient way to calculate distance function

Question

I have a 3D matrix (dimension nx,nz,ny) which corresponds to a physical domain. This matrix contains a continuous field from -1 (phase 1) to +1 (phase 2); the interface between the two phases is the level 0 of this field.

Now, I want to calculate efficiently the signed distance function from the interface for every point in the domain.

I tried two possibilities (sgn is the sign of my field, with values +1,0,-1, xyz contains the grid as triplets of x,y,z at each point and dist is the signed distance function I want to calculate).

double precision, dimension(nx,nz,ny) :: dist,sgn,eudist
integer :: i,j,k
double precision :: seed,posit,tmp(nx)

do j=1,ny  
do k=1,nz  
 do i=1,nx  
  seed=sgn(i,k,j)  
    ! look for interface  
    eudist=(xyz(:,:,:,1)-x(i))**2+(xyz(:,:,:,2)-z(k))**2+(xyz(:,:,:,3)-y(j))**2  
    ! find min within mask  
    posit=minval(eudist,seed*sgn.le.0)  
    ! tmp fits in cache, small speed-up  
    tmp(i)=-seed*dsqrt(posit)  
  enddo  
  dist(:,k,j)=tmp  
enddo  
enddo  

I also tried a second version, which is quite similar to the above one but it calculates the Euclidean distance only in a subset of the whole matrix. With this second version there is some speed up, but it is still too slow. I would like to know whether there is a more efficient way to calculate the distance function.

Second version:

double precision, dimension(nx,nz,ny) :: dist,sgn
double precision, allocatable, dimension(:,:,:) :: eudist
integer :: i,j,k  , ii,jj,kk
integer :: il,iu,jl,ju,kl,ku
double precision :: seed, deltax,deltay,deltaz,tmp(nx)

deltax=max(int(nx/4),1)
deltay=max(int(ny/4),1)
deltaz=max(int(nz/2),1)

allocate(eudist(2*deltax+1,2*deltaz+1,2*deltay+1))

do j=1,ny
 do k=1,nz
  do i=1,nx
   ! look for closest point in box 2*deltax+1,2*deltaz+1,2*deltay+1
   il=max(1,i-deltax)
   iu=min(nx,i+deltax)
   jl=max(1,j-deltay)
   ju=min(ny,j+deltay)
   kl=max(1,k-deltaz)
   ku=min(nz,k+deltaz)

   eudist(:,1:ku-kl+1,:)=(xyz(il:iu,kl:ku,jl:ju,1)-x(i))**2 &
 &                      +(xyz(il:iu,kl:ku,jl:ju,2)-z(k))**2 &
 &                      +(xyz(il:iu,kl:ku,jl:ju,3)-y(j))**2

   seed=sgn(i,k,j)
   tmp(i)=minval(eudist(:,1:ku-kl+1,:),seed*sgn(il:iu,kl:ku,jl:ju).le.0)
   tmp(i)=-seed*dsqrt(tmp(i))

  enddo
  dist(:,k,j)=tmp
 enddo
enddo

eudist: Euclidean distance between the point i,k,j and any other point in a box 2*deltax+1,2*deltaz+1,2*deltay+1 centered in i,k,j. This reduces computational cost, as the distance is calculated only in a subset of the whole grid (here I am assuming that the subset is large enough to contain an interfacial point).

After Vladimir suggestion (x,y,z are the axes determining grid position, xyz(i,k,j)=(x(i),z(k),y(j)) ):

double precision, dimension(nx,nz,ny) :: dist,sgn
double precision :: x(nx), y(ny), z(nz)
double precision, allocatable, dimension(:,:,:) :: eudist
double precision, allocatable, dimension(:) :: xd,yd,zd
integer :: i,j,k  , ii,jj,kk
integer :: il,iu,jl,ju,kl,ku
double precision :: seed, deltax,deltay,deltaz,tmp(nx)

deltax=max(int(nx/4),1)
deltay=max(int(ny/4),1)
deltaz=max(int(nz/2),1)

allocate(eudist(2*deltax+1,2*deltaz+1,2*deltay+1))
allocate(xd(2*deltax+1))
allocate(yd(2*deltay+1))
allocate(zd(2*deltaz+1))

do j=1,ny
 do k=1,nz
  do i=1,nx
   ! look for closest point in box 2*deltax+1,2*deltaz+1,2*deltay+1
   il=max(1,i-deltax)
   iu=min(nx,i+deltax)
   jl=max(1,j-deltay)
   ju=min(ny,j+deltay)
   kl=max(1,k-deltaz)
   ku=min(nz,k+deltaz)

   do ii=1,iu-il+1
     xd(ii)=(xyz(il+ii-1)-x(i))**2
   end do

   do jj=1,ju-jl+1
     yd(jj)=(y(jj+jl-1)-y(j))**2
   end do

   do kk=1,ku-kl+1
     zd(kk)=(z(kk+kl-1)-z(k))**2
   end do

   do jj=1,ju-jl+1
     do kk=1,ku-kl+1
       do ii=1,iu-il+1
         eudist(ii,kk,jj)=xd(ii)+yd(jj)+zd(kk)
       enddo
     enddo
   enddo

   seed=sgn(i,k,j)
   tmp(i)=minval(eudist(:,1:ku-kl+1,:),seed*sgn(il:iu,kl:ku,jl:ju).le.0)
   tmp(i)=-seed*dsqrt(tmp(i))

  enddo
  dist(:,k,j)=tmp
 enddo
enddo

EDIT: more information on the problem at hand. The grid is an orthogonal grid mapped to a matrix. The number of points of this grid is of the order of 1000 in each direction (in total about 1 billion points).

My goal is switching from a sign function (+1,0,-1) to a signed distance function in the entire grid in an efficient way.

Answer 1

I would still do what I suggested, no matter if you do that on a subset or across the whole plane. Take advantage of the orthogonal grid, it is a great thing to have

do j=1,ny
 do k=1,nz
  do i=1,nx
   ! look for closest point in box 2*deltax+1,2*deltaz+1,2*deltay+1
   il=max(1,i-deltax)
   iu=min(nx,i+deltax)
   jl=max(1,j-deltay)
   ju=min(ny,j+deltay)
   kl=max(1,k-deltaz)
   ku=min(nz,k+deltaz)


   do ii = il,iu
     xd(i) = (xyz(ii,kl:ku,jl:ju,1)-x(i))**2
   end do

   do jj = jl,ju
     yd(i) = (xyz(il:iu,kl:ku,jj,2)-y(j))**2
   end do

   do kk = kl,ku
     zd(k) = (xyz(il:iu,kk,jl:ju,3)-z(k))**2
   end do

   do jj = jl,ju
     do kk = kl,ku
       do ii = il,iu
         eudist(il:iu,kl:ku,jl:ju) = xd(ii) + yd(jj) + zd(kk)
       end do
     end do
   end do
   ....

  enddo
  dist(:,k,j)=tmp
 enddo
enddo

Consider separating the whole thing that is inside the outer triple loop into a subroutine or a function. It would not be faster, but it would be much more readable. Especially for us here, It would be enough for us here to only deal with that function, the outer loop is just a confusing extra layer.