Python遥感图像处理应用篇039 GDAL遥感图像角点检测（Harris、Shi-Tomasi、FAST、SIFT、SURF、ORB、MinEngen、Brisk）更新中……

最新推荐文章于 2023-12-16 22:51:18 发布
最新推荐文章于 2023-12-16 22:51:18 发布
阅读量1k 收藏 8
点赞数 5
文章标签： python 图像处理 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/soderayer/article/details/132528830
版权
本文详细讲解了遥感图像角点检测方法，如Harris、Shi-Tomasi、FAST等，并提供了相应的Python实现代码。适合计算机视觉和遥感领域的技术实践。

摘要生成于 C知道 ，由 DeepSeek-R1 满血版支持， 前往体验 >

1.角点检测方法

遥感图像角点检测是一种用于识别图像中角点（也称为兴趣点）的方法，角点通常代表着图像中的显著特征点。以下是一些常用的遥感图像角点检测方法：

  1. Harris 角点检测算法：Harris 算法通过计算图像中每个像素的灰度变化量，并结合局部邻域的自相关矩阵，来判断该点是否为角点。

  2. Shi-Tomasi 角点检测算法：Shi-Tomasi 算法是对 Harris 算法的改进，它选取图像中最显著的角点，通过设置角点响应函数的阈值。

  3. FAST（Features from Accelerated Segment Test）角点检测算法：FAST 角点检测算法是一种基于像素强度的高速角点检测方法，通过利用圆周上像素的灰度变化快速确定角点。

  4. SIFT（Scale-Invariant Feature Transform）角点检测算法：SIFT 是一种局部特征描述符，其中包含了角点检测步骤，它可以在不同尺度和旋转下对角点进行检测。

  5. SURF（Speeded-Up Robust Features）角点检测算法：SURF 是一种快速且具有尺度不变性的图像特征描述符，其中包含了用于角点检测的步骤。

  6. ORB（Oriented FAST and Rotated BRIEF）角点检测算法：ORB 是一种基于 FAST 角点检测器和 BRIEF 描述符的算法，用于检测和描述图像中的角点。

  7. MinEigen （角点检测器）：全称为 "Minimum Eigenvalue" 角点检测器。

    MinEigen 角点检测器是一种基于像素灰度值的方法，它使用了图像中的 Hessian 矩阵来检测角点。Hessian 矩阵描述了图像的局部结构信息，而最小特征值代表了灰度变化最明显的方向。MinEigen 角点检测器的工作原理如下：（1）对于图像中的每个像素点，构建一个局部窗口（通常是一个固定大小的正方形或圆）。（2）在窗口内计算 Hessian 矩阵。（3）计算 Hessian 矩阵的特征值，其中最小特征值表示该像素点的角点响应强度。（4）根据设定的阈值或其他规则来确定是否将该点标识为角点。（5）MinEigen 角点检测器相对于一些其他方法具有较好的响应速度，在计算上较为高效。它在许多计算机视觉应用中被广泛使用，如目标跟踪、图像拼接和三维重建等。

  8. BRISK算法是2011年ICCV上《BRISK:Binary Robust Invariant Scalable Keypoints》文章中，提出来的一种特征提取算法，也是一种二进制的特征描述算子。它具有较好的旋转不变性、尺度不变性，较好的鲁棒性。在图像配准应用中，速度比较：SIFT>SURF>BRISK>FREAK>ORB，在对有较大模糊的图像配准时，BRISK算法在其中表现最为出色。

这些方法在遥感图像处理中常被广泛应用，可以用于目标检测、图像匹配、图像拼接等任务。根据实际需求和图像特征，选择适合的角点检测方法进行处理。不同的角点检测方法有各自的优缺点和适用条件，在实际应用中需要根据具体情况进行选择和调整参数

2.python实现代码及效果展示

2.1 Harris角点检测

代码函数：

get_img_Harris(inputimgfile,k=0.04,window_size=3,IsAllBands=True，band_num=1)

参数：k=0.04,window_size = 3,IsAllBands=False,band_num=1

IsAllBands=False，只计算band_num波段

IsAllBands=True，计算所有波段。

from osgeo import gdal
import numpy as np
from scipy.ndimage import convolve

