形态学轮廓提取

1. 流程分析
   
   原图像imgA：
   
   噪声滤除图像imgB:
   
   腐蚀处理图像imgC:
   
   相减操作图像imgD:
   
   二值化处理结果imgE：
   

2. 原理分析：
   2.1膨胀：求像素的局部最大值
   将图像（或图像的一部分区域，我们称之为A）与核（我们称之为B）进行如下卷积操作:
   
   2.2腐蚀：求像素的局部最小值
   将图像（或图像的一部分区域，我们称之为A）与核（我们称之为B）进行如下卷积操作:
   
   （膨胀腐蚀原理详细可参见博客：http://blog.sina.com.cn/s/blog_6f57a7150100ooin.html 非常感谢博主的分享。）
   2.3噪声滤除：
   开运算： A∘ B=(A⊝B)⊕B
   闭运算： A·B=(A⊕B)⊝B
   噪声滤除： {[(A⊝B)⊕B]⊕B}⊝B=（A∘ B）·B
   
   在上图中，（a）是原图像，外部有噪声块，内部有噪声孔，（b）为结构元素，尺寸大于所有噪声块和噪声孔，（c）是用（b）去腐蚀（a）的结果，可见，外部的噪声块被去除；之后用（b）去膨胀(c)两次，得到（e）,此时已经去除了内部的噪声孔，在进行一次腐蚀操作，得到和原图一样大小的去噪图像(f)。

3. 实验对比：
   3.1图像格式对比：
   参数设置：不进行噪声滤除，腐蚀2次，十字核，尺寸3*3
   二值图像;
   
   灰度图像：
   
   RGB图像：
   

   3.2核样式对比：
   参数设置：RGB图像，不进行”噪声滤除”，腐蚀1次，尺寸5*5，二值化阈值100.
   

   3.3核尺寸对比：
   参数设置：RGB图像，不进行”噪声滤除”，十字型核样式，腐蚀1次，二值化阈值100
   

   3.4腐蚀次数对比：
   参数设置：RGB图像，不进行”噪声滤除”，十字型核样式，尺寸5*5，二值化阈值100
   

   3.5噪声滤除强度对比：
   参数设置： RGB图像，十字型，腐蚀1次，核样式，尺寸5*5，二值化阈值100
   

   3.6对比分析：
   图像格式：RGB格式提取轮廓效果较好，灰度和二值图像信息相对有限。
   核样式：十字型效果最佳，矩形和圆形容易造成噪点。
   核尺寸：尺寸过小，轮廓提取不完整；尺寸过大，容易出现噪点。
   腐蚀次数：腐蚀次数过小，轮廓提取不完整；次数过多，容易出现噪点，且轮廓边缘强度过大。
   噪声滤除：影响把图像轮廓深度信息提取。（标红区域）

4. 代码：

morphology.h

#pragma once
#include "stdafx.h"
#include "opencv2/highgui/highgui.hpp"
#include <opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
using namespace std;

Mat getKernelMatrix(int shape, Size ksize);//内核矩阵生成
void erodeFun(Mat &src, Mat &dst, Mat kernel, int iterations);//腐蚀操作
void dilateFun(Mat &src, Mat &dst, Mat kernel, int iterations);//膨胀操作
void openFun(Mat &src, Mat &dst, Mat kernel, int iterations);//开操作
void closeFun(Mat &src, Mat &dst, Mat kernel, int iterations);//闭操作
void subtractFun(Mat src, Mat erode_ouput, Mat &externalGradientImg, int threshold);//减操作

morphology.cpp

