显变化

代码实现：

Mat RGB2GRAY(const Mat &src){
    if(!src.data || src.channels()!=3){
        exit(1);
    }
    int row = src.rows;
    int col = src.cols;
    Mat gray = Mat(row, col, CV_8UC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            gray.at<uchar>(i, j) = (uchar)(0.114 * src.at<Vec3b>(i, j)[0] + 0.587 * src.at<Vec3b>(i, j)[1] + 0.299 * src.at<Vec3b>(i, j)[2]);
        }
    }
    return gray;
}

void SobelGradDirction(Mat &src, Mat &sobelX, Mat &sobelY){
    int row = src.rows;
    int col = src.cols;
    sobelX = Mat::zeros(src.size(), CV_32SC1);
    sobelY = Mat::zeros(src.size(), CV_32SC1);

    for(int i = 0; i < row-1; i++){
        for(int j = 0; j < col-1; j++){
            double gradY = src.at<uchar>(i+1, j-1) + src.at<uchar>(i+1, j) * 2 + src.at<uchar>(i+1, j+1) - src.at<uchar>(i-1, j-1) - src.at<uchar>(i-1, j) * 2 - src.at<uchar>(i-1, j+1);
            sobelY.at<float>(i, j) = abs(gradY);
            double gradX = src.at<uchar>(i-1, j+1) + src.at<uchar>(i, j+1) * 2 + src.at<uchar>(i+1, j+1) - src.at<uchar>(i-1, j-1) - src.at<uchar>(i, j-1) * 2 - src.at<uchar>(i+1, j-1);
            sobelX.at<float>(i, j) = abs(gradX);
        }
    }
    //将梯度数组转换为8位无符号整形
    convertScaleAbs(sobelX, sobelX);
    convertScaleAbs(sobelY, sobelY);
}

Mat SobelXX(const Mat src){
    int row = src.rows;
    int col = src.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src.at<uchar>(i, j) * src.at<uchar>(i, j);
        }
    }
    return dst;
}

Mat SobelYY(const Mat src){
    int row = src.rows;
    int col = src.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src.at<uchar>(i, j) * src.at<uchar>(i, j);
        }
    }
    return dst;
}

Mat SobelXY(const Mat src1, const Mat &src2){
    int row = src1.rows;
    int col = src1.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            dst.at<float>(i, j) = src1.at<uchar>(i, j) * src2.at<uchar>(i, j);
        }
    }
    return dst;
}

void separateGaussianFilter(Mat_<float> &src, Mat_<float> &dst, int ksize, double sigma){
    CV_Assert(src.channels()==1 || src.channels() == 3); //只处理单通道或者三通道图像
    //生成一维的
    double *matrix = new double[ksize];
    double sum = 0;
    int origin = ksize / 2;
    for(int i = 0; i < ksize; i++){
        double g = exp(-(i-origin) * (i-origin) / (2 * sigma * sigma));
        sum += g;
        matrix[i] = g;
    }
    for(int i = 0; i < ksize; i++) matrix[i] /= sum;
    int border = ksize / 2;
    copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT);
    int channels = dst.channels();
    int rows = dst.rows - border;
    int cols = dst.cols - border;
    //水平方向
    for(int i = border; i < rows - border; i++){
        for(int j = border; j < cols - border; j++){
            float sum[3] = {0};
            for(int k = -border; k<=border; k++){
                if(channels == 1){
                    sum[0] += matrix[border + k] * dst.at<float>(i, j+k);
                }else if(channels == 3){
                    Vec3f rgb = dst.at<Vec3f>(i, j+k);
                    sum[0] += matrix[border+k] * rgb[0];
                    sum[1] += matrix[border+k] * rgb[1];
                    sum[2] += matrix[border+k] * rgb[2];
                }
            }
            for(int k = 0; k < channels; k++){
                if(sum[k] < 0) sum[k] = 0;
                else if(sum[k] > 255) sum[k] = 255;
            }
            if(channels == 1)
                dst.at<float>(i, j) = static_cast<float>(sum[0]);
            else if(channels == 3){
                Vec3f rgb = {static_cast<float>(sum[0]), static_cast<float>(sum[1]), static_cast<float>(sum[2])};
                dst.at<Vec3f>(i, j) = rgb;
            }
        }
    }
    //竖直方向
    for(int i = border; i < rows - border; i++){
        for(int j = border; j < cols - border; j++){
            float sum[3] = {0};
            for(int k = -border; k<=border; k++){
                if(channels == 1){
                    sum[0] += matrix[border + k] * dst.at<float>(i+k, j);
                }else if(channels == 3){
                    Vec3f rgb = dst.at<Vec3f>(i+k, j);
                    sum[0] += matrix[border+k] * rgb[0];
                    sum[1] += matrix[border+k] * rgb[1];
                    sum[2] += matrix[border+k] * rgb[2];
                }
            }
            for(int k = 0; k < channels; k++){
                if(sum[k] < 0) sum[k] = 0;
                else if(sum[k] > 255) sum[k] = 255;
            }
            if(channels == 1)
                dst.at<float>(i, j) = static_cast<float>(sum[0]);
            else if(channels == 3){
                Vec3f rgb = {static_cast<float>(sum[0]), static_cast<float>(sum[1]), static_cast<float>(sum[2])};
                dst.at<Vec3f>(i, j) = rgb;
            }
        }
    }
    delete [] matrix;
}