#定义角点计算函
def get_Harris(image,k,window_size):
    # 运行 Harris 角点检测算法
    # 计算图像的梯度
    dx = np.gradient(image, axis=1)
    dy = np.gradient(image, axis=0)

    # 计算角点响应函数
    Ix2 = dx ** 2
    Iy2 = dy ** 2
    Ixy = dx * dy
    Sx2 = convolve(Ix2, np.ones((window_size, window_size)))
    Sy2 = convolve(Iy2, np.ones((window_size, window_size)))
    Sxy = convolve(Ixy, np.ones((window_size, window_size)))
    det = Sx2 * Sy2 - Sxy ** 2
    trace = Sx2 + Sy2
    response = det - k * trace ** 2

    # 阈值化角点响应函数
    threshold = 0.01 * np.max(response)
    corner_map = np.zeros_like(image)
    corner_map[response > threshold] = 255
    return  corner_map
#k是角点响应参数
def get_img_Harris(inputimgfile,k=0.04,window_size = 3,IsAllBands=False,band_num=1):
    gdal.UseExceptions()
    # 读取影像数据
    dataset = gdal.Open(inputimgfile)
    band_nums=dataset.RasterCount  #波段数
    cols = dataset.RasterXSize
    rows = dataset.RasterYSize

    if IsAllBands==False:  #只计算指定第band_num波段的结果
        band = dataset.GetRasterBand(band_num)
        # 将波段数据加载到一个 NumPy 数组中
        image = band.ReadAsArray()
        corner_map = get_Harris(image, k, window_size)

        # 将结果保存为影像文件
        driver = gdal.GetDriverByName('GTiff')
        outimgfile = inputimgfile.replace('.tif', '_Harris_Band'+str(band_num)+'.img')
        output_dataset = driver.Create(outimgfile, cols, rows, 1, gdal.GDT_Byte)
        output_dataset.GetRasterBand(1).WriteArray(corner_map)

        # 设置投影和地理转换信息
        output_dataset.SetProjection(dataset.GetProjection())
        output_dataset.SetGeoTransform(dataset.GetGeoTransform())

        output_dataset.FlushCache()
        # 关闭数据集
        output_dataset = None
        dataset = None
        print("计算Harris角点已保存为：" + outimgfile)
    else:    #计算所有波段的Harris角点
        # 将结果保存为影像文件
        driver = gdal.GetDriverByName('GTiff')
        outimgfile = inputimgfile.replace('.tif', '_HarrisAllBands.img')
        output_dataset = driver.Create(outimgfile, cols, rows, band_nums, gdal.GDT_Byte)
        for i in range(1,band_nums+1):
            band = dataset.GetRasterBand(i)
            # 将波段数据加载到一个 NumPy 数组中
            image = band.ReadAsArray()
            corner_map = get_Harris(image, k, window_size)
            output_dataset.GetRasterBand(i).WriteArray(corner_map)

        # 设置投影和地理转换信息
        output_dataset.SetProjection(dataset.GetProjection())
        output_dataset.SetGeoTransform(dataset.GetGeoTransform())

        output_dataset.FlushCache()
        # 关闭数据集
        output_dataset = None
        dataset = None
        print("计算Harris角点已保存为：" + outimgfile)

if __name__=='__main__':
    inputimgfile=r'E:\2022桌面\lena_rs.tif'
    get_img_Harris(inputimgfile,k=0.04,window_size=3,IsAllBands=True)

效果：所有波段和第一波段效果图

原图

所有波段计算结果

 第一波段计算结果

2.2 Shi-Tomasi角点检测

和Harris方法类似

代码如下：

from osgeo import gdal
import numpy as np
from scipy.ndimage import convolve

#定义角点计算函
def get_Shi_Tomasi(image,window_size,threshold):
    # 运行 Harris 角点检测算法
    # 计算图像的梯度
    dx = np.gradient(image, axis=1)
    dy = np.gradient(image, axis=0)

    # 计算角点响应函数（Shi-Tomasi）
    Ix2 = dx ** 2
    Iy2 = dy ** 2
    Ixy = dx * dy
    Sx2 = convolve(Ix2, np.ones((3, 3)))
    Sy2 = convolve(Iy2, np.ones((3, 3)))
    Sxy = convolve(Ixy, np.ones((3, 3)))
    det = Sx2 * Sy2 - Sxy ** 2
    trace = Sx2 + Sy2
    # trace[trace < threshold] = threshold  # threshold
    trace[trace == 0] = np.finfo(float).eps  # 防止除以零
    response = np.minimum(det / trace, threshold)

    # 阈值化角点响应函数
    corner_map = np.zeros_like(image)
    corner_map[response > threshold] = 255
    return  corner_map