#include "stdafx.h"
#include "morphology.h"
/*
函数： 获取内核矩阵
参数： shape内核样式（矩形: MORPH_RECT 交叉形十字: MORPH_CROSS 椭圆形: MORPH_ELLIPSE）
        ksize内核大小
*/
Mat getKernelMatrix(int shape, Size ksize) {
    Point center((int)ksize.width / 2, (int)ksize.height / 2);
    int i, j;
    int r = 0, c = 0;//圆参数
    double helperR = 0;//圆参数

    if (ksize == Size(1, 1))//单位内核都是矩形内核
        shape = MORPH_RECT;

    if (shape == MORPH_ELLIPSE)//椭圆形内核参数处理
    {
        r = ksize.height / 2;
        c = ksize.width / 2;
        helperR = r ? 1. / ((double)r*r) : 0;
    }

    Mat kernel(ksize, CV_8U);
    //行遍历
    for (i = 0; i < ksize.height; i++)
    {
        uchar* ptr = kernel.data + i*kernel.step;
        int j1 = 0, j2 = 0;

        //确定每一行1的区间[j1,j2]
        if (shape == MORPH_RECT || (shape == MORPH_CROSS && i == center.y))
            j2 = ksize.width;
        else if (shape == MORPH_CROSS)
            j1 = center.x, j2 = j1 + 1;
        else
        {
            int dy = i - r;
            if (std::abs(dy) <= r)//圆
            {
                int dx = saturate_cast<int>(c*sqrt((r*r - dy*dy)*helperR));
                j1 = std::max(c - dx, 0);
                j2 = std::min(c + dx + 1, ksize.width);
            }
        }
        //赋值操作
        for (j = 0; j < j1; j++)
            ptr[j] = 0;
        for (; j < j2; j++)
            ptr[j] = 1;
        for (; j < ksize.width; j++)
            ptr[j] = 0;
    }

    return kernel;
}

/*
函数： 腐蚀操作
参数： src:原图像
        dst:结果图像
        kernel：核矩阵
        iterations：操作次数
*/
void erodeFun(Mat &src, Mat &dst, Mat kernel, int iterations) {
    //参数判定
    if (iterations == 0) { return; }
    if (kernel.cols == 1 && kernel.rows == 1){ return; }


    Mat_<uchar> kernelMatrix = kernel;
    int HALFKERLEN = kernelMatrix.cols / 2;//正方形
    int KERLEN = kernelMatrix.cols;
    Mat pre = src.clone();
    Mat next = src.clone();

    for (int iterator = 0; iterator < iterations; iterator++) {
   //腐蚀次数
        for (int i = 0; i<pre.rows; i++)
        {
            for (int j = 0; j < pre.cols; j++)
            {
                //处理每一个像素值的每一个通道
                Vec3b minPixel = pre.at<Vec3b>(i, j);
                //遍历核
                for (int ki = 0; ki < KERLEN; ki++) {
                    for (int kj = 0; kj < KERLEN; kj++) {
                        if ((int)kernelMatrix.at<uchar>(ki, kj) == 1) {
   //核为1，处理
                            int pi = ki + i - HALFKERLEN;
                            int pj = kj + j - HALFKERLEN;
                            if (pi >= 0 && pi < pre.rows && pj >= 0 && pj < pre.cols)//防止超限
                            {
                                if (pre.at<Vec3b>(pi, pj)[0] < minPixel[0]) {
                                    minPixel[0] = pre.at<Vec3b>(pi, pj)[0];
                                }
                                if (pre.at<Vec3b>(pi, pj)[1] < minPixel[1]) {
                                    minPixel[1] = pre.at<Vec3b>(pi, pj)[1];
                                }
                                if (pre.at<Vec3b>(pi, pj)[2] < minPixel[2]) {
                                    minPixel[2] = pre.at<Vec3b>(pi, pj)[2];
                                }
                            }
                        }
                    }
                }
                next.at<Vec3b>(i, j) = minPixel;
            }
        }
        pre = next.clone();//更新
    }
    dst = next.clone();//更新
}