Mat harrisResponse(Mat_<float> GaussXX, Mat_<float> GaussYY, Mat_<float> GaussXY, float k){
    int row = GaussXX.rows;
    int col = GaussXX.cols;
    Mat_<float> dst(row, col, CV_32FC1);
    for(int i = 0; i < row; i++){
        for(int j = 0; j < col; j++){
            float a = GaussXX.at<float>(i, j);
            float b = GaussYY.at<float>(i, j);
            float c = GaussXY.at<float>(i, j);
            dst.at<float>(i, j) = a * b - c * c - k * (a + b) * (a + b);
        }
    }
    return dst;
}

int dir[8][2] = {0, 1, 0, -1, 1, 0, -1, 0, 1, 1, 1, -1, -1, 1, -1, -1};

void MaxLocalValue(Mat_<float>&resultData, Mat srcGray, Mat &resultImage, int kSize){
    int r = kSize / 2;
    resultImage = srcGray.clone();
    int row = resultImage.rows;
    int col = resultImage.cols;
    for(int i = r; i < row - r; i++){
        for(int j = r; j < col - r; j++){
            bool flag = true;
            for(int k = 0; k < 8; k++){
                int tx = i + dir[k][0];
                int ty = j + dir[k][1];
                if(resultData.at<float>(i, j) < resultData.at<float>(tx, ty)){
                    flag = false;
                }
            }
            if(flag && (int)resultData.at<float>(i, j) > 18000){
                circle(resultImage, Point(i, j), 5, Scalar(0, 0, 255), 2, 8, 0);
            }
        }
    }
}

效果图：

?

?1.3数学方法刻画角点特征

?

?

?

?

?

共可分为三种情况：

? ?（1） 图像中的直线：一个特征值大，另一个特征值小，λ1>λ2或λ2>λ1。自相关函数值在某一方向上大，在其他方向上小。

? ?（2）图像中的平面：两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。

? ?（3）图像中的角点：两个特征值都大，且近似相等，自相关函数在所有方向都增大。

?1.4角点响应函数

由于我们是通过M的两个特征值的大小对图像进行分类，所以，定义角点相应函数R：?

?

?其中k为经验常数，一般取k=0.04~0.06。增大k的值将减小响应值R，降低检测的灵敏性，减少被检测角点的数量；减小k值，将增大角点响应值R，增加被检测角点的数量。为了去除加权常数κ，我们通常使用商数detM/(traceM)2作为指示器。所以，上图可以转化为：

?

? ?R 只与M的特征值有关

? ? ? ? ? ? ? ?? 角点：R 为大数值正数（λ1和λ2都是很大的正数）

? ? ? ? ? ? ? ?? 边缘：R 为大数值负数

? ? ? ? ? ? ? ?? 平坦区：R 为小数值（如果λ1≈λ2≈0，该区域为空）

