n维数字图像欧氏距离变换算法

算法简介

该算法主要用于在三维图像中计算有效体素之间的最短欧几里得距离。

算法出处：NEW ALGORITHMS FOR EUCLIDEAN DISTANCE TRANSFORMATION OF AN n-DIMENSIONAL DIGITIZED PICTURE WITH APPLICATIONS
算法在体渲染加速中的应用：Accelerated Volume Rendering with Chebyshev Distance Maps

算法过程

公式表示

符号说明
图像维度：$L\times M\times N$，其中$1\le i\le L,1\le j\le M,1\le k\le N$
输入图像： F = { f i j k } F=\{f_{ijk}\} F={fijk​}
中间图像： G = { g i j k } G=\{g_{ijk}\} G={gijk​}， H = { h i j k } H=\{h_{ijk}\} H={hijk​}
输出图像： S = { s i j k } S=\{s_{ijk}\} S={sijk​}

步骤1

输入：F
输出：G
公式：当$f_{xjk}=0,1\le x\le L$，求$g_{ijk}=min\{(i-x{)}^2\}$

$f_{xjk}=0$ 表示在图像中体素值为0是有效体素

计算示例图：

步骤2

输入：G
输出：H
公式：当$1\le y\le M$，求 h i j k = m i n { g i j k + ( j − y ) 2 } h_{ijk}=min\{ g_{ijk} + (j-y)^2 \} hijk​=min{gijk​+(j−y)2}

计算示例图：

步骤3

输入：H
输出：S
公式：当$1\le z\le N$，求 s i j k = m i n { h i j k + ( k − z ) 2 } s_{ijk}=min\{ h_{ijk} + (k-z)^2 \} sijk​=min{hijk​+(k−z)2}

计算过程和步骤2类似，只是换了维度，因此不再赘述

算法示例

为了描述方便，这里使用二维图像作为操作目标，且每个像素是二值的。由于是二维的，所以步骤会少一步，即无步骤3，最终的输出结果是H

原始图像F，值为0表示有效像素

中间图像G，因为最后一行无有效像素，所以保持默认值

最终图像H，像素内存储最近有效像素的距离的平方

我们将图像H中的值开方，即可得到每个像素的最近有效像素的欧式距离

代码实现

在计算机上实现

步骤1

前向扫描
输入：F
输出：G’
设$g_{0jk}^{\prime }=L^2$
要求$i$从$1$至$L$，如果$f_{ijk}\ne 0$，则$g_{ijk}^{\prime }=(\sqrt{g_{(i-1)jk}^{\prime }}+1{)}^2$；否则$g_{ijk}^{\prime }=0$

反向扫描
输入：G’
输出：G
设$g_{(L+1)jk}=L^2$
要求$i$从$L$至$1$，如果$g_{ijk}^{\prime }\ne 0$，则$g_{ijk}=min\{(\sqrt{g_{(i+1)jk}}+1{)}^2,g_{ijk}^{\prime }\}$；否则$g_{ijk}=0$

步骤2

输入：G
输出：H

要求$n$从$-r_1$至$r_2$，其中$r_1\le r,r_2\le r,r=\sqrt{g_{ijk}}$，求 h i j k = m i n { g i ( j + n ) k + n 2 } h_{ijk}=min\{ g_{i(j+n)k} + n^2 \} hijk​=min{gi(j+n)k​+n2}
注意$1\le j+n\le M$，所以要保证$j-r_1\ge 1,j+r_2\le M$，从而$r_1\le j-1,r_2\le M-j$
所以 r 1 = m i n { r , j − 1 } , r 2 = m i n { r , M − j } r_1 = min \{ r, j-1 \},r_2 = min\{ r,M-j \} r1​=min{r,j−1},r2​=min{r,M−j}

步骤3

输入：H
输出：S