/*
函数： 膨胀操作
参数： src:原图像
        dst:结果图像
        kernel：核矩阵
        iterations：操作次数
*/
void dilateFun(Mat &src, Mat &dst, Mat kernel, int iterations) {

    Mat_<uchar> kernelMatrix = kernel;
    int HALFKERLEN = kernelMatrix.cols / 2;//正方形
    int KERLEN = kernelMatrix.cols;
    Mat pre = src.clone();
    Mat next = src.clone();

    for (int iterator = 0; iterator < iterations; iterator++) {
   //膨胀次数
        for (int i = 0; i<pre.rows; i++)
        {
            for (int j = 0; j < pre.cols; j++)
            {
                //处理每一个像素值的每一个通道
                Vec3b maxPixel = pre.at<Vec3b>(i, j);
                //遍历核
                for (int ki = 0; ki < KERLEN; ki++) {
                    for (int kj = 0; kj < KERLEN; kj++) {
                        if ((int)kernelMatrix.at<uchar>(ki, kj) == 1) {
   //核为1，处理
                            int pi = ki + i - HALFKERLEN;
                            int pj = kj + j - HALFKERLEN;
                            if (pi >= 0 && pi < pre.rows && pj >= 0 && pj < pre.cols)//不超限
                            {
                                if (pre.at<Vec3b>(pi, pj)[0] > maxPixel[0]) {
                                    maxPixel[0] = pre.at<Vec3b>(pi, pj)[0];
                                }
                                if (pre.at<Vec3b>(pi, pj)[1] > maxPixel[1]) {
                                    maxPixel[1] = pre.at<Vec3b>(pi, pj)[1];
                                }
                                if (pre.at<Vec3b>(pi, pj)[2] > maxPixel[2]) {
                                    maxPixel[2] = pre.at<Vec3b>(pi, pj)[2];
                                }
                            }
                        }
                    }
                }
                next.at<Vec3b>(i, j) = maxPixel;
            }
        }
        pre = next.clone();
    }
    dst = next.clone();
}


/*
函数： 开操作
参数： src:原图像
        dst:结果图像
        kernel：核矩阵
        iterations：操作次数
*/
void openFun(Mat &src, Mat &dst, Mat kernel, int iterations) {
    Mat erode = src.clone();
    erodeFun(src, erode, kernel, iterations);   //腐蚀
    dilateFun(erode, dst, kernel, iterations);//膨胀
}

/*
函数： 闭操作
参数： src:原图像
        dst:结果图像
        kernel：核矩阵
        iterations：操作次数
*/
void closeFun(Mat &src, Mat &dst, Mat kernel, int iterations) {
    if (iterations == 0) { return; }
    Mat dilate = src.clone();
    dilateFun(src, dilate, kernel, iterations); //腐蚀
    erodeFun(dilate, dst, kernel, iterations);//膨胀
}


void subtractFun(Mat src, Mat erode_ouput, Mat &externalGradientImg, int threshold) {

    for (int i = 0; i < src.rows; i++) {
        for (int j = 0; j < src.cols; j++) {
            /*cout << (int)src.at<Vec3b>(i, j)[0] << " ";
            cout << (int)src.at<Vec3b>(i, j)[1] << " ";
            cout << (int)src.at<Vec3b>(i, j)[2] << " " << endl;*/
            if (src.at<Vec3b>(i, j)[0] - erode_ouput.at<Vec3b>(i, j)[0] > threshold ||
                src.at<Vec3b>(i, j)[1] - erode_ouput.at<Vec3b>(i, j)[1] > threshold ||
                src.at<Vec3b>(i, j)[2] - erode_ouput.at<Vec3b>(i, j)[2] > threshold) {
                externalGradientImg.at<Vec3b>(i, j)[0] = 255;
                externalGradientImg.at<Vec3b>(i, j)[1] = 255;
                externalGradientImg.at<Vec3b>(i, j)[2] = 255;

            }
            else {
                externalGradientImg.at<Vec3b>(i, j)[0] = 0;
                externalGradientImg.at<Vec3b>(i, j)[1] = 0;
                externalGradientImg.at<Vec3b>(i, j)[2] = 0;
            }

        }
    }
}

main.cpp

#include "stdafx.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <stdlib.h>
#include <stdio.h>
#include "morphology.h"

using namespace cv;
/// 全局变量
Mat src, outLineResult;