? ? ? ? ?在判断角点的时候，对角点响应函数R进行阈值处理：R > threshold，提取R的局部极大值。
?

1.5Harris角点检测步骤

Harris角点检测可以分为5个步骤：

（1）计算图像I(x,y)I(x,y)在xx和yy两个方向的梯度IxIx,IyIy;

（2）计算图像两个方向梯度的乘积;

?

（3）使用高斯函数对Ix2、Iy2、IxIy进行高斯加权(取σ=2,ksize=3)，计算中心点为(x,y)(x,y)的窗口W对应的矩阵M；

?

（4）计算每个像素点(x,y)处的(x,y)处的Harris响应值R；?

?（5）过滤大于某一阈值t的R值；

?1.6角点检测实例

（1）代码实现:

# -*- coding: utf-8 -*-
from pylab import *
from PIL import Image
from PCV.localdescriptors import harris
from scipy.ndimage import filters

def compute_harris_response(im, sigma=3):
    # 在一幅灰度图像中，对每个像素计算Harris角点检测器响应函数
    
    # 计算导数
    k = 0.04
    imx = zeros(im.shape)
    filters.gaussian_filter(im, (sigma, sigma), (0, 1), imx)
    imy = zeros(im.shape)
    filters.gaussian_filter(im, (sigma, sigma), (1, 0), imy)
    
    # 计算harris矩阵分量
    Wxx = filters.gaussian_filter(imx * imx, sigma)
    Wxy = filters.gaussian_filter(imx * imy, sigma)
    Wyy = filters.gaussian_filter(imy * imy, sigma)
    
    # 计算矩阵的特征值和迹
    Wdet = Wxx * Wyy - Wxy ** 2
    Wtr = Wxx + Wyy
    
    # 返回像素值为 Harris 响应函数值的一幅图像
    print(Wdet-(k*Wtr)) # 输出角点响应函数R
    return Wdet / Wtr  # 此处可消除参数k的影响
    
# 读入图像(读取图像到数组中)
im = array(Image.open('p3.jpg').convert('L'))  # 括号内：读取一幅图像，并将其转换成灰度图像

# 检测harris角点
harrisim = compute_harris_response(im)
                # compute_harris_response(im, sigma):
                #     在一幅灰度图像中，对每一个像素计算Harris角点检测器响应函数
                #     im:(数组图像）  sigma:标准差默认为3
                
# Harris响应函数
# print (harrisim)
harrisim1 = 255 - harrisim
figure()
gray()

#画出Harris响应图
subplot(141)
imshow(harrisim1)
print harrisim1.shape
axis('off')
axis('equal')

threshold = [0.01, 0.05, 0.1]
for i, thres in enumerate(threshold):
    # enumerate()函数用于将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标，一般用在for 循环当中
        filtered_coords = harris.get_harris_points(harrisim, 6, thres)
                                # get_harris_points(harrisim, min_dist=10, threshold=0.1):
                                #       从一幅Harris 响应图像harrisim中返回角点。min_dist 为分割角点和图像边界的最少像素数目（默认为10）。
                                #       threshold 为阀值，大于阀值的harris响应函数值被认为是可能的角点，并在harrism_t矩阵中相应的位置1,其余地方置0
        subplot(1, 4, i+2)
        imshow(im)
        print im.shape
        plot([p[1] for p in filtered_coords], [p[0] for p in filtered_coords], '*') # 标出角点
        axis('off')
        
show()

上面代码中重写compute_harris_response()函数，以便将角点响应函数R输出。

截图：

?（2）输入图片特征为纹理平坦丰富：

?

?1.7图像中寻找对应点

Harris 角点检测器仅仅能够检测出图像中的兴趣点，但是没有给出通过比较图像间的兴趣点来寻找匹配角点的方法。我们需要在每个点上加入描述子信息，并给出一 个比较这些描述子的方法。兴趣点描述子是分配给兴趣点的一个向量，描述该点附近的图像的表观信息。描述子越好，寻找到的对应点越好。我们用对应点或者点的对应来描述相同物体和场景点在不同图像上形成的像素点。使用归一化的互相关矩阵进行两幅图像的比较，代码如下：
?

