Euclidean distance map（EDM)这个概念可能听过的人也很少，其主要是用在二值图像中，作为一个很有效的中间处理手段存在。一般的处理都是将灰度图处理成二值图或者一个二值图处理成另外一个二值图，而EDM算法确是由一幅二值图生成一幅灰度图。其核心定义如下：
The definition is simple enough: each point in the foreground is assigned a brightness value equal to its straight line (hence “Euclidean”) distance from the nearest point in the background.
或者用一句更简洁的话说就是：The value of each pixel is the distance to the nearest background pixel (for background pixels, the EDM is 0)。
EDM由很多用处，比如由于像素都是分布在矩形格上，造成一些形态学上的操作存在一些方向偏移，以及在距离计算上的一些偏差都可以使用EDM技术克服，使用EDM来实现的腐蚀、膨胀或者开和闭操作相比普通的逐像素（模板操作）的结果来的更加各向同性化，如下图所示（很明显，基于EDM处理的图结果更加光滑，没有任何的方向性，而传统的逐像素的算法则有点棱角分明）：
图一： 使用普通的逐像素的处理（从做至右依次是：二值原图、闭操作、开操作，找到的边界）
图二： 使用基于EDM的处理（从做至右依次是：二值原图、闭操作、开操作，找到的边界）
直接从图像的前景点搜索和其最接近的背景点的距离是一个极其低效和耗时的操作。一些研究者通过只取几个方向距离也实现了一些不同的EDM效果，比如只有45或者90度的角度的距离，但是这些距离的效果和传统的一些处理方法同样存在这没有足够的各向同性的性质，因此，也不具有很大的意义，如下图所示：
为了快速的实现EDM值的计算，不少作者提出了一些快速算法，其中比较经典的是Danielsson在1980提出的算法，其通过两次遍历图像给出了和原始计算一样的效果，具体步骤如下所示：
1. Assign the brightness value of 0 to each pixel in the background and a large positive value (greater than the maximum feature width) to each pixel in a feature.
2. Proceeding from left to right and top to bottom, assign each pixel within a feature a brightness value one greater than the smallest value of any of its neighbors.
3. Repeat step 2, proceeding from right to left and bottom to top.
在相关的资料中寻找参考的代码，可以找到如下一些资料：
The EDM algorithm is similar to the 8SSEDT in F. Leymarie, M. D. Levine, in: CVGIP Image Understanding, vol. 55 (1992), pp 84-94
https://linkinghub.elsevier.com/retrieve/pii/104996609290008Q
最直接的就是ImageJ了，在其开源的目录：ImageJ\source\ij\plugin\filter下有Edm.java文件分享EDM算法的实现，不过那个代码写得比较恶心，我后面也不知道在什么地方找到了一个比较清洁的代码（也许是我自己改清洁的），核心思想是和这里的一致的：
public static int ONE = ; // 这个值和一个像素的距离对应
public static int SQRT2 = ; // 这个值和Sqr(2)个像素的距离对应
public static int SQRT5 = ; // 这个值和Sqr(5)个像素的距离对应 public static void CalaEuclideanDistanceMap(FastBitmap Bmp)
{
int X, Y, Index, Value,Max;
int Width, Height, Xmax, Ymax ;
byte* Pointer; Width = Bmp.Width; Height = Bmp.Height; Xmax = Width - ; Ymax = Height - ; // 1. Assign the brightness value of 0 to each pixel in the background and a large positive
// value (greater than the maximum feature width) to each pixel in a feature. int[] ImageData = new int[Width * Height];
for (Y = , Index = ; Y < Height; Y++)
{
Pointer = Bmp.Pointer + Y * Bmp.Stride;
for (X = ; X < Width; X++, Index++) ImageData[Index] = Pointer[X]<<;
} // 2. Proceeding from left to right and top to bottom, assign each pixel within a feature a
// brightness value one greater than the smallest value of any of its neighbors. for (Y = ,Index=; Y < Height; Y++)
{
Pointer = Bmp.Pointer + Y * Bmp.Stride;
for (X = ; X < Width; X++,Index++)
{
if (Pointer[X] > )
{
if ((X <= ) || (X >= Xmax) || (Y <= ) || (Y >= Ymax))
SetEdgeValue(Index, Width, ImageData, X, Y, Xmax, Ymax);
else
SetValue(Index, Width, ImageData);
}
}
} // 3. Repeat step 2, proceeding from right to left and bottom to top.
for (Y = Height - ,Index=Width*Height-; Y >= ; Y--)
{
Pointer = Bmp.Pointer + Y * Bmp.Stride;
for (X = Width - ; X >= ; X--,Index--)
{
if (Pointer[X] > )
if ((X <= ) || (X >= Xmax) || (Y <= ) || (Y >= Ymax))
SetEdgeValue(Index, Width, ImageData, X, Y, Xmax, Ymax);
else
SetValue(Index, Width, ImageData);
}
} // Find the max value of the data
for (Y = ,Max=, Index = ; Y < Height; Y++)
for (X = ; X < Width; X++, Index++)
if (Max < ImageData[Index]) Max = ImageData[Index]; for (Y = , Index = ; Y < Height; Y++)
{
Pointer = Bmp.Pointer + Y * Bmp.Stride;
for (X = ; X < Width; X++)
{
Value = (ImageData[Index]) * / Max;
Pointer[X] = (byte)Value;
Index++;
}
}
} private static void SetValue(int Offset, int Width, int[] ImageData)
{
int Value;
int Index1 = Offset - Width - Width - ;
int Index2 = Index1 + Width;
int Index3 = Index2 + Width;
int Index4 = Index3 + Width;
int Index5 = Index4 + Width;
int Min = int.MaxValue;
// *
// * x *
// *
Value = ImageData[Index2 + ] + ONE;
if (Value < Min) Min = Value;
Value = ImageData[Index3 + ] + ONE;
if (Value < Min) Min = Value;
Value = ImageData[Index3 + ] + ONE;
if (Value < Min) Min = Value;
Value = ImageData[Index4 + ] + ONE;
if (Value < Min) Min = Value; // * *
// x
// * *
Value = ImageData[Index2 + ] + SQRT2;
if (Value < Min) Min = Value;
Value = ImageData[Index2 + ] + SQRT2;
if (Value < Min) Min = Value;
Value = ImageData[Index4 + ] + SQRT2;
if (Value < Min) Min = Value;
Value = ImageData[Index4 + ] + SQRT2;
if (Value < Min) Min = Value; // * *
// * *
// x
// * *
// * *
Value = ImageData[Index1 + ] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index1 + ] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index2 + ] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index4 + ] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index5 + ] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index5 + ] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index4] + SQRT5;
if (Value < Min) Min = Value;
Value = ImageData[Index2] + SQRT5;
if (Value < Min) Min = Value; ImageData[Offset] = Min; } private static void SetEdgeValue(int Offset, int Width, int[] ImageData, int X, int Y, int xmax, int ymax)
{ }
我删除了边缘部分的代码（不希望懒人直接使用，直接拷贝过去就能用），代码也很简单，第一步是将原始图像数据放大，放大的尺度文章中讲的一定要大于最大的特征的尺度，其实就是图像中前景和背景最远的那个距离值，实际上这个值不可能超过图像的对角线长度，一般我们随便放大一个倍数就可以了，上述代码放大了65536倍，完全足够了，实际还可以缩小的，对于背景色，原始值为0，放大后还是0。
那么第二步，是从做到右，从上到下进行处理，处理时的原则先计算半径为1时周边的4个点，为了避免浮点计算，我们把距离的定义放大了100倍，比如距离为1，就变为100，距离为根号2，就变为141。这四个点如果他们的值有一个黑色，则肯定就是这四个点中最小的值了，也是靠近改点的最小的距离了，如果四个点都是白色，我们在把半径扩展到根号2，查看这个时候四个点的颜色情况，注意或者时候求最小值时需要加上141了，接着在看看根号5时的情况了。
为什么只需要求这三个距离的最小值，我还一时没有想清楚，但是这里有一点就是，这个处理过程对每个点是有先后依赖顺序的，我们看到在SetValue函数中当前点ImageData的更新会影响到下一个点的状况的，这也是有道理的，因为当前点的值处理完后就表示这个点距离最近的背景的距离，那么下一个点计算时其和当前的距离在加上当前点和背景的距离就反应了下一个点和背景的距离。而这种前后依赖的关系对于我们的后续的SSE造成了一定的影响。
我们看下这段代码的显示结果：
原图 二值化 EDM图
EDM图中越是白色的地方说明此处距离背景越远，而途中的红色圆圈则表示一个小小的孤立黑点也会对EDM图产生不利的影响，比如，我们如果把二值化后的图进行简单的去除白色和黑色孤点处理后，在进行EDM处理，则效果会平滑很多，如下所示：
二值化 去除孤点 EDM图
此时的EDM图中清晰的可以看到5根明显的分界线，这对于后续的一些分割等等识别工作很有裨益。
为测试速度，我们选择了一幅3000*2000的二值图进行测试， 因为这个算法的时间是和图像内容有一定关系的（黑色的部分不需要处理），所以我们的二值图中有一半的颜色是白色，大概处理时间在160ms左右（上述C#的代码)，而且未做太多的代码优化，应该说时间还是很快的。
时间没有极限，我还是希望能尝试进一步加速，于是我还是构思了实现了改算法的SSE优化。
先从简单的搞起，第一，我前面说过放大的系数没有必要到65536倍，只要这里的值大于了对角线大小就OK了，比如我们放大64倍，那白色的值就变为255 * 64 = 16320，足够平常使用了，这个时候有个好处就是我们可以不用int类型数据来记录中间的距离了，而可以直接使用unsigned short类型。此时最后一部分的量化我们可飞快的用SSE实现：
void IM_GetMinMaxValue(unsigned short *Src, int Length, unsigned short &Min, unsigned short &Max)
{
int BlockSize = , Block = Length / BlockSize;
__m128i MinValue = _mm_set1_epi16();
__m128i MaxValue = _mm_set1_epi16();
for (int Y = ; Y < Block * BlockSize; Y += BlockSize)
{
__m128i SrcV = _mm_loadu_si128((__m128i *)(Src + Y));
MinValue = _mm_min_epu16(MinValue, SrcV);
MaxValue = _mm_max_epu16(MaxValue, SrcV);
}
Min = _mm_extract_epi16(_mm_minpos_epu16(MinValue), );
Max = - _mm_extract_epi16(_mm_minpos_epu16(_mm_subs_epu16(_mm_set1_epi16(), MaxValue)), ); // 经过测试结果正确
for (int Y = Block * BlockSize; Y < Length; Y++)
{
if (Min > Src[Y]) Min = Src[Y];
if (Max < Src[Y]) Max = Src[Y];
}
}
int IM_Normalize(unsigned short *Src, unsigned char*Dest, int Width, int Height, bool BmpFormat = false)
{
if ((Src == NULL) || (Dest == NULL)) return IM_STATUS_NULLREFRENCE;
if ((Width <= ) || (Height <= )) return IM_STATUS_INVALIDPARAMETER; int Stride = BmpFormat == false ? Width : WIDTHBYTES(Width);
unsigned short Min = , Max = ;
IM_GetMinMaxValue(Src, Width * Height, Min, Max);
if (Max == Min)
{
memset(Dest, , Height * Stride);
}
else
{
float Inv = 255.0f / (Max - Min);
int BlockSize = , Block = Width / BlockSize;
__m128i Zero = _mm_setzero_si128();
for (int Y = ; Y < Height; Y++)
{
unsigned short *LinePS = Src + Y * Width;
unsigned char *LinePD = Dest + Y * Stride;
for (int X = ; X < Block * BlockSize; X += BlockSize) // 正规化
{
__m128i SrcV = _mm_loadu_si128((__m128i *)(LinePS + X));
__m128i Diff = _mm_subs_epu16(SrcV, _mm_set1_epi16(Min));
__m128i Value1 = _mm_cvtps_epi32(_mm_mul_ps(_mm_cvtepi32_ps(_mm_cvtepu16_epi32(Diff)), _mm_set1_ps(Inv)));
__m128i Value2 = _mm_cvtps_epi32(_mm_mul_ps(_mm_cvtepi32_ps(_mm_cvtepu16_epi32(_mm_srli_si128(Diff, ))), _mm_set1_ps(Inv)));
_mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_packus_epi32(Value1, Value2), Zero));
}
for (int X = Block * BlockSize; X < Width; X++)
{
LinePD[X] = (int)((LinePS[X] - Min) * Inv + 0.5f);
}
}
}
return IM_STATUS_OK;
}
因为在原始的代码中，最小值必然为0，所以原始的C#代码没有计算最小值，而本代码为了通用性，也计算了最小值，其过程也相当简单，就是_mm_max_epu16和_mm_min_epu16的调用，而由于数据是unsigned short的，也可以快捷的使用_mm_minpos_epu16（其他数据类型都没有这个函数哦）获得8个ushort类型的最小值。后续的量化到unsigned char过程的也基本就是普通函数的调用，不存在其他技巧，主要是要注意_mm_cvtepu16_epi32这种数据类型的转换函数，要用的恰当。
注意到前面讲的ImageData的更新是有前后依赖的，这就非常不利于我一次性计算多个像素的值，从而类似于SSE图像算法优化系列九：灵活运用SIMD指令16倍提升Sobel边缘检测的速度（4000*3000的24位图像时间由480ms降低到30ms）这篇文章的优化方式就无法直接复现。但是我们仔细观察，发现在SetValue函数中，除了在半径为1的领域求最小值的计算中涉及到了本行值，其他的在一行值的更新过程中是不相互干涉的，那么我们可以把这些单独提取出来计算作为中间值存储，然后在用普通的C代码和领域为1的左右两个位置的量进行比较。这样一样能起到不小的加速作用。
当然要注意一点，由于这种依赖关系，在跳到下行计算时，更新过的这一行数据必须得以保留，否则又会得到错误的结果。
我们采用类似Sobel优化博文中的优化方式，采用了5行临时数据缓冲区来解决边缘处的处理问题，这样就无需为边缘处在单独写代码，同时，我们没有必要先对扩大后的数据进行计算，而时在5行数据更新的同时更新此部分数据，这样在整个过程中可以少一部分拷贝工作。
// 1. Assign the brightness value of 0 to each pixel in the background and a large positive
// value (greater than the maximum feature width) to each pixel in a feature.
__forceinline void IM_GetExpandAndZoomInData(unsigned char *Src, unsigned short *Dest, int Width, int Height)
{
const int BlockSize = , Block = Width / BlockSize;
for (int X = ; X < Block * BlockSize; X += BlockSize)
{
_mm_storeu_si128((__m128i *)(Dest + X + ), _mm_slli_epi16(_mm_cvtepu8_epi16(_mm_loadl_epi64((__m128i *)(Src + X))), ));
}
for (int X = Block * BlockSize; X < Width; X++) Dest[X + ] = Src[X] << ;
Dest[] = Dest[]; Dest[] = Dest[Width > ? : ]; // 镜像数据
Dest[Width + ] = Dest[Width + ]; Dest[Width + ] = Dest[Width > ? Width : ];
}
从左到右，从右到左方向的更新。
int IM_ShowEuclideanDistanceMap(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
int Channel = Stride / Width;
if ((Src == NULL) || (Dest == NULL)) return IM_STATUS_NULLREFRENCE;
if ((Width <= ) || (Height <= )) return IM_STATUS_INVALIDPARAMETER;
if (Channel != ) return IM_STATUS_INVALIDPARAMETER;
if (IM_IsBinaryImage(Src, Width, Height, Stride) == false) return IM_STATUS_INVALIDPARAMETER;
int Status = IM_STATUS_OK;
int ExpandW = + Width + ; const int ONE = ; // 这个值和一个像素的距离对应
const int SQRT2 = ; // 这个值和Sqr(2)个像素的距离对应
const int SQRT5 = ; // 这个值和Sqr(5)个像素的距离对应
const int BlockSize = , Block = Width / BlockSize;
unsigned short *Distance = (unsigned short *)malloc(Width * Height * sizeof(unsigned short));
unsigned short *Buffer = (unsigned short *)malloc((ExpandW * + Width) * sizeof(unsigned short)); // 5行缓冲用于记录相邻5行取样数据，1行用于记录中间结果
if ((Distance == NULL) || (Buffer == NULL))
{
Status = IM_STATUS_OUTOFMEMORY;
goto FreeMemory;
} unsigned short *First = Buffer, *Second = First + ExpandW, *Third = Second + ExpandW;
unsigned short *Fourth = Third + ExpandW, *Five = Fourth + ExpandW, *Temp = Five + ExpandW; // 2. Proceeding from left to right and top to bottom, assign each pixel within a feature a
// brightness value one greater than the smallest value of any of its neighbors. IM_GetExpandAndZoomInData(Src, First, Width, Height); // 把他们写在这个过程中，可以减少一次赋值
IM_GetExpandAndZoomInData(Src, Second, Width, Height);
IM_GetExpandAndZoomInData(Src, Third, Width, Height); // 前面三行数据相同
IM_GetExpandAndZoomInData(Src + IM_ClampI(, , Height - ) * Stride, Fourth, Width, Height);
IM_GetExpandAndZoomInData(Src + IM_ClampI(, , Height - ) * Stride, Five, Width, Height); for (int Y = ; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
if (Y != )
{
unsigned short *Temp = First; First = Second; Second = Third; Third = Fourth; Fourth = Five; Five = Temp; // 交换指针
}
if (Y == Height - ) // 倒数第二行
{
memcpy(Five, Fourth, ExpandW * sizeof(unsigned short));
}
else if (Y == Height - ) // 最后一行，不能用第三行的数据，因为第三行是修改后的
{
IM_GetExpandAndZoomInData(Src + (Height - ) * Stride, Fourth, Width, Height);
IM_GetExpandAndZoomInData(Src + IM_ClampI(Height - , , Height - ) * Stride, Five, Width, Height);
}
else
{
IM_GetExpandAndZoomInData(Src + (Y + ) * Stride, Five, Width, Height);
}
for (int X = ; X < Block * BlockSize; X += BlockSize)
{
__m128i SrcV = _mm_loadl_epi64((__m128i *)(LinePS + X));
if (_mm_movemask_epi8(SrcV) != ) // 8个字节全是白色，则不需要继续处理
{
__m128i FirstP1 = _mm_loadu_si128((__m128i *)(First + X + ));
__m128i FirstP3 = _mm_loadu_si128((__m128i *)(First + X + )); __m128i SecondP0 = _mm_loadu_si128((__m128i *)(Second + X + ));
__m128i SecondP1 = _mm_loadu_si128((__m128i *)(Second + X + ));
__m128i SecondP2 = _mm_loadu_si128((__m128i *)(Second + X + ));
__m128i SecondP3 = _mm_loadu_si128((__m128i *)(Second + X + ));
__m128i SecondP4 = _mm_loadu_si128((__m128i *)(Second + X + )); __m128i FourthP0 = _mm_loadu_si128((__m128i *)(Fourth + X + ));
__m128i FourthP1 = _mm_loadu_si128((__m128i *)(Fourth + X + ));
__m128i FourthP2 = _mm_loadu_si128((__m128i *)(Fourth + X + ));
__m128i FourthP3 = _mm_loadu_si128((__m128i *)(Fourth + X + ));
__m128i FourthP4 = _mm_loadu_si128((__m128i *)(Fourth + X + )); __m128i FiveP1 = _mm_loadu_si128((__m128i *)(Five + X + ));
__m128i FiveP3 = _mm_loadu_si128((__m128i *)(Five + X + ));
// *
// * x *
// *
__m128i Min0 = _mm_min_epu16(SecondP2, FourthP2); // 因为第三行是前后相关的，所以不能在这里参与计算 // * *
// x
// * *
__m128i Min1 = _mm_min_epu16(_mm_min_epu16(SecondP1, SecondP3), _mm_min_epu16(FourthP1, FourthP3)); // * *
// * *
// x
// * *
// * *
__m128i Min2 = _mm_min_epu16( _mm_min_epu16(_mm_min_epu16(FirstP1, FirstP3), _mm_min_epu16(SecondP0, SecondP4)),
_mm_min_epu16(_mm_min_epu16(FourthP0, FourthP4), _mm_min_epu16(FiveP1, FiveP3))); __m128i Min = _mm_min_epu16(_mm_min_epu16(_mm_adds_epu16(Min0, _mm_set1_epi16(ONE)), _mm_adds_epu16(Min1, _mm_set1_epi16(SQRT2))), _mm_adds_epu16(Min2, _mm_set1_epi16(SQRT5))); _mm_storeu_si128((__m128i *)(Temp + X), Min);
}
else
{
memset(Temp + X, , * sizeof(unsigned short));
}
}
for (int X = Block * BlockSize; X < Width; X++)
{
if (LinePS[X] == )
{
unsigned short Min0 = IM_Min(Second[X + ], Fourth[X + ]);
unsigned short Min1 = IM_Min(IM_Min(Second[X + ], Second[X + ]), IM_Min(Fourth[X + ], Fourth[X + ]));
unsigned short Min2 = IM_Min(IM_Min(IM_Min(First[X + ], First[X + ]), IM_Min(Second[X + ], Second[X + ])),
IM_Min(IM_Min(Fourth[X + ], Fourth[X + ]), IM_Min(Five[X + ], Five[X + ])));
Temp[X] = IM_Min(IM_Min(Min0 + ONE, Min1 + SQRT2), Min2 + SQRT5);
}
else
{
Temp[X] = ;
}
}
for (int X = ; X < Width; X++)
{
Third[X + ] = IM_Min(IM_Min(Third[X + ] + ONE, Third[X + ] + ONE), Temp[X]);
}
Third[] = Third[]; Third[] = Third[Width > ? : ]; // 镜像数据
Third[Width + ] = Third[Width + ]; Third[Width + ] = Third[Width > ? Width : ]; memcpy(Distance + Y * Width, Third + , Width * sizeof(unsigned short)); // 复制回去
} // 3. Repeat step 2, proceeding from right to left and bottom to top.
//
// 此处代码可自行添加
//
//
Status = IM_Normalize(Distance, Dest, Width, Height, true); FreeMemory:
if (Distance != NULL) free(Distance);
if (Buffer != NULL) free(Buffer);
return Status;
}
看起来代码量增加了不少，不过为了效率都值得啊，经过测试，SSE优化后的速度从之前的160ms提升至60ms左右，加速还是相当可观的。
上面代码有一处有的BUG，留待有兴趣研究的朋友自行查找。
那么我们来看看用EDM怎么实现腐蚀或者膨胀这一类操作。在上述代码中，最后的Distance数据中其实保存的就是某个点和最近的背景的距离，当然整个距离根据前面的代码是放大了100倍的，而且注意到距离的最小值必然为1（100），即在Distance数据中除了0之外的最小值为100。如果我们要进行腐蚀操作（即黑色部分向白色扩展），我们可以指定一个距离（即所谓的半径），在Distance中如果数据小于这个距离，则变成了背景，而只有大于整个距离的部分才会依旧是前景。注意到这里所谓的距离可以是小数（除以100后的结果），这就比传统的半径只能为整数的需求进了一步。而膨胀操作，只需要反色后在进行距离变换，然后在和腐蚀一样的操作即可。
我们来做个测试，如下图所示，很明显的EDM版本的腐蚀是各向同性的，原来是个圆，处理后还是个圆，而普通的算法，则逐渐的在想圆角矩形转变。
原图 半径为10的腐蚀（EDM版本） 半径为10的腐蚀（普通版本）
当然，普通版本的腐蚀也可以实现类似EDM那个效果的，比如我们使用真正圆形的半径的腐蚀（普通的半径的意义是都是矩形半径），但是一个核心的问题是随着半径的增加，圆形半径的计算量会成倍上升，虽然圆形半径也有快速算法，但是他始终无法做到矩形矩形半径那种O(1)算法的。而EDM算法确是和半径无关的。
EDM还有很多其他的用处，我们可以在ImageJ的代码里找到其更多的功能，这部分有待于读者自己去研究。
Demo下载地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar
