下载文档原格式
PyCharm是一种流行的集成开发环境（IDE），用于编写Python代码。
SIFT（尺度不变特征转换）是一种计算机视觉算法，用于检测和描述图像中的局部特征。
在本文中，我们将介绍如何使用PyCharm编写SIFT特征提取的代码。
SIFT特征提取是一种在计算机视觉领域广泛应用的技术，它可以检测图像中的角点和边缘，并生成具有旋转和尺度不变性的特征描述子。
这使得SIFT特征在图像匹配、目标识别和三维重构等领域具有很高的实用价值。
在PyCharm中编写SIFT特征提取的代码需要以下步骤：1. 安装OpenCV库：OpenCV是一个开源的计算机视觉库，它包含了许多用于图像处理和计算机视觉领域的函数和工具。
在PyCharm中，我们可以使用pip命令来安装OpenCV库。
2. 导入OpenCV库：在PyCharm中，我们使用import语句来导入OpenCV库，以便在代码中调用OpenCV中的函数和类。
3. 加载图像：使用OpenCV库中的imread函数来加载要进行特征提取的图像，将其存储在一个变量中以备后续使用。
4. 创建SIFT对象：通过调用OpenCV库中的SIFT_create函数来创建一个SIFT对象，以便对图像进行特征提取操作。
5. 提取特征点和描述子：使用SIFT对象的detectAndCompute函数来提取图像中的关键点和对应的描述子，这些描述子将用于后续的特征匹配和识别操作。
6. 显示特征点：通过调用OpenCV库中的drawKeypoints函数来在图像上绘制出提取出的特征点，以便可视化和检查特征提取的效果。
7. 保存结果图像：使用OpenCV库中的imwrite函数将包含特征点的图像保存到指定的文件路径。
通过上述步骤，我们可以在PyCharm中编写出SIFT特征提取的代码，并在图像上成功提取出关键点和描述子。
这样的代码可以被应用到图像匹配、目标识别和其他计算机视觉任务中，为实际应用提供了便利。
sift算法代码
以下是SIFT算法的伪代码：
1. 尺度空间的构建:
a. 对图像进行高斯滤波，生成不同尺度的图像金字塔。
b. 使用拉普拉斯金字塔计算图像的尺度空间。
2. 关键点检测:
a. 对每个尺度空间进行极值点检测，找到图像中的候选关键点。
b. 对候选关键点进行精确定位，排除不稳定的关键点。
3. 方向计算:
a. 对每个关键点邻域内的像素计算梯度幅值和方向。
b. 构建梯度方向直方图，找到主要方向。
4. 描述子生成:
a. 将关键点邻域划分为若干个子区域。
b. 在每个子区域内计算梯度方向直方图，得到局部特征描述子。
5. 关键点匹配:
a. 使用特征描述子之间的距离进行关键点匹配。
b. 使用最近邻筛选和次近邻比率筛选进行匹配点的筛选。
其中，步骤2-4都是针对单个关键点进行操作的，可以进行并
行计算。
具体的代码实现可以参考开源的SIFT库，如OpenCV中的SIFT实现。
一、介绍SIFT算法SIFT（Scale-Invariant Feature Transform）算法是一种用于图像处理和计算机视觉领域的特征提取算法，由David Lowe在1999年提出。
SIFT算法具有旋转、尺度、光照等方面的不变性，能够对图像进行稳健的特征点提取，被广泛应用于物体识别、图像匹配、图像拼接、三维重建等领域。
二、SIFT算法原理SIFT算法的主要原理包括尺度空间极值点检测、关键点定位、关键点方向确定、关键点描述等步骤。
其中，尺度空间极值点检测通过高斯差分金字塔来检测图像中的极值点，关键点定位则利用DoG响应函数进行关键点细化，关键点方向确定和关键点描述部分则通过梯度方向直方图和关键点周围区域的梯度幅度信息来完成。
三、使用Matlab实现SIFT算法在Matlab中实现SIFT算法，需要对SIFT算法的每个步骤进行详细的编程和调试。
需要编写代码进行图像的高斯金字塔和高斯差分金字塔的构建，计算尺度空间极值点，并进行关键点定位。
需要实现关键点的方向确定和描述子生成的算法。
将所有步骤整合在一起，完成SIFT算法的整体实现。
四、SIFT算法复杂代码的编写SIFT算法涉及到的步骤较多，需要编写复杂的代码来实现。
在编写SIFT算法的Matlab代码时，需要考虑到算法的高效性、可扩展性和稳定性。
具体来说，需要注意以下几点：1. 高斯差分金字塔和高斯金字塔的构建：在构建高斯差分金字塔时，需要编写代码实现图像的高斯滤波和图像的降采样操作，以得到不同尺度空间的图像。
还需要实现高斯差分金字塔的构建，以检测图像中的极值点。
2. 尺度空间极值点检测：在检测图像中的极值点时，需要编写代码实现对高斯差分金字塔的极值点检测算法，以找到图像中的潜在关键点。
3. 关键点的定位：关键点定位阶段需要编写代码实现对尺度空间极值点的精确定位，消除低对比度点和边缘响应点，并进行关键点的精细化操作。
4. 关键点的方向确定和描述子生成：在这一步骤中，需要编写代码实现对关键点周围区域的梯度幅度信息的计算和关键点方向的确定，以及生成关键点的描述子。
% [descriptors, locs] = sift(img)%% This function returns IMAGE's SIFT keypoints.% Input parameters:% img: the image.%% Returned:% descriptors: a K-by-128 matrix, where each row gives an invariant% descriptor for one of the K keypoints. The descriptor is a vector % of 128 values normalized to unit length.% locs: K-by-4 matrix, in which each row has the 4 values for a% keypoint location (row, column, scale, orientation). The% orientation is in the range [-PI, PI] radians.%% Credits: Thanks for initial version of this program to D. Alvaro and% J.J. Guerrero, Universidad de Zaragoza (modified by D. Lowe) function [descriptors, locs] = sift(img)% If you have the Image Processing Toolbox, you can uncomment the following % lines to allow input of color images, which will be converted to grayscale. if isrgb(img)img = rgb2gray(img);end[rows, cols] = size(img);% Convert into PGM imagefile, readable by "keypoints" executablef = fopen('tmp.pgm', 'w');if f == -1error('Could not create file tmp.pgm.');endfprintf(f, 'P5\n%d\n%d\n255\n', cols, rows);fwrite(f, img', 'uint8');fclose(f);% Call keypoints executableif isunixcommand = '!./sift ';elsecommand = '!siftWin32 ';endcommand = [command ' <tmp.pgm >tmp.key'];eval(command);% Open tmp.key and check its headerg = fopen('tmp.key', 'r');if g == -1error('Could not open file tmp.key.');end[header, count] = fscanf(g, '%d %d', [1 2]);if count ~= 2error('Invalid keypoint file beginning.');endnum = header(1);len = header(2);if len ~= 128error('Keypoint descriptor length invalid (should be 128).');end% Creates the two output matrices (use known size for efficiency) locs = double(zeros(num, 4));descriptors = double(zeros(num, 128));% Parse tmp.keyfor i = 1:num[vector, count] = fscanf(g, '%f %f %f %f', [1 4]); %row col scale ori if count ~= 4error('Invalid keypoint file format');endlocs(i, :) = vector(1, :);[descrip, count] = fscanf(g, '%d', [1 len]);if (count ~= 128)error('Invalid keypoint file value.');end% Normalize each input vector to unit lengthdescrip = descrip / sqrt(sum(descrip.^2));descriptors(i, :) = descrip(1, :);endfclose(g);delete('tmp.pgm');。
sift 特征点检测matlab函数SIFT特征点检测MATLAB函数引言：SIFT（Scale-Invariant Feature Transform）是一种常用的图像特征点检测算法，它具有尺度不变性和旋转不变性等优点，在计算机视觉和图像处理领域得到广泛应用。
本文将介绍SIFT特征点检测的MATLAB函数及其使用方法。
一、什么是SIFT特征点检测？SIFT特征点检测是一种基于局部特征的图像处理方法，它通过在图像中检测出具有稳定尺度和旋转不变性的关键点，从而实现图像的特征匹配、目标识别等应用。
SIFT算法主要包括关键点检测和特征描述两个步骤。
二、SIFT特征点检测的MATLAB函数在MATLAB中，我们可以使用vlfeat工具箱提供的函数来实现SIFT 特征点的检测。
具体来说，vl_sift函数可以用于检测图像中的SIFT特征点，并计算出每个特征点的主要特征描述子。
三、vl_sift函数的使用方法1. 函数调用格式[f, d] = vl_sift(I)其中，I为待检测的图像，返回值f为检测到的特征点的位置和尺度信息，d为特征点的主要描述子。
2. 示例代码下面是一个简单的示例代码，展示了如何使用vl_sift函数进行SIFT特征点检测。
```matlabI = imread('image.jpg'); % 读取图像I = single(rgb2gray(I)); % 转为灰度图像[f, d] = vl_sift(I); % SIFT特征点检测imshow(I); % 显示原图像hold on;vl_plotframe(f); % 显示特征点hold off;```通过上述代码，我们可以将图像中的SIFT特征点检测出来，并在原图像上用红色圆圈标出。
四、SIFT特征点检测的应用SIFT特征点检测在计算机视觉和图像处理领域有着广泛的应用。
以下是一些常见的应用场景：1. 特征匹配：通过比较两幅图像的SIFT特征点，可以进行图像的特征匹配，用于图像拼接、目标跟踪等任务。
sift matlab 代码Sift Matlab 代码Sift算法（Scale-invariant feature transform）是一种用于图像处理和计算机视觉中的特征提取算法。
Matlab是一种常用的编程语言和软件工具，用于实现各种算法和程序。
本文将介绍如何使用Matlab 编写Sift算法的代码，并对其原理和实现进行详细说明。
一、Sift算法原理Sift算法是一种基于局部特征的图像匹配算法，其主要思想是在图像中寻找关键点（keypoints），并对这些关键点进行描述，以便在不同图像之间进行匹配。
Sift算法具有尺度不变性和旋转不变性，能够稳定地提取图像的特征。
Sift算法的主要步骤包括尺度空间极值检测、关键点定位、关键点描述和关键点匹配等。
在尺度空间极值检测中，算法会在图像的不同尺度下检测局部极值点，这些点被认为是关键点的候选。
然后，通过对这些候选点进行精确定位和剔除不稳定点，最终得到真正的关键点。
接着，算法会对每个关键点周围的图像区域进行描述，生成描述子向量。
最后，通过比较不同图像的关键点描述子，实现图像匹配。
二、Matlab实现Sift算法在Matlab中实现Sift算法，可以使用现成的开源库或者自己编写代码。
一种常见的做法是使用vlfeat工具包，该工具包提供了Sift算法的Matlab接口，并包含了各种图像处理和特征提取的函数。
需要安装vlfeat工具包，并将其添加到Matlab的路径中。
然后，可以使用vl_sift函数来提取图像的Sift特征。
该函数会返回关键点的位置、尺度、方向以及描述子等信息，可以根据需要进行进一步处理和分析。
除了使用vlfeat工具包，也可以根据Sift算法的原理，编写自己的Sift代码。
在Matlab中，可以通过图像梯度计算、高斯金字塔构建和关键点描述等步骤，实现Sift算法的各个功能模块。
三、总结本文介绍了Sift算法的原理和在Matlab中的实现方法。
标题：MATLAB中的SIFT角点提取一、介绍SIFT角点提取SIFT（Scale-Invariant Feature Transform）是一种用于计算图像局部特征的算法，它能够在不同尺度和旋转下保持稳定性，并且对于光照变化具有一定的鲁棒性。
SIFT角点提取算法是计算机视觉领域中十分重要的一部分，能够对图像中的关键点进行稳健提取，因此在图像配准、目标检测和三维重建等领域有着广泛的应用。
二、MATLAB中的SIFT角点提取函数在MATLAB中，可以通过vlfeat工具箱来实现SIFT角点提取。
vlfeat 是一个视觉库，其中包括了SIFT、SURF等经典的图像特征提取算法。
具体来说，vlfeat中提供的函数vl_sift可以用于在MATLAB环境下实现SIFT特征提取。
通过调用vl_sift函数，在图像中提取出关键点的坐标、尺度和方向信息，并计算出各个关键点的描述子。
三、MATLAB中如何进行SIFT角点提取1. 准备图像数据在使用MATLAB进行SIFT角点提取之前，首先需要准备好图像数据。
可以使用imread函数读取图像，将其转换为灰度图像或彩色图像均可。
2. 加载vlfeat库通过MATLAB的addpath函数将vlfeat库添加到MATLAB的工作路径中，以便调用其中的函数。
3. 调用vl_sift函数使用vl_sift函数可以对图像进行SIFT角点提取，其调用方式如下：[f, d] = vl_sift(I)其中I为输入的图像，f为提取出的特征点信息，d为各个特征点的描述子。
4. 可视化可以通过vl_plotframe函数将提取出的关键点在图像上进行可视化显示，以便于观察提取结果。
四、SIFT角点提取的应用示例SIFT角点提取在图像配准中有着广泛的应用，下面以图像配准为例进行说明。
1. 加载图像加载两幅待配准的图像，通常我们会选择一幅作为基准图像，另一幅作为待配准图像。
2. SIFT角点提取对两幅图像分别调用vl_sift函数，提取出它们的SIFT特征点。
同名点配准算法 python同名点配准算法是一种在计算机视觉和图像处理中常用的技术，用于将两个或多个图像中的对应点进行匹配。
在Python中，可以使用OpenCV库来实现同名点配准算法。
以下是使用Python和OpenCV实现同名点配准算法的示例代码：pythonimport cv2import numpy as np# 读取两个需要配准的图像img1 = cv2.imread('image1.jpg')img2 = cv2.imread('image2.jpg')# 定义特征检测器feature_detector = cv2.SIFT_create()# 检测并计算关键点和描述符keypoints1, descriptors1 = feature_detector.detectAndCompute(img1, None) keypoints2, descriptors2 = feature_detector.detectAndCompute(img2, None)# 匹配描述符matcher = cv2.BFMatcher()matches = matcher.knnMatch(descriptors1, descriptors2, k=2)# 按照距离排序并去除不良匹配good_matches = []for m, n in matches:if m.distance < 0.75*n.distance:good_matches.append([m])# 提取匹配点的坐标src_pts = np.float32([keypoints1[m[0].queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)dst_pts = np.float32([keypoints2[m[0].trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)# 计算仿射变换矩阵M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)# 应用仿射变换height, width = img1.shape[:2]warped_img = cv2.warpPerspective(img2, M, (width, height))在上述代码中，首先使用SIFT算法检测并计算两个图像的关键点和描述符，然后使用BFMatcher进行匹配。
经典算法SIFT实现即代码解释：以下便是sift源码库编译后的效果图：为了给有兴趣实现sift算法的朋友提供个参考，特整理此文如下。
要了解什么是sift算法，请参考：九、图像特征提取与匹配之SIFT算法。
ok，咱们下面，就来利用Rob Hess维护的sift 库来实现sift算法：首先，请下载Rob Hess维护的sift 库：/hess/code/sift/下载Rob Hess的这个压缩包后，如果直接解压缩，直接编译，那么会出现下面的错误提示：编译提示:error C1083: Cannot open include file: 'cxcore.h': No such file or directory，找不到这个头文件。
这个错误，是因为你还没有安装opencv，因为：cxcore.h和cv.h是开源的OPEN CV头文件,不是VC++的默认安装文件,所以你还得下载OpenCV并进行安装。
然后，可以在OpenCV文件夹下找到你所需要的头文件了。
据网友称，截止2010年4月4日，还没有在VC6.0下成功使用opencv2.0的案例。
所以，如果你是VC6.0的用户请下载opencv1.0版本。
vs的话，opencv2.0,1.0任意下载。
以下，咱们就以vc6.0为平台举例，下载并安装opencv1.0版本、gsl等。
当然，你也可以用vs编译，同样下载opencv（具体版本不受限制）、gsl等。
请按以下步骤操作：一、下载opencv1.0/projects/opencvlibrary/files/opencv-win/1.0/OpenCV_1.0.exe/download二、安装opencv1.0，配置Windows环境变量1、安装注意：假如你是将OpenCV安装到C:/Program Files/OpenCV（如果你安装的时候选择不是安装在C盘，则下面所有对应的C盘都改为你所安装在的那个“X盘”，即可），在安装时选择"将/OpenCV/bin加入系统变量"，打上“勾”。
SIFT特征匹配处理⼀、SIFT算法特征原理SIFT即尺度不变特征转换，它⽤来检测图像的局部性特征，在空间尺度中寻找极值点，提取这点的位置、尺度、旋转不变量。
这些关键点是⼀些⼗分突出，不会因光照和噪⾳等因素⽽变化的点，如⾓点、边缘点、暗区的亮点及亮区的暗点等，所以与影像的⼤⼩和旋转⽆关，对光线、噪声、视⾓改变的容忍度也很⾼。
SIFT特征检测有四步：1.尺度空间的极值检测：搜索所有尺度空间上的图像，通过⾼斯微分函数来识别潜在的对尺度和选择不变的兴趣点。
2.特征点定位：在每个候选的位置上，通过⼀个拟合精细模型来确定位置尺度，关键点的选取依据他们的稳定程度。
3.特征⽅向赋值：基于图像局部的梯度⽅向，分配给每个关键点位置⼀个或多个⽅向，后续的所有操作都是对于关键点的⽅向、尺度和位置进⾏变换，从⽽提供这些特征的不变性。
4.特征点描述：在每个特征点周围的邻域内，在选定的尺度上测量图像的局部梯度，这些梯度被变换成⼀种表⽰，这种表⽰允许⽐较⼤的局部形状的变形和光照变换。
⼆、SIFT特征匹配处理 对两张图⽚进⾏SIFT特征匹配处理 1.两张差异较⼤的集美⼤学尚⼤楼 原图: 图1 图2 SIFT特征匹配： 图3 2. 两张差异较⼩的集美⼤学尚⼤楼 原图： 图4 图5 SIFT特征匹配：代码：1import io2from PIL import Image, ImageTk3import tkinter as tk45import cv26import numpy as np7 MIN_MATCH_COUNT = 489 img1 = cv2.imread("C:/Users/w/PycharmProjects/sift/picture/1.jpg")10 img2 = cv2.imread("C:/Users/w/PycharmProjects/sift/picture/16.jpg")11 g1 = cv2.cvtColor(img1, cv2.COLOR_BGR2GRAY)12 g2 = cv2.cvtColor(img2, cv2.COLOR_BGR2GRAY)13 sift = cv2.xfeatures2d.SIFT_create()14 match = cv2.FlannBasedMatcher(dict(algorithm =2, trees =1), {})15 kp1, de1 = sift.detectAndCompute(g1,None)16 kp2, de2 = sift.detectAndCompute(g2,None)17 m = match.knnMatch(de1, de2, 2)18 m = sorted(m,key = lambda x:x[0].distance)19 ok = [m1 for (m1, m2) in m if m1.distance < 0.7 * m2.distance]20 med = cv2.drawMatches(img1, kp1, img2, kp2, ok, None)2122 cv2.imwrite("C:/Users/w/PycharmProjects/sift/picture/b.jpg", med)23#24# cv2.imshow("0", med)25# cv2.waitKey()26# cv2.destroyAllWindows()272829def resize(w, h, w_box, h_box, pil_image):30 f1 = 1.0 * w_box / w # 1.0 forces float division in Python231 f2 = 1.0 * h_box / h32 factor = min([f1, f2])33 width = int(w * factor)34 height = int(h * factor)35return pil_image.resize((width, height), Image.ANTIALIAS)3637 root = ()38# size of image display box you want39# 期望图像显⽰的⼤⼩40 w_box = 80041 h_box = 10004243# 以⼀个PIL图像对象打开44 pil_image = Image.open(r'C:/Users/w/PycharmProjects/sift/picture/b.jpg')4546# get the size of the image47# 获取图像的原始⼤⼩48 w, h = pil_image.size4950# resize the image so it retains its aspect ration51# but fits into the specified display box52# 缩放图像让它保持⽐例，同时限制在⼀个矩形框范围内53 pil_image_resized = resize(w, h, w_box, h_box, pil_image)5455# convert PIL image object to Tkinter PhotoImage object56# 把PIL图像对象转变为Tkinter的PhotoImage对象57 tk_image = ImageTk.PhotoImage(pil_image_resized)5859# put the image on a widget the size of the specified display box60# Label: 这个⼩⼯具，就是个显⽰框，⼩窗⼝，把图像⼤⼩显⽰到指定的显⽰框61 label = bel(root, image=tk_image, width=w_box, height=h_box)62# padx,pady是图像与窗⼝边缘的距离63 label.pack(padx=5, pady=5)64 root.mainloop()SIFT Code三、SIFT和Harris特征匹配处理的对⽐ 例如，下图是福建省厦门市集美⼤学的尚⼤楼的图⽚：⽤sift算法和Harris算法找到关键点并绘制关键点： SIFT 算法 Harris 算法SIFT 代码：1 import cv22 import numpy as np34 img = cv2.imread('C:/Users/w/PycharmProjects/untitled2/11.jpg')5 gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)6 sift = cv2.xfeatures2d.SIFT_create()78 kp = sift.detect(gray, None) # 找到关键点910 img = cv2.drawKeypoints(gray, kp, img) # 绘制关键点1112 cv2.imshow('sp', img)13 cv2.waitKey(0)SIFT CodeHarris 代码：1 import cv22 import numpy as np34 filename = 'C:/Users/w/PycharmProjects/sift/picture/11.jpg'56 img = cv2.imread(filename)7 gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)8 gray = np.float32(gray)9 #图像转换为float3210 dst = cv2.cornerHarris(gray,2,3,0.04)11 #result is dilated for marking the corners, not important12 dst = cv2.dilate(dst,None)#图像膨胀13 # Threshold for an optimal value, it may vary depending on the image.14 #print(dst)15 #img[dst>0.00000001*dst.max()]=[0,0,255] #可以试试这个参数，⾓点被标记的多余了⼀些16 img[dst>0.01*dst.max()]=[0,0,255]#⾓点位置⽤红⾊标记17 #这⾥的打分值以⼤于0.01×dst 中最⼤值为边界1819 cv2.imshow('dst',img)20 if cv2.waitKey(0) & 0xff == 27:21 cv2.destroyAllWindows()Harris Code SIFT 特征检测：SIFT 从理论上说是⼀种相似不变量，即对图像尺度变化和旋转是不变量。
有关SIFT函数的关键代码(2010-06-08 13:35:18)转载分类：杂感标签：sift代码关键点杂谈function [ pos, scale, orient, desc ] = SIFT( im, octaves, intervals, object_mask, contrast_threshold, curvature_threshold, interactive )% [ pos, scale, orient, desc ] = SIFT( im, octaves, intervals, object_mask, contrast_threshold, curvature_threshold, interactive )if ~exist('octaves')octaves = 4;endif ~exist('intervals')intervals = 2;endif ~exist('object_mask')object_mask = ones(size(im));endif size(object_mask) ~= size(im)object_mask = ones(size(im));endif ~exist('contrast_threshold')contrast_threshold = 0.02;endif ~exist('curvature_threshold')curvature_threshold = 10.0;endif ~exist('interactive')interactive = 1;end% Check that the image is normalized to [0,1]if( (min(im(:)) < 0) | (max(im(:)) > 1) )fprintf( 2, 'Warning: image not normalized to [0,1].\n' );end% Blur the image with a standard deviation of 0.5 to prevent aliasing% and then upsample the image by a factor of 2 using linear interpolation.% Lowe claims that this increases the number of stable keypoints by% a factor of 4.if interactive >= 1fprintf( 2, 'Doubling image size for first octave...\n' );endtic;antialias_sigma = 0.5;if antialias_sigma == 0signal = im;elseg = gaussian_filter( antialias_sigma );%至少比0.5大if exist('corrsep') == 3signal = corrsep( g, g, im );elsesignal = conv2( g, g, im, 'same' );endendsignal = im;[X Y] = meshgrid( 1:0.5:size(signal,2), 1:0.5:size(signal,1) );signal = interp2( signal, X, Y, '*linear' );subsample = [0.5]; % subsampling rate for doubled image is 1/2if interactive >= 1fprintf( 2, 'Prebluring image...\n' );endpreblur_sigma = sqrt(sqrt(2)^2 - (2*antialias_sigma)^2);if preblur_sigma == 0gauss_pyr{1,1} = signal;elseg = gaussian_filter( preblur_sigma );if exist('corrsep') == 3gauss_pyr{1,1} = corrsep( g, g, signal );elsegauss_pyr{1,1} = conv2( g, g, signal, 'same' );endendclear signalpre_time = toc;if interactive >= 1fprintf( 2, 'Preprocessing time %.2f seconds.\n', pre_time );end% The initial blurring for the first image of the first octave of the pyramid.initial_sigma = sqrt( (2*antialias_sigma)^2 + preblur_sigma^2 );% Keep track of the absolute sigma for the octave and scaleabsolute_sigma = zeros(octaves,intervals+3);absolute_sigma(1,1) = initial_sigma * subsample(1);% Keep track of the filter sizes and standard deviations used to generate the pyramid filter_size = zeros(octaves,intervals+3);filter_sigma = zeros(octaves,intervals+3);% Generate the remaining levels of the geometrically sampled gaussian and DOG pyramidsif interactive >= 1fprintf( 2, 'Expanding the Gaussian and DOG pyramids...\n' );endtic;for octave = 1:octavesif interactive >= 1fprintf( 2, '\tProcessing octave %d: image size %d x %d subsample %.1f\n', octave, size(gauss_pyr{octave,1},2), size(gauss_pyr{octave,1},1), subsample(octave) );fprintf( 2, '\t\tInterval 1 sigma %f\n', absolute_sigma(octave,1) );endsigma = initial_sigma;g = gaussian_filter( sigma );filter_size( octave, 1 ) = length(g);filter_sigma( octave, 1 ) = sigma;DOG_pyr{octave} = zeros(size(gauss_pyr{octave,1},1),size(gauss_pyr{octave,1},2),intervals+2);%M*N*Zfor interval = 2:(intervals+3)%从第二层开始方便差分DOG计算% Compute the standard deviation of the gaussian filter needed to% produce the% next level of the geometrically sampled pyramid. Here, sigma_i+1 = k*sigma.% By definition of successive convolution, the required blurring sigma_f to% produce sigma_i+1 from sigma_i is:%% sigma_i+1^2 = sigma_f,i^2 + sigma_i^2% (k*sigma_i)^2 = sigma_f,i^2 + sigma_i^2%% therefore:%% sigma_f,i = sqrt(k^2 - 1)sigma_i%% where k = 2^(1/intervals) to span the octave, so:%% sigma_f,i = sqrt(2^(2/intervals) - 1)sigma_isigma_f = sqrt(2^(2/intervals) - 1)*sigma;g = gaussian_filter( sigma_f );sigma = (2^(1/intervals))*sigma;%得到下一个sigma% Keep track of the absolute sigmaabsolute_sigma(octave,interval) = sigma * subsample(octave);% Store the size and standard deviation of the filter for later usefilter_size(octave,interval) = length(g);filter_sigma(octave,interval) = sigma;if exist('corrsep') == 3gauss_pyr{octave,interval} = corrsep( g, g, gauss_pyr{octave,interval-1} );else%卷积后获得当前层的高斯金子塔gauss_pyr{octave,interval} = conv2( g, g, gauss_pyr{octave,interval-1}, 'same' );end%获得DOG金字塔DOG_pyr{octave}(:,:,interval-1) = gauss_pyr{octave,interval} - gauss_pyr{octave,interval-1};if interactive >= 1fprintf( 2, '\t\tInterval %d sigma %f\n', interval, absolute_sigma(octave,interval) );endendif octave < octaves% The gaussian image 2 images from the top of the stack for% this octave have be blurred by 2*sigma. Subsample this image by a% factor of 2 to procuduce the first image of the next octave.sz = size(gauss_pyr{octave,intervals+1});[X Y] = meshgrid( 1:2:sz(2), 1:2:sz(1) );gauss_pyr{octave+1,1} = interp2(gauss_pyr{octave,intervals+1},X,Y,'*nearest');absolute_sigma(octave+1,1) = absolute_sigma(octave,intervals+1);subsample = [subsample subsample(end)*2];endendpyr_time = toc;if interactive >= 1fprintf( 2, 'Pryamid processing time %.2f seconds.\n', pyr_time );end% Display the gaussian pyramid when in interactive modeif interactive >= 2sz = zeros(1,2);sz(2) = (intervals+3)*size(gauss_pyr{1,1},2);for octave = 1:octavessz(1) = sz(1) + size(gauss_pyr{octave,1},1);endpic = zeros(sz);y = 1;for octave = 1:octaves%显示所有组所有层的图像x = 1;sz = size(gauss_pyr{octave,1});for interval = 1:(intervals + 3)pic(y:(y+sz(1)-1),x:(x+sz(2)-1)) = gauss_pyr{octave,interval};x = x + sz(2);endy = y + sz(1);endfig = figure;clf;imshow(pic);resizeImageFig( fig, size(pic), 0.25 );fprintf( 2, 'The gaussian pyramid (0.25 scale).\nPress any key to continue.\n' ); pause;close(fig)end% Display the DOG pyramid when in interactive modeif interactive >= 2sz = zeros(1,2);sz(2) = (intervals+2)*size(DOG_pyr{1}(:,:,1),2);for octave = 1:octavessz(1) = sz(1) + size(DOG_pyr{octave}(:,:,1),1);endpic = zeros(sz);y = 1;for octave = 1:octavesx = 1;sz = size(DOG_pyr{octave}(:,:,1));for interval = 1:(intervals + 2)pic(y:(y+sz(1)-1),x:(x+sz(2)-1)) = DOG_pyr{octave}(:,:,interval);x = x + sz(2);endy = y + sz(1);endfig = figure;clf;imshow(pic);resizeImageFig( fig, size(pic), 0.25 );fprintf( 2, 'The DOG pyramid (0.25 scale).\nPress any key to continue.\n' ); pause;close(fig)end%已经可能到这边xx = [ 1 -2 1 ];yy = xx';xy = [ 1 0 -1; 0 0 0; -1 0 1 ]/4;%| Dxx Dxy |%| |%| Dxy Dyy |%| |% Coordinates of keypoints after each stage of processing for display% in interactive mode.raw_keypoints = [];contrast_keypoints = [];curve_keypoints = [];% Detect local maxima in the DOG pyramidif interactive >= 1fprintf( 2, 'Locating keypoints...\n' );endtic;loc = cell(size(DOG_pyr)); % boolean maps of keypointsfor octave = 1:octavesif interactive >= 1fprintf( 2, '\tProcessing octave %d\n', octave );endfor interval = 2:(intervals+1)keypoint_count = 0;contrast_mask = abs(DOG_pyr{octave}(:,:,interval)) >= contrast_threshold;loc{octave,interval} = zeros(size(DOG_pyr{octave}(:,:,interval)));if exist('corrsep') == 3edge = 1;elseedge = ceil(filter_size(octave,interval)/2);endfor y=(1+edge):(size(DOG_pyr{octave}(:,:,interval),1)-edge)for x=(1+edge):(size(DOG_pyr{octave}(:,:,interval),2)-edge)% Only check for extrema where the object mask is 1if object_mask(round(y*subsample(octave)),round(x*subsample(octave))) == 1if( (interactive >= 2) | (contrast_mask(y,x) == 1) )tmp = DOG_pyr{octave}((y-1):(y+1),(x-1):(x+1),(interval-1):(interval+1));pt_val = tmp(2,2,2);if( (pt_val == min(tmp(:))) | (pt_val == max(tmp(:))) )% The point is a local extrema of the DOG image. Store its coordinates for% displaying keypoint location in interactive mode.raw_keypoints = [raw_keypoints; x*subsample(octave) y*subsample(octave)];if abs(DOG_pyr{octave}(y,x,interval)) >= contrast_threshold% The DOG image at the extrema is above the contrast threshold. Store% its coordinates for displaying keypoint locations in interactive mode.contrast_keypoints = [contrast_keypoints; raw_keypoints(end,:)];%最后一行加入对比度% Compute the entries of the Hessian matrix at the extrema location.Dxx = sum(DOG_pyr{octave}(y,x-1:x+1,interval) .* xx);Dyy = sum(DOG_pyr{octave}(y-1:y+1,x,interval) .* yy);Dxy = sum(sum(DOG_pyr{octave}(y-1:y+1,x-1:x+1,interval) .* xy));% Compute the trace and the determinant of the Hessian.Tr_H = Dxx + Dyy;Det_H = Dxx*Dyy - Dxy^2;% Compute the ratio of the principal curvatures.curvature_ratio = (Tr_H^2)/Det_H;if ((Det_H >= 0) & (curvature_ratio < curvature_threshold))% The ratio of principal curvatures is below the threshold (i.e.,% it is not an edge point). Store its coordianates for displaying% keypoint locations in interactive mode.curve_keypoints = [curve_keypoints; raw_keypoints(end,:)];%最后一行加入去边缘关键点矩阵% Set the loc map to 1 to at this point to indicate a keypoint.loc{octave,interval}(y,x) = 1;keypoint_count = keypoint_count + 1;endendendendendendendif interactive >= 1fprintf( 2, '\t\t%d keypoints found on interval %d\n', keypoint_count, interval );endendendkeypoint_time = toc;if interactive >= 1fprintf( 2, 'Keypoint location time %.2f seconds.\n', keypoint_time );end% Display results of extrema detection and keypoint filtering in interactive mode.if interactive >= 2fig = figure;clf;imshow(im);hold on;plot(raw_keypoints(:,1),raw_keypoints(:,2),'y+');resizeImageFig( fig, size(im), 2 );fprintf( 2, 'DOG extrema (2x scale).\nPress any key to continue.\n' );pause;close(fig);fig = figure;clf;imshow(im);hold on;plot(contrast_keypoints(:,1),contrast_keypoints(:,2),'y+');resizeImageFig( fig, size(im), 2 );fprintf( 2, 'Keypoints after removing low contrast extrema (2x scale).\nPress any key to continue.\n' );pause;close(fig);fig = figure;clf;imshow(im);hold on;plot(curve_keypoints(:,1),curve_keypoints(:,2),'y+');resizeImageFig( fig, size(im), 2 );fprintf( 2, 'Keypoints after removing edge points using principal curvature filtering (2x scale).\nPress any key to continue.\n' );pause;close(fig);endclear raw_keypoints contrast_keypoints curve_keypointsg = gaussian_filter( 1.5 * absolute_sigma(1,intervals+3) / subsample(1) );zero_pad = ceil( length(g) / 2 );if interactive >= 1fprintf( 2, 'Computing gradient magnitude and orientation...\n' );endtic;mag_thresh = zeros(size(gauss_pyr));mag_pyr = cell(size(gauss_pyr));grad_pyr = cell(size(gauss_pyr));for octave = 1:octavesfor interval = 2:(intervals+1)diff_x = 0.5*(gauss_pyr{octave,interval}(2:(end-1),3:(end))-gauss_pyr{octave,interval}(2:(end-1),1:(end-2) ));diff_y =0.5*(gauss_pyr{octave,interval}(3:(end),2:(end-1))-gauss_pyr{octave,interval}(1:(end-2),2:(end-1) ));mag = zeros(size(gauss_pyr{octave,interval}));mag(2:(end-1),2:(end-1)) = sqrt( diff_x .^ 2 + diff_y .^ 2 );mag_pyr{octave,interval} = zeros(size(mag)+2*zero_pad);mag_pyr{octave,interval}((zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad)) = mag;grad = zeros(size(gauss_pyr{octave,interval}));grad(2:(end-1),2:(end-1)) = atan2( diff_y, diff_x );grad(find(grad == pi)) = -pi;grad_pyr{octave,interval} = zeros(size(grad)+2*zero_pad);grad_pyr{octave,interval}((zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad)) = grad; endendclear mag gradgrad_time = toc;if interactive >= 1fprintf( 2, 'Gradient calculation time %.2f seconds.\n', grad_time );endnum_bins = 36;hist_step = 2*pi/num_bins;hist_orient = [-pi:hist_step:(pi-hist_step)];% Initialize the positions, orientations, and scale information% of the keypoints to emtpy matrices.pos = [];orient = [];scale = [];if interactive >= 1fprintf( 2, 'Assigining keypoint orientations...\n' );endtic;for octave = 1:octavesif interactive >= 1fprintf( 2, '\tProcessing octave %d\n', octave );endfor interval = 2:(intervals + 1)if interactive >= 1fprintf( 2, '\t\tProcessing interval %d ', interval );endkeypoint_count = 0;g = gaussian_filter( 1.5 * absolute_sigma(octave,interval)/subsample(octave) );hf_sz = floor(length(g)/2);g = g'*g;loc_pad = zeros(size(loc{octave,interval})+2*zero_pad);loc_pad((zero_pad+1):(end-zero_pad),(zero_pad+1):(end-zero_pad)) = loc{octave,interval};[iy ix]=find(loc_pad==1);for k = 1:length(iy)x = ix(k);y = iy(k);wght = g.*mag_pyr{octave,interval}((y-hf_sz):(y+hf_sz),(x-hf_sz):(x+hf_sz));grad_window = grad_pyr{octave,interval}((y-hf_sz):(y+hf_sz),(x-hf_sz):(x+hf_sz));orient_hist=zeros(length(hist_orient),1);for bin=1:length(hist_orient)% Compute the diference of the orientations mod pidiff = mod( grad_window - hist_orient(bin) + pi, 2*pi ) - pi;% Accumulate the histogram binsorient_hist(bin)=orient_hist(bin)+sum(sum(wght.*max(1 - abs(diff)/hist_step,0)));%orient_hist(bin)=orient_hist(bin)+sum(sum(wght.*(abs(diff) <= hist_step)));end% Find peaks in the orientation histogram using nonmax suppression.peaks = orient_hist;rot_right = [ peaks(end); peaks(1:end-1) ];rot_left = [ peaks(2:end); peaks(1) ];peaks( find(peaks < rot_right) ) = 0;peaks( find(peaks < rot_left) ) = 0;% Extract the value and index of the largest peak.[max_peak_val ipeak] = max(peaks);peak_val = max_peak_val;while( peak_val > 0.8*max_peak_val )A = [];b = [];for j = -1:1A = [A; (hist_orient(ipeak)+hist_step*j).^2 (hist_orient(ipeak)+hist_step*j) 1];bin = mod( ipeak + j + num_bins - 1, num_bins ) + 1;b = [b; orient_hist(bin)];endc = pinv(A)*b;max_orient = -c(2)/(2*c(1));while( max_orient < -pi )max_orient = max_orient + 2*pi;endwhile( max_orient >= pi )max_orient = max_orient - 2*pi;end% Store the keypoint position, orientation, and scale informationpos = [pos; [(x-zero_pad) (y-zero_pad)]*subsample(octave) ];orient = [orient; max_orient];scale = [scale; octave interval absolute_sigma(octave,interval)];keypoint_count = keypoint_count + 1;% Get the next peakpeaks(ipeak) = 0;[peak_val ipeak] = max(peaks);endendif interactive >= 1fprintf( 2, '(%d keypoints)\n', keypoint_count );endendendclear loc loc_padorient_time = toc;if interactive >= 1fprintf( 2, 'Orientation assignment time %.2f seconds.\n', orient_time );end% Display the keypoints with scale and orientation in interactive mode.if interactive >= 2fig = figure;clf;imshow(im);hold on;display_keypoints( pos, scale(:,3), orient, 'y' );resizeImageFig( fig, size(im), 2 );fprintf( 2, 'Final keypoints with scale and orientation (2x scale).\nPress any key to continue.\n' ); pause;close(fig);endorient_bin_spacing = pi/4;orient_angles = [-pi:orient_bin_spacing:(pi-orient_bin_spacing)];grid_spacing = 4;[x_coords y_coords] = meshgrid( [-6:grid_spacing:6] );feat_grid = [x_coords(:) y_coords(:)]';[x_coords y_coords] = meshgrid( [-(2*grid_spacing-0.5):(2*grid_spacing-0.5)] );feat_samples = [x_coords(:) y_coords(:)]';feat_window = 2*grid_spacing;desc = [];if interactive >= 1fprintf( 2, 'Computing keypoint feature descriptors for %d keypoints', size(pos,1) );endfor k = 1:size(pos,1)x = pos(k,1)/subsample(scale(k,1));y = pos(k,2)/subsample(scale(k,1));% Rotate the grid coordinates.M = [cos(orient(k)) -sin(orient(k)); sin(orient(k)) cos(orient(k))];feat_rot_grid = M*feat_grid + repmat([x; y],1,size(feat_grid,2));feat_rot_samples = M*feat_samples + repmat([x; y],1,size(feat_samples,2));feat_desc = zeros(1,128);for s = 1:size(feat_rot_samples,2)x_sample = feat_rot_samples(1,s);y_sample = feat_rot_samples(2,s);% Interpolate the gradient at the sample position[X Y] = meshgrid( (x_sample-1):(x_sample+1), (y_sample-1):(y_sample+1) );Z = gauss_pyr{scale(k,1),scale(k,2)};[m,n] = size(Z);G = interp2( 1:m,1:n,Z, X, Y, '*linear' );G(find(isnan(G))) = 0;diff_x = 0.5*(G(2,3) - G(2,1));diff_y = 0.5*(G(3,2) - G(1,2));mag_sample = sqrt( diff_x^2 + diff_y^2 );grad_sample = atan2( diff_y, diff_x );if grad_sample == pigrad_sample = -pi;end% Compute the weighting for the x and y dimensions.x_wght = max(1 - (abs(feat_rot_grid(1,:) - x_sample)/grid_spacing), 0);y_wght = max(1 - (abs(feat_rot_grid(2,:) - y_sample)/grid_spacing), 0);pos_wght = reshape(repmat(x_wght.*y_wght,8,1),1,128);diff = mod( grad_sample - orient(k) - orient_angles + pi, 2*pi ) - pi;orient_wght = max(1 - abs(diff)/orient_bin_spacing,0);orient_wght = repmat(orient_wght,1,16);% Compute the gaussian weighting.g = exp(-((x_sample-x)^2+(y_sample-y)^2)/(2*feat_window^2))/(2*pi*feat_window^2);% Accumulate the histogram bins.feat_desc = feat_desc + pos_wght.*orient_wght*g*mag_sample;endfeat_desc = feat_desc / norm(feat_desc);feat_desc( find(feat_desc > 0.2) ) = 0.2;feat_desc = feat_desc / norm(feat_desc);desc = [desc; feat_desc];if (interactive >= 1) & (mod(k,25) == 0)fprintf( 2, '.' );endenddesc_time = toc;sample_offset = -(subsample - 1);for k = 1:size(pos,1)pos(k,:) = pos(k,:) + sample_offset(scale(k,1));endif size(pos,1) > 0scale = scale(:,3);endif interactive >= 1fprintf( 2, '\nDescriptor processing time %.2f seconds.\n', desc_time );fprintf( 2, 'Processing time summary:\n' );fprintf( 2, '\tPreprocessing:\t%.2f s\n', pre_time );fprintf( 2, '\tPyramid:\t%.2f s\n', pyr_time );fprintf( 2, '\tKeypoints:\t%.2f s\n', keypoint_time );fprintf( 2, '\tGradient:\t%.2f s\n', grad_time );fprintf( 2, '\tOrientation:\t%.2f s\n', orient_time );fprintf( 2, '\tDescriptor:\t%.2f s\n', desc_time );fprintf( 2, 'Total processing time %.2f seconds.\n', pre_time + pyr_time + keypoint_time + grad_time + orient_time + desc_time );end。
