[转载]距离变换算法

出处：http://blog.csdn.net/lichengyu/article/details/8009246

《图像处理分析与机器视觉》

 

计算全局图像中各个像素点对子图像的距离。

AL	AL
AL	P
AL

 

                          　　　　　　　　

Mask  1

 

	BR
P	BR
BR	BR

 

                       

Mask 2

 

1. 将图像进行二值化，子图像值为0，背景为255；

2. 利用Mask 1从左向右，从上到下扫描，p点是当前像素点，q点是Mask 1中AL邻域中的点，D()为距离计算，包括棋盘距离、城市距离和欧式距离。F(p)为p点的像素值，计算

F(p) = min( F(p),  F(q)+D(p,q) ), 其中，q属于AL.

3. 再利用Mask 2从右向左，从下向上扫描，计算

F(p) = min( F(p),  F(q)+D(p,q) ), 其中，q属于BR

4. F(p) 则为距离变换后的图像。

 

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "math.h"
#include <stdio.h>


using namespace cv;

Mat& DistTran(Mat& I);
Mat& NormImage(Mat& I);

static float Round(float f)
{
    return ( ceil(f)-f > f-floor(f) ) ? floor(f) : ceil(f);
}

int ChessBoardDist(int x1, int y1, int x2, int y2)
{
    return (abs(x1-x2) > abs(y1-y2)) ? abs(x1-x2) : abs(y1-y2);
}

int CityBlockDist(int x1, int y1, int x2, int y2)
{
    return ( abs(x1-x2) + abs(y1-y2) );
}

float EuclideanDist(int x1, int y1, int x2, int y2)
{
    return sqrt( (float)((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)) );
}

int MyMin(int x, int y)
{
    return (x < y) ? x : y;
}

float MyMin(float x, float y)
{
    return (x < y) ? x : y;
}

/**
 * @function main
 */
int main( int argc, char** argv )
{
  char* imageName = "test_wushuang.bmp";

  Mat image;
  image = imread(imageName,1);
  if(!image.data)
  {
      printf("No image data\n");
  }

  Mat gray_image;
  cvtColor( image, gray_image, CV_RGB2GRAY );

  
  DistTran( gray_image );
  NormImage( gray_image );

  
  imwrite("EuclideanDist_wushuang.bmp", gray_image);

  namedWindow( imageName, CV_WINDOW_AUTOSIZE );
  namedWindow( "Gray image", CV_WINDOW_AUTOSIZE );

  imshow( imageName, image );
  imshow( "Gray image", gray_image );

  waitKey(0);


  
  return 0;
}


Mat& DistTran(Mat& I)
{
    // accept only char type matrices
    CV_Assert(I.depth() != sizeof(uchar));
    int channels = I.channels();
    int nRows = I.rows * channels;
    int nCols = I.cols;
    //if (I.isContinuous())
    //{
    //    nCols *= nRows;
    //    nRows = 1;
    //}

    int i,j;
    uchar* p;
    uchar* q;
    //int min = 0;
    //int dis = 0;
    float fMin = 0.0;
    float fDis = 0.0;

    // pass throuth from top to bottom, left to right
    for( i = 1; i < nRows-1; ++i)
    {
        p = I.ptr<uchar>(i);
        for ( j = 1; j < nCols; ++j)
        {
            /*q = I.ptr<uchar>(i-1);
            dis = CityBlockDist(i, j, i-1, j-1);
            min = MyMin( p[j], dis+q[j-1] );

            dis = CityBlockDist(i, j, i-1, j);
            min = MyMin( min, dis+q[j] );

            q = I.ptr<uchar>(i);
            dis = CityBlockDist(i, j, i, j-1);
            min = MyMin( min, dis+q[j-1] );

            q = I.ptr<uchar>(i+1);
            dis = CityBlockDist(i, j, i+1, j-1);
            min = MyMin( min, dis+q[j-1] );

            p[j] = min;*/

            q = I.ptr<uchar>(i-1);
            fDis = EuclideanDist(i, j, i-1, j-1);
            fMin = MyMin( (float)p[j], fDis+q[j-1] );

            fDis = EuclideanDist(i, j, i-1, j);
            fMin = MyMin( fMin, fDis+q[j] );

            q = I.ptr<uchar>(i);
            fDis = EuclideanDist(i, j, i, j-1);
            fMin = MyMin( fMin, fDis+q[j-1] );

            q = I.ptr<uchar>(i+1);
            fDis = EuclideanDist(i, j, i+1, j-1);
            fMin = MyMin( fMin, fDis+q[j-1] );

            p[j] = (uchar)Round(fMin);
        }
    }

    // pass throuth from bottom to top, right to left
    for( i = nRows-2; i > 0; i-- )
    {
        p = I.ptr<uchar>(i);
        for( j = nCols-1; j >= 0; j-- )
        {
            /*q = I.ptr<uchar>(i+1);
            dis = CityBlockDist(i, j, i+1, j);
            min = MyMin( p[j], dis+q[j] );

            dis = CityBlockDist(i, j, i+1, j+1);
            min = MyMin( min, dis+q[j+1] );

            q = I.ptr<uchar>(i);
            dis = CityBlockDist(i, j, i, j+1);
            min = MyMin( min, dis+q[j+1] );

            q = I.ptr<uchar>(i-1);
            dis = CityBlockDist(i, j, i-1, j+1);
            min = MyMin( min, dis+q[j+1] );

            p[j] = min;*/

            q = I.ptr<uchar>(i+1);
            fDis = EuclideanDist(i, j, i+1, j);
            fMin = MyMin( (float)p[j], fDis+q[j] );

            fDis = EuclideanDist(i, j, i+1, j+1);
            fMin = MyMin( fMin, fDis+q[j+1] );

            q = I.ptr<uchar>(i);
            fDis = EuclideanDist(i, j, i, j+1);
            fMin = MyMin( fMin, fDis+q[j+1] );

            q = I.ptr<uchar>(i-1);
            fDis = EuclideanDist(i, j, i-1, j+1);
            fMin = MyMin( fMin, fDis+q[j+1] );

            p[j] = (uchar)Round(fMin);
        }
    }

    return I;
}


Mat& NormImage(Mat& I)
{
    // accept only char type matrices
    CV_Assert(I.depth() != sizeof(uchar));
    int channels = I.channels();
    int nRows = I.rows * channels;
    int nCols = I.cols;
    //if (I.isContinuous())
    //{
    //    nCols *= nRows;
    //    nRows = 1;
    //}

    int i,j;
    uchar* p;
    int min = 256;
    int max = -1;

    // Do not count the outer boundary
    for( i = 1; i < nRows-1; ++i)
    {
        p = I.ptr<uchar>(i);
        for ( j = 1; j < nCols-1; ++j)
        {
            if( min > p[j] )  min = p[j];
            if( max < p[j] )  max = p[j];
        }
    }

    for( i = 1; i < nRows-1; ++i)
    {
        p = I.ptr<uchar>(i);
        for ( j = 1; j < nCols-1; ++j)
        {
            p[j] = (p[j] - min) * 255 / (max - min);
        }
    }

    return I;
}