要求$n$从$-r$至$r$，其中$r=\sqrt{h_{ijk}}$，求 s i j k = m i n { h i j ( k + n ) + n 2 } s_{ijk}=min\{ h_{ij(k+n)}+ n^2 \} sijk​=min{hij(k+n)​+n2}
注意$1\le k+n\le N$，所以要保证$k-r_1\ge 1,k+r_2\le N$，从而$r_1\le k-1,r_2\le M-k$
所以 r 1 = m i n { r , k − 1 } , r 2 = m i n { r , M − k } r_1 = min \{ r, k-1 \},r_2 = min\{ r,M-k \} r1​=min{r,k−1},r2​=min{r,M−k}

编写C++代码

输入数据是三维向量

void edt(vector<vector<vector<float>>>& f)
{
	int N = f.size();
	int M = f[0].size();
	int L = f[0][0].size();

	///步骤1
	///前向扫描
	for (int k = 0; k < N; k++)
	{
		for (int j = 0; j < M; j++)
		{
			float df = L;
			for (int i = 0; i < L; i++)
			{
				if (f[k][j][i] != 0)
				{
					df = df + 1;
				}
				else
				{
					df = 0;
				}
				f[k][j][i] = df * df;
			}
		}
		///反向扫描
		for (int j = 0; j < M; j++)
		{
			float db = L;
			for (int i = L - 1; i >= 0; i--)
			{
				if (f[k][j][i] != 0)
				{
					db = db + 1;
				}
				else
				{
					db = 0;
				}
				f[k][j][i] = std::fmin(f[k][j][i], db * db);
			}
		}
	}
	

	///步骤2
	for (int k = 0; k < N; k++)
	{		
		for (int i = 0; i < L; i++)
		{
			vector<float> buff(M);
			for (int j = 0; j < M; j++)
			{
				buff[j] = f[k][j][i];
			}

			for (int j = 0; j < M; j++)
			{
				float d = buff[j];
				if (d != 0)
				{
					float rMax = int(std::sqrt(d)) + 1;
					float rStart = std::min(rMax, float(j));
					float rEnd = std::min(rMax, float(M - 1 - j));

					for (int n = -rStart; n < rEnd; n++)
					{
						float w = buff[j + n] + n * n;
						if (w < d)
						{
							d = w;
						}
					}
				}
				f[k][j][i] = d;
			}
		}
	}
	
	///步骤3，与步骤2类似
	for (int j = 0; j < M; j++)
	{
		for (int i = 0; i < L; i++)
		{
			vector<float> buff(N);
			for (int k = 0; k < N; k++)
			{
				buff[k] = f[k][j][i];
			}

			for (int k = 0; k < N; k++)
			{
				float d = buff[k];
				if (d != 0)
				{
					float rMax = int(std::sqrt(d)) + 1;
					float rStart = std::min(rMax, float(k));
					float rEnd = std::min(rMax, float(N - 1 - k));

					for (int n = -rStart; n < rEnd; n++)
					{
						float w = buff[k + n] + n * n;
						if (w < d)
						{
							d = w;
						}
					}
				}
				f[k][j][i] = d;
			}
		}
	}
}

测试代码

#include<iostream>
#include<cmath>
#include<vector>
using namespace std;

int main(int argc, char *argv[])
{
	//图像维度4x4x1
	vector<vector<vector<float>>> F(1, vector<vector<float>>(4, vector<float>(4)));
	
	//输入图像数据
	F[0][0][0] = 1; F[0][0][1] = 1; F[0][0][2] = 1; F[0][0][3] = 0;
	F[0][1][0] = 1; F[0][1][1] = 0; F[0][1][2] = 1; F[0][1][3] = 0;
	F[0][2][0] = 0; F[0][2][1] = 1; F[0][2][2] = 1; F[0][2][3] = 1;
	F[0][3][0] = 1; F[0][3][1] = 1; F[0][3][2] = 1; F[0][3][3] = 1;

	edt(F);

	//输出计算数据
	for (int k = 0; k < F.size(); k++)
	{
		for (int j = 0; j < F[k].size(); j++)
		{
			for (int i = 0; i < F[k][j].size(); i++)
			{
				cout << F[k][j][i] << " ";
			}
			cout << endl;
		}		
	}

	return 0;
}

结果