def get_descriptors(image,filtered_coords,wid=5):
    """对于每个返回的点，返回点周围2*wid+1个像素的值（假设选取点的min_distance > wid）"""
    desc = []
    for coords in filtered_coords:
        patch = image[coords[0] - wid:coords[0]+wid+1,coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)
    return desc
 
def match(desc1,desc2,threshold=0.5):
    """对于第一幅图像中的每个角点描述子，使用归一化互相关，选取它在第二幅图像中的匹配角点"""
    n = len(desc1[0])
    # 点对的距离
    d = -np.ones((len(desc1), len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            # mean:求平均值，std:求标准差
            d1 = (desc1[i] - np.mean(desc1[i])) / np.std(desc1[i])
            d2 = (desc2[j] - np.mean(desc2[j])) / np.std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1)
            if ncc_value > threshold:  # 数值较高的距离代表两个点
                d[i,j] = ncc_value
    # np.argsort():将数组中的元素从小到大排列，提取其对应的index(索引)，然后输出
    # 数值较高的距离代表两个点能够更好的匹配，因此对距离取相反数
    ndx = np.argsort(-d)
    matchscores = ndx[:0]
    return matchscores
 
# 从第二幅图像向第一幅图像匹配，然后过滤掉在两种方法种不都是最好的匹配
def match_twosided(desc1,desc2,threshold=0.5):
    """两边对称版本的match()"""
    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)
    # np.where:找到n维数组中特定数值的索引
    ndx_12 = np.where(matches_12 >= 0)[0]
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1
    return matches_12
 
# 实现匹配点的可视化
def appendimages(im1,im2):
    """返回将两幅图像并排拼接成一幅图像"""
    # 选取具有最少行数的图像，然后填充足够的空行
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]
 
    if rows1 < rows2:
        im1 = np.concatenate((im1, np.zeros((rows2 - rows1, im1.shape[1]))), axis=0)
    elif rows1 > rows2:
        im2 = np.concatenate((im2, np.zeros((rows1 - rows2, im2.shape[1]))), axis=0)
    # 如果无此两种情况，则他们行数相同，无需填充
    return np.concatenate((im1,im2),axis=1)
 
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
    """显示一幅带有连接匹配之间连线的图片"""
    """
    im1,im2:数组图像，locs1.locs2:特征位置，matchscores(match())的输出，
    show_below(如果图像应该显示在匹配的下方)
    """
    im3 = appendimages(im1,im2)
    if show_below:
        im3 = np.vstack((im3, im3))
    plt.imshow(im3)
    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0:
            plt.plot([locs1[i][1],locs2[m][1]+cols1],[locs1[i][0],locs2[m][0]],'c')
    plt.axis('off')

测试代码：

# -*- coding: utf-8 -*-
"""
@author: RRJ
@software: PyCharm
@file: harrisdemo3.py
@time: 2022/3/27 15:19
"""
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from Localdescriptors import harris
from tools.imtools import imresize
 
# 读入图像
im1 = np.array(Image.open("../JMU_img/JMU_ZS1.jpg").convert("L"))
im2 = np.array(Image.open("../JMU_img/JMU_ZS2.jpg").convert("L"))
 
# resize加快匹配速度
im1 = imresize(im1, (im1.shape[1]//2, im1.shape[0]//2))
im2 = imresize(im2, (im2.shape[1]//2, im2.shape[0]//2))
 
wid = 5
harrisim = harris.compute_harris_response(im1, 5)
filtered_coords1 = harris.get_harris_points(harrisim, wid+1)
d1 = harris.get_descriptors(im1, filtered_coords1, wid)
 
harrisim = harris.compute_harris_response(im2, 5)
filtered_coords2 = harris.get_harris_points(harrisim, wid+1)
d2 = harris.get_descriptors(im2, filtered_coords2, wid)
 
print('starting matching')
matches = harris.match_twosided(d1, d2)
 
plt.figure()
plt.gray()
harris.plot_matches(im1, im2, filtered_coords1, filtered_coords2, matches)
plt.show()

运行截图：

剩下的明天再写！！！！?