我有一个3D矩阵（维度nx，nz，ny），它对应一个物理域。该矩阵包含从-1（阶段1）到+1（阶段2）的连续字段;两个阶段之间的接口是该字段的级别0。
现在，我想有效地计算域中每个点的接口的带符号距离函数。
我尝试了两种可能性（sgn是我的字段的符号，值为+ 1,0，-1，xyz包含网格作为每个点的x，y，z的三元组，dist是我想要计算的符号距离函数） 。
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
我还尝试了第二个版本，它与上面的版本非常相似，但它仅在整个矩阵的子集中计算欧氏距离。有了这个第二个版本，有一些加速，但它仍然太慢。我想知道是否有更有效的方法来计算距离函数。
第二版：
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：点i，k，j与框中任何其他点之间的欧几里德距离2 * deltax + 1,2 * deltaz + 1,2 * deltay + 1以i，k，j为中心。这降低了计算成本，因为距离仅在整个网格的子集中计算（这里我假设子集足够大以包含界面点）。
在Vladimir建议之后（x，y，z是确定网格位置的轴，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
编辑：有关手头问题的更多信息。网格是映射到矩阵的正交网格。该网格的点数在每个方向上大约为1000（总共大约10亿个点）。
我的目标是以有效的方式从整个网格中的符号函数（+ 1,0，-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 = 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
考虑将外部三重循环内的整个事物分成子例程或函数。它不会更快，但它会更具可读性。特别是对我们来说，这里仅仅处理该函数就足够了，外部循环只是一个令人困惑的额外层。