#k是角点响应参数
def get_img_Shi_Tomasi(inputimgfile,window_size = 3,threshold=0.01,IsAllBands=False,band_num=1):
    gdal.UseExceptions()
    # 读取影像数据
    dataset = gdal.Open(inputimgfile)
    band_nums=dataset.RasterCount  #波段数
    cols = dataset.RasterXSize
    rows = dataset.RasterYSize

    if IsAllBands==False:  #只计算指定第band_num波段的结果
        band = dataset.GetRasterBand(band_num)
        # 将波段数据加载到一个 NumPy 数组中
        image = band.ReadAsArray()
        corner_map = get_Shi_Tomasi(image,threshold, window_size)

        # 将结果保存为影像文件
        driver = gdal.GetDriverByName('GTiff')
        outimgfile = inputimgfile.replace('.tif', '_Shi_Tomasi_Band'+str(band_num)+'.img')
        output_dataset = driver.Create(outimgfile, cols, rows, 1, gdal.GDT_Byte)
        output_dataset.GetRasterBand(1).WriteArray(corner_map)

        # 设置投影和地理转换信息
        output_dataset.SetProjection(dataset.GetProjection())
        output_dataset.SetGeoTransform(dataset.GetGeoTransform())

        output_dataset.FlushCache()
        # 关闭数据集
        output_dataset = None
        dataset = None
        print("计算Harris角点已保存为：" + outimgfile)
    else:    #计算所有波段的Harris角点
        # 将结果保存为影像文件
        driver = gdal.GetDriverByName('GTiff')
        outimgfile = inputimgfile.replace('.tif', '_Shi_TomasiAllBands.img')
        output_dataset = driver.Create(outimgfile, cols, rows, band_nums, gdal.GDT_Byte)
        for i in range(1,band_nums+1):
            band = dataset.GetRasterBand(i)
            # 将波段数据加载到一个 NumPy 数组中
            image = band.ReadAsArray()
            corner_map = get_Shi_Tomasi(image,threshold, window_size)
            output_dataset.GetRasterBand(i).WriteArray(corner_map)

        # 设置投影和地理转换信息
        output_dataset.SetProjection(dataset.GetProjection())
        output_dataset.SetGeoTransform(dataset.GetGeoTransform())

        output_dataset.FlushCache()
        # 关闭数据集
        output_dataset = None
        dataset = None
        print("计算Harris角点已保存为：" + outimgfile)

if __name__=='__main__':
    inputimgfile=r'E:\2022桌面\lena_rs.tif'
    get_img_Shi_Tomasi(inputimgfile,threshold=0.01,window_size=3,IsAllBands=False)

2.3 FAST角点检测

代码

"""
Edward Rosten和Tom Drummond在2006年发表的“Machine learning for high-speed corner detection”文章中提出了一种FAST特征，
并在2010年对这篇论文作了小幅度的修改后重新发表。FAST的全称为Features From Accelerated Segment Test。Rosten等人将FAST角点定义为：
若某像素点与其周围领域内足够多的像素点处于不同的区域，则该像素点可能为角点。也就是某些属性与众不同，考虑灰度图像，
即若该点的灰度值比其周围领域内足够多的像素点的灰度值大或者小，则该点可能为角点。
FAST算法比其他已知的角点检测算法要快很多倍，但是当图片中的噪点较多时，它的健壮性并不好，而且算法的效果还依赖于一个阈值tt。
而且FAST不产生多尺度特征而且FAST特征点没有方向信息，这样就会失去旋转不变性。
"""
from osgeo import gdal
import numpy as np
from skimage.feature import corner_fast

def get_img_FAST(inputimgfile,outimgfile,band_num=1):
    gdal.UseExceptions()
    # 读取影像数据
    dataset = gdal.Open(inputimgfile)
    band = dataset.GetRasterBand(band_num)
    cols = dataset.RasterXSize
    rows = dataset.RasterYSize
    # 将波段数据加载到一个 NumPy 数组中
    data = band.ReadAsArray(0, 0, cols, rows)

    # 运行 FAST 角点检测算法
    corners = corner_fast(data,12,threshold=0.025)

    # 创建一个和原始影像相同尺寸的二值图像，将角点标记为白色
    result = np.zeros((rows, cols), dtype=np.uint8)
    # result[corners] = 255
    corner_indices = np.transpose(np.nonzero(corners))  # 转换为整数索引
    for idx in corner_indices:
        result[idx[0], idx[1]] = 255

    # 将结果保存为影像文件
    driver = gdal.GetDriverByName('GTiff')
    output_dataset = driver.Create(outimgfile, cols, rows, 1, gdal.GDT_Byte)
    output_dataset.GetRasterBand(1).WriteArray(result)

    # 设置投影和地理转换信息
    output_dataset.SetProjection(dataset.GetProjection())
    output_dataset.SetGeoTransform(dataset.GetGeoTransform())

    output_dataset.FlushCache()
    # 关闭数据集
    output_dataset = None
    dataset = None
    print("计算FAST已保存为：" + outimgfile)