int kernel_type = 0;//核样式
int const kernel_type_max = 2;
int kernel_size = 0;//核尺寸
int const kernel_size_max= 9;
int erode_times = 0;//腐蚀次数
int const erode_times_max = 5;
int removedetails_times = 0;//去除细节次数
int const removedetails_times_max = 2;

void outLine(int, void*);

int main()
{
    src = imread("test.png");

    if (!src.data)
    {
        return -1;
    }

    /// 创建显示窗口
    namedWindow("Erosion Demo", CV_WINDOW_AUTOSIZE);
    namedWindow("Dilation Demo", CV_WINDOW_AUTOSIZE);
    cvMoveWindow("Dilation Demo", src.cols, 0);

    createTrackbar("核样式:", "Erosion Demo",
        &kernel_type, kernel_type_max,
        outLine);
    createTrackbar("核尺寸:", "Erosion Demo",
        &kernel_size, kernel_size_max,
        outLine);
    createTrackbar("腐蚀次数:", "Erosion Demo",
        &erode_times, erode_times_max,
        outLine);
    createTrackbar("去噪强度:", "Erosion Demo",
        &removedetails_times, removedetails_times_max,
        outLine);
    /// Default start
    outLine(0, 0);


    waitKey(0);
    return 0;
}

/**  @function Erosion  */
void outLine(int, void*)
{
    int type;
    string tpyeStr = "";
    string removeDetailType = "";
    if (kernel_type == 0) { type = MORPH_RECT; tpyeStr = "正方形"; }
    else if (kernel_type == 1) { type = MORPH_CROSS; tpyeStr = "十字形"; }
    else if (kernel_type == 2) { type = MORPH_ELLIPSE; tpyeStr = "圆形"; }

    //输出
    cout << endl<<"----------------------" << endl;
    cout << "  核样式："<< tpyeStr;
    cout << "  核尺寸：" << kernel_size;
    cout << "  轮廓强度：" << erode_times;
    cout << "  细节去除样式：" << removeDetailType;
    cout << "  细节去除强度：" << removedetails_times;

    Mat element = getKernelMatrix(type, Size(kernel_size, kernel_size));
    Mat erode_ouput,closeResult,openResult;
    erode_ouput = src.clone();
    outLineResult = src.clone();
    openResult = src.clone();
    closeResult = src.clone();
    //噪声滤除
    if (removedetails_times == 1) {
        closeFun(src, closeResult, element, 1);
        openFun(closeResult, openResult, element, 1);
    }
    else if (removedetails_times == 2) {
        openFun(src, openResult, element, 1);//去噪
        closeFun(openResult, closeResult, element, 1);
        openFun(closeResult, openResult, element, 1);//去噪
        closeFun(openResult, closeResult, element, 1);  
    }

    // 腐蚀操作
    erodeFun(closeResult, erode_ouput, element , erode_times);

    //减操作
    subtractFun(closeResult, erode_ouput, outLineResult,50);

    imshow("Erosion Demo", outLineResult);
    cout << "  完成";

}

人脸边缘轮廓提取（下巴、嘴唇）

使用上述方法对人脸进行轮廓提取，调整参数，得到相对最佳的效果，但是会出现下巴和嘴唇提取失败的现象，所以针对这些边缘信息，要自行提取。

1. 流程分析：
   

2. 嘴唇定位&内外轮廓及中心点确认：

   使用人脸特征提取和模式匹配方法去定位嘴唇并确认中心点以及内外轮廓。
   

3. 自适应算法确定轮廓点
   

Step1:确定路径L(i)的起始点start(i)和终止点end(i) 。
Step2:从start(i)向end(i)做步长为1的像素点搜索，每次搜索计算像素梯度值，当梯度值大于设定的阈值时，停止该次搜索，并保存该点为轮廓点。
Step3:确定是否完成搜索：如完成，进行下一步轮廓提取。如果未完成，更新搜索路径L(i+1),转至Step1.

提取结果：

改变搜索方向，即可对嘴唇进行轮廓点的提取：