二值图像距离图生成算法

算法主要步骤

初始化

• 将图像分为内部点和非内部点（内部点本身为1且4领域皆为1）
• 初始化距离图像，内部点值为INF，非内部点为0

更新距离

- 由上到下，由左至右迭代计算，伪代码如下:

def forward_conv(image):
    for i in range(hight):
        for j in range(width):
            for pixel in Neighbour:
                dist = template[pixel]+dist[pixel]
            find the minimun dist
            dist[i][j] = min(dist[i][j],dist)

• 后向更新距离图

代码

代码是基于ITK写的，将就看吧

#include <iostream>
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImage.h"
#include <itkImageRegionIteratorWithIndex.h>

using namespace std;

typedef     itk::Image< unsigned char, 2 >      ImageType;
typedef     itk::Image< float, 2 >              FltImageType;

#define INF 1e8

FltImageType::Pointer InitEDT(ImageType::Pointer image)
{
    FltImageType::Pointer edt = FltImageType::New();
    edt->SetRegions(image->GetBufferedRegion());
    edt->Allocate();

    itk::ImageRegionIteratorWithIndex<FltImageType> it1(edt, edt->GetBufferedRegion());
    itk::ImageRegionIteratorWithIndex<ImageType>    it2(image, image->GetBufferedRegion());
    it1.GoToBegin();
    it2.GoToBegin();
    int N4[4][2] = { {-1,0},{1,0},{0,-1},{0,1} };
    while (!it1.IsAtEnd())
    {
        if (it2.Get() != 0)
        {
            ImageType::IndexType idx = it2.GetIndex();
            int cnt = 0;
            for (int i = 0; i < 4; i++)
            {
                ImageType::IndexType tmp;
                tmp[0] = idx[0] + N4[i][0];
                tmp[1] = idx[1] + N4[i][1];
                if (image->GetBufferedRegion().IsInside(tmp) && image->GetPixel(tmp) != 0)
                    cnt++;
            }
            if (cnt == 4)
                it1.Set(INF);
            else
                it1.Set(0);
        }
        else
            it1.Set(0);
        ++it1;
        ++it2;
    }
    return edt;
}

void ForwardConv(FltImageType::Pointer edt)
{
    int N4[4][2] = { { 1,-1 },{ 0,-1 },{ -1,-1 },{ -1,0 } };
    int D4[4] = { 4,3,4,3 };
    for (int i = 0; i < edt->GetBufferedRegion().GetSize()[0]; i++)
    {
        for (int j = 0; j < edt->GetBufferedRegion().GetSize()[1]; j++)
        {
            FltImageType::IndexType idx;
            idx[0] = i;
            idx[1] = j;
            if (edt->GetBufferedRegion().IsInside(idx) && edt->GetPixel(idx) > 0)
            {
                float val = edt->GetPixel(idx);
                for (int k = 0; k < 4; k++)
                {
                    FltImageType::IndexType nei;
                    nei[0] = idx[0] + N4[k][0];
                    nei[1] = idx[1] + N4[k][1];
                    if (edt->GetBufferedRegion().IsInside(nei))
                        val = (edt->GetPixel(nei) + D4[k]) < val ? (edt->GetPixel(nei) + D4[k]) : val;
                }
                edt->SetPixel(idx, val);
            }
        }
    }
}

void BackwardConv(FltImageType::Pointer edt)
{
    int N4[4][2] = { { -1,1 },{ 0,1 },{ 1,1 }, { 1,0 } };
    int D4[4] = { 4,3,4,3 };
    for (int i = edt->GetBufferedRegion().GetSize()[0] - 1; i >= 0; i--)
    {
        for (int j = edt->GetBufferedRegion().GetSize()[1] - 1; j >= 0; j--)
        {
            FltImageType::IndexType idx;
            idx[0] = i;
            idx[1] = j;
            if (edt->GetBufferedRegion().IsInside(idx) && edt->GetPixel(idx) > 0)
            {
                float val = edt->GetPixel(idx);
                for (int k = 0; k < 4; k++)
                {
                    FltImageType::IndexType nei;
                    nei[0] = idx[0] + N4[k][0];
                    nei[1] = idx[1] + N4[k][1];
                    if (edt->GetBufferedRegion().IsInside(nei))
                        val = (edt->GetPixel(nei) + D4[k]) < val ? (edt->GetPixel(nei) + D4[k]) : val;
                }
                edt->SetPixel(idx, val);
            }
        }
    }
}

void main()
{
    typedef     itk::ImageFileReader< ImageType >   ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName("1.mha");
    reader->Update();

    FltImageType::Pointer edt = InitEDT(reader->GetOutput());
    ForwardConv(edt);
    BackwardConv(edt);

    typedef itk::ImageFileWriter< FltImageType >    WriterType;
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput(edt);
    writer->SetFileName("2.mha");
    writer->Update();
    system("pause");
}

算法效果

输入：

输出：

参考

[1].https://blog.csdn.net/trent1985/article/details/18081761