if __name__=='__main__':
    inputimgfile=r'E:\2022桌面\lena_grey.jpg'
    #检测图片格式
    outimgfile = inputimgfile.replace('.jpg', '_3FAST.img')
    get_img_FAST(inputimgfile, outimgfile,band_num=1)

2.4 SIFT角点检测

2.5 SURF角点检测

2.6 ORB角点检测

Python遥感图像处理应用篇(三十一)：Python+GDAL实现遥感图像平移、旋转
soderayer的博客
03-20 1650
根据该函数的参数，要实现图像的平移，其实我们只需要修改图像的左上角起始点坐标就能够实现了，比如我们将图像向东易懂5000m，向南移动5000m，那么，我们是需要将图像六元组中的第一个和第四个要素平移一定距离就行了。Y坐标方向是负值，所以我们向南移动就是-5000m，向东移动就是+5000m。元组是不可变的，所以我们需要将元组转换为列表，更改列表，然后将列表转换回元组。运行环境：windows10 pycharm python3.7.7。GDAL 中，GeoTransform是一个六个元素的元组。
Python遥感图像处理应用篇037 GDAL+Scikit-image遥感图像主成分分析PCA
soderayer的博客
08-25 1283
使用多波段遥感图像进行主成分分析，这里使用了6个波段的数据计算和显示效果图如下：左边是原图像IR+R+G显示图，右边是计算得到的3个主成分组合显示的RGB图像。
python读取遥感影像
12-06
遥感影像读取
Python opencv用笔记本摄像头检测道路拐角
xiangxianghehe的博客
12-05 1170
这里利用了Harris角点来检测物体边缘和路口拐角，其原理是：人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果在各个方向上移动这个特征的小窗口，窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生
Python遥感图像处理指南(2)-卫星影像分析MNDWI提取
最新发布
gis收藏家的博客
12-16 1301
在第一篇文章中，我们使用了我们导入的软件包中已定义的函数，如 rasterio 中的 open 和 read 以及 matplotlib 中的 imshow。现在，我们可以编写自己的函数了。函数是一段代码，它接收一些参数作为输入并返回输出（经过一些计算），可以在代码中多次重复使用，也可以被其他代码重用。在 Python 中定义一个新函数非常简单。return v在此语法示例中，arg1和arg2强制参数，应该传递给函数，否则会引发错误，而arg3是一个可选参数，如果是，则假定默认值 1从函数调用中省略。
Python遥感图像处理指南(1)-环境准备及卫星影像简单处理
gis收藏家的博客
12-16 1547
这是我打算发布的关于如何使用 Python 编程语言处理地理空间数据的系列中的第一篇。在该领域工作几年后，我注意到有许多来自地球科学和相关领域的专业人士不是程序员，但在某种程度上需要使用编程语言来执行分析、数学计算或简单地自动化工作流程。通常，这些专业人士很难实现这一目标。如果没有适当的知识，他们最终会开发出复杂的代码，而不遵循容易出错的最佳实践，并且无法有效地处理地理空间领域的大量数据。在另一个极端，我们注意到 IT 专业人员掌握了编程概念，但不具备正确处理地理空间数据的知识。
Python遥感图像处理应用篇036 GDAL+Scikit-image计算遥感图像梯度直方图HOG
soderayer的博客
08-24 902
方向梯度直方图HOG（Histogram of Oriented Gradient, HOG）特征是一种在计算机视觉和图像处理中用来进行物体检测的特征描述子。它通过计算和统计图像局部区域的梯度方向直方图来构成特征。Hog特征结合SVM分类器已经被广泛应用于图像识别中，尤其在行人检测中获得了极大的成功。与其他的特征描述方法相比，HOG有很多优点。首先，由于HOG是在图像的局部方格单元上操作，所以它对图像几何的和光学的形变都能保持很好的不变性，这两种形变只会出现在更大的空间领域上。
Python遥感图像处理应用篇(三十四)：GDAL+Scikit-image+GLCM计算遥感图像纹理特征
soderayer的博客
08-13 1301
在遥感图像处理中，计算纹理特征是一种常见的方法来描述图像中不同区域之间的纹理差异。本文主要是用GLCM来计算遥感图像的纹理特征。
python+gdal+遥感图像拼接(mosaic)的实例
09-17
Python结合GDAL库进行遥感图像拼接，也称为镶嵌，是遥感图像处理中的常见操作，主要用于将多个影像合并成一幅大图。这个过程在处理卫星数据或无人机航拍图像时尤其有用，因为这些数据通常由多个小块图像组成。GDAL...
Python遥感图像处理应用篇(二十二)：Python+GDAL 批量等距离裁剪影像-续
soderayer的博客
11-04 1762
本篇对上一篇博客代码进行少量修改，实现遥感图像按照指定大小等距裁剪。比如我们需要将原始遥感图像按照60km*60km大小标准块裁剪，我们将上一篇中的裁剪函数参数修改一下，将行列号修改为裁剪大小。
opencv - Python 实现角点检测
04-16
基于Opencv的cornerHarris函数实现角点检测
python gdal + skimage实现基于遥感影像的传统图像分割及合并外加矢量化
12-21
根据我前述博客中对图像传分割算法及图像块合并方法的实验探究，在此将这些方法用于遥感影像并尝试矢量化。 这个过程中我自己遇到了一个棘手的问题，在最后的结果那里有描述，希望知道的朋友帮忙解答一下，谢谢！ 直接上代码： # -*- coding: utf-8 -*- import os import cv2 import gdal from osgeo import ogr,osr import numpy as np from skimage import morphology, color from skimage.segmentation import felzenszwalb, slic,
基于GDAL的Python实现遥感影像PCA的代码
09-10
PCA基本步骤： 对数据进行归一化处理（代码中并非这么做的，而是直接减去均值） 计算归一化后的数据集的协方差矩阵 计算协方差矩阵的特征值和特征向量 保留最重要的k个特征（通常k要小于n），也可以自己制定，也可以选择一个阈值，然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值，则选择这个k 找出k个特征值对应的特征向量 将m * n的数据集乘以k个n维的特征向量的特征向量（n * k）,得到最后降维的数据。
遥感图像检测
12-09
遥感图像处理AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
基于Python手动实现Harris角点检测
weixin_46585836的博客
12-01 4227
基于Python手动实现Harris角点检测算法
python 角点检测
buaixvexi的博客
02-25 1852
一. Harris角点 1.Harris角点原理 角点原理来源于人对角点的感性判断，即图像在各个方向灰度有明显变化。算法的核心是利用局部窗口在图像上进行移动判断灰度发生较大的变化。 2.Harris算法数学模型 根据算法思想，构建数学模型，计算移动窗口的的灰度差值 上图中W函数表示窗口函数，M矩阵为偏导数矩阵。对于矩阵可以进行对称矩阵的变化，假设利用两个特征值进行替代，其几何含义类似下图中的表...
【西安】Python-GEE遥感云大数据分析、管理与可视化技术及多领域案例实践应用
WangYan2022的博客
02-22 794
帮助科研工作者掌握Earth Engine的实际应用能力，将以Python编程语言为基础，结合案例从平台搭建、影像数据分析、本地和云端数据管理，以及云端数据论文出版级可视化等方面进行讲解和进阶训练。
遥感图像处理中常用的python操作
热门推荐
QQ704630835的博客
10-24 3万+
本章节主要参考《python地理空间分析指南》第六章。 文章中的所有操作都可以在ENVI中完成，这里只提供一个图像处理的思路。 一、图像波段变换 波段变换最常用的地方就是进行图像显示，例如使用假彩色图像以凸显植被信息。图像波段变换即将图像波段的组合顺序重新排列并显示的方法，主要使用GDAL库。下面以一个例子进行实验，首先给出数据的下载地址： http://git.io/vqs41 打...
Python实现Harris角点检测算法及完整源码
BUG？不存在的！
04-24 959
Harris角点检测算法的核心思想是：对于一幅图像中的局部区域，在不同方向上移动一个小的位移量，如果这个位移量引起的图像灰度变化大，则说明该区域存在角点。具体地，我们可以根据一定的公式计算每个像素点的响应值，然后选取响应值最大的像素点作为角点。该算法可以在图像中检测到角点，这些角点可以用于图像匹配、跟踪、三维重建等领域。本文将介绍如何使用Python实现Harris角点检测算法，并提供完整的源代码。本文介绍了如何使用Python实现Harris角点检测算法，并提供了完整的源代码。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

空中旋转篮球

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值