图像处理之角点检测算法(Harris Corner Detection)

简介: 图像处理之角点检测算法(Harris Corner Detection)

图像处理之角点检测算法(Harris Corner Detection)

Harris角点检测是通过数学计算在图像上发现角点特征的一种算法,而且其具有旋转不

性的特质。OpenCV中的Shi-Tomasi角点检测就是基于Harris角点检测改进算法。

基本原理:

角点是一幅图像上最明显与重要的特征,对于一阶导数而言,角点在各个方向的变化是

最大的,而边缘区域在只是某一方向有明显变化。一个直观的图示如下:


数学原理:

基本数学公式如下:

其中W(x, y)表示移动窗口,I(x, y)表示像素灰度值强度,范围为0~255。根据泰勒级数

计算一阶到N阶的偏导数,最终得到一个Harris矩阵公式:


根据Harris的矩阵计算矩阵特征值 ,然后计算Harris角点响应值:

其中K为系数值,通常取值范围为0.04 ~ 0.06之间。

算法详细步骤

第一步:计算图像X方向与Y方向的一阶高斯偏导数Ix与Iy

第二步:根据第一步结果得到Ix^2 , Iy^2与Ix*Iy值

第三步:高斯模糊第二步三个值得到Sxx, Syy, Sxy

第四部:定义每个像素的Harris矩阵,计算出矩阵的两个特质值

第五步:计算出每个像素的R值

第六步:使用3X3或者5X5的窗口,实现非最大值压制

第七步:根据角点检测结果计算,最提取到的关键点以绿色标记,显示在原图上。

程序关键代码解读

第一步计算一阶高斯偏导数的Ix与Iy值代码如下:

  filter.setDirectionType(GaussianDerivativeFilter.X_DIRECTION);
    BufferedImage xImage = filter.filter(grayImage, null);
    getRGB( xImage, 0, 0, width, height, inPixels );
    extractPixelData(inPixels, GaussianDerivativeFilter.X_DIRECTION, height, width);
    
    filter.setDirectionType(GaussianDerivativeFilter.Y_DIRECTION);
    BufferedImage yImage = filter.filter(grayImage, null);
    getRGB( yImage, 0, 0, width, height, inPixels );
    extractPixelData(inPixels, GaussianDerivativeFilter.Y_DIRECTION, height, width);

关于如何计算高斯一阶与二阶偏导数请看这里:

http://blog.csdn.net/jia20003/article/details/16369143

http://blog.csdn.net/jia20003/article/details/7664777

第三步:分别对第二步计算出来的三个值,单独进行高斯

模糊计算,代码如下:

  private void calculateGaussianBlur(int width, int height) {
        int index = 0;
        int radius = (int)window_radius;
        double[][] gw = get2DKernalData(radius, sigma);
        double sumxx = 0, sumyy = 0, sumxy = 0;
        for(int row=0; row<height; row++) {
          for(int col=0; col<width; col++) {            
            for(int subrow =-radius; subrow<=radius; subrow++)
            {
              for(int subcol=-radius; subcol<=radius; subcol++)
              {
                int nrow = row + subrow;
                int ncol = col + subcol;
                if(nrow >= height || nrow < 0)
                {
                  nrow = 0;
                }
                if(ncol >= width || ncol < 0)
                {
                  ncol = 0;
                }
                int index2 = nrow * width + ncol;
                HarrisMatrix whm = harrisMatrixList.get(index2);
                sumxx += (gw[subrow + radius][subcol + radius] * whm.getXGradient());
                sumyy += (gw[subrow + radius][subcol + radius] * whm.getYGradient());
                sumxy += (gw[subrow + radius][subcol + radius] * whm.getIxIy());
              }
            }
            index = row * width + col;
            HarrisMatrix hm = harrisMatrixList.get(index);
            hm.setXGradient(sumxx);
            hm.setYGradient(sumyy);
            hm.setIxIy(sumxy);
            
            // clean up for next loop
            sumxx = 0;
            sumyy = 0;
            sumxy = 0;
          }
        }   
  }

第六步:非最大信号压制(non-max value suppression)

这个在边源检测中是为了得到一个像素宽的边缘,在这里则

是为了得到准确的一个角点像素,去掉非角点值。代码如下:

  /***
   * we still use the 3*3 windows to complete the non-max response value suppression
   */
  private void nonMaxValueSuppression(int width, int height) {
        int index = 0;
        int radius = (int)window_radius;
        for(int row=0; row<height; row++) {
          for(int col=0; col<width; col++) {
            index = row * width + col;
            HarrisMatrix hm = harrisMatrixList.get(index);
            double maxR = hm.getR();
            boolean isMaxR = true;
            for(int subrow =-radius; subrow<=radius; subrow++)
            {
              for(int subcol=-radius; subcol<=radius; subcol++)
              {
                int nrow = row + subrow;
                int ncol = col + subcol;
                if(nrow >= height || nrow < 0)
                {
                  nrow = 0;
                }
                if(ncol >= width || ncol < 0)
                {
                  ncol = 0;
                }
                int index2 = nrow * width + ncol;
                HarrisMatrix hmr = harrisMatrixList.get(index2);
                if(hmr.getR() > maxR)
                {
                  isMaxR = false;
                }
              }             
            }
            if(isMaxR)
            {
              hm.setMax(maxR);
            }
          }
        }
    
  }

运行效果:


程序完整源代码:

package com.gloomyfish.image.harris.corner;
 
import java.awt.image.BufferedImage;
import java.util.ArrayList;
import java.util.List;
 
import com.gloomyfish.filter.study.GrayFilter;
 
public class HarrisCornerDetector extends GrayFilter {
  private GaussianDerivativeFilter filter;
  private List<HarrisMatrix> harrisMatrixList;
  private double lambda = 0.04; // scope : 0.04 ~ 0.06
  
  // i hard code the window size just keep it' size is same as 
  // first order derivation Gaussian window size
  private double sigma = 1; // always
  private double window_radius = 1; // always
  public HarrisCornerDetector() {
    filter = new GaussianDerivativeFilter();
    harrisMatrixList = new ArrayList<HarrisMatrix>();
  }
 
  @Override
  public BufferedImage filter(BufferedImage src, BufferedImage dest) {
    int width = src.getWidth();
        int height = src.getHeight();
        initSettings(height, width);
        if ( dest == null )
            dest = createCompatibleDestImage( src, null );
        
        BufferedImage grayImage = super.filter(src, null);
        int[] inPixels = new int[width*height];
        
    // first step  - Gaussian first-order Derivatives (3 × 3) - X - gradient, (3 × 3) - Y - gradient
    filter.setDirectionType(GaussianDerivativeFilter.X_DIRECTION);
    BufferedImage xImage = filter.filter(grayImage, null);
    getRGB( xImage, 0, 0, width, height, inPixels );
    extractPixelData(inPixels, GaussianDerivativeFilter.X_DIRECTION, height, width);
    
    filter.setDirectionType(GaussianDerivativeFilter.Y_DIRECTION);
    BufferedImage yImage = filter.filter(grayImage, null);
    getRGB( yImage, 0, 0, width, height, inPixels );
    extractPixelData(inPixels, GaussianDerivativeFilter.Y_DIRECTION, height, width);
        
    // second step - calculate the Ix^2, Iy^2 and Ix^Iy
    for(HarrisMatrix hm : harrisMatrixList)
    {
      double Ix = hm.getXGradient();
      double Iy = hm.getYGradient();
      hm.setIxIy(Ix * Iy);
      hm.setXGradient(Ix*Ix);
      hm.setYGradient(Iy*Iy);
    }
    
    // 基于高斯方法,中心点化窗口计算一阶导数和,关键一步 SumIx2, SumIy2 and SumIxIy, 高斯模糊
    calculateGaussianBlur(width, height);
 
    // 求取Harris Matrix 特征值 
    // 计算角度相应值R R= Det(H) - lambda * (Trace(H))^2
    harrisResponse(width, height);
    
    // based on R, compute non-max suppression
    nonMaxValueSuppression(width, height);
    
    // match result to original image and highlight the key points
    int[] outPixels = matchToImage(width, height, src);
    
    // return result image
    setRGB( dest, 0, 0, width, height, outPixels );
    return dest;
  }
  
  
  private int[] matchToImage(int width, int height, BufferedImage src) {
    int[] inPixels = new int[width*height];
        int[] outPixels = new int[width*height];
        getRGB( src, 0, 0, width, height, inPixels );
        int index = 0;
        for(int row=0; row<height; row++) {
          int ta = 0, tr = 0, tg = 0, tb = 0;
          for(int col=0; col<width; col++) {
            index = row * width + col;
            ta = (inPixels[index] >> 24) & 0xff;
                tr = (inPixels[index] >> 16) & 0xff;
                tg = (inPixels[index] >> 8) & 0xff;
                tb = inPixels[index] & 0xff;
                HarrisMatrix hm = harrisMatrixList.get(index);
                if(hm.getMax() > 0)
                {
                  tr = 0;
                  tg = 255; // make it as green for corner key pointers
                  tb = 0;
                  outPixels[index] = (ta << 24) | (tr << 16) | (tg << 8) | tb;
                }
                else
                {
                  outPixels[index] = (ta << 24) | (tr << 16) | (tg << 8) | tb;                  
                }
                
          }
        }
    return outPixels;
  }
  /***
   * we still use the 3*3 windows to complete the non-max response value suppression
   */
  private void nonMaxValueSuppression(int width, int height) {
        int index = 0;
        int radius = (int)window_radius;
        for(int row=0; row<height; row++) {
          for(int col=0; col<width; col++) {
            index = row * width + col;
            HarrisMatrix hm = harrisMatrixList.get(index);
            double maxR = hm.getR();
            boolean isMaxR = true;
            for(int subrow =-radius; subrow<=radius; subrow++)
            {
              for(int subcol=-radius; subcol<=radius; subcol++)
              {
                int nrow = row + subrow;
                int ncol = col + subcol;
                if(nrow >= height || nrow < 0)
                {
                  nrow = 0;
                }
                if(ncol >= width || ncol < 0)
                {
                  ncol = 0;
                }
                int index2 = nrow * width + ncol;
                HarrisMatrix hmr = harrisMatrixList.get(index2);
                if(hmr.getR() > maxR)
                {
                  isMaxR = false;
                }
              }             
            }
            if(isMaxR)
            {
              hm.setMax(maxR);
            }
          }
        }
    
  }
  
  /***
   * 计算两个特征值,然后得到R,公式如下,可以自己推导,关于怎么计算矩阵特征值,请看这里:
   * http://www.sosmath.com/matrix/eigen1/eigen1.html
   * 
   *  A = Sxx;
   *  B = Syy;
   *  C = Sxy*Sxy*4;
   *  lambda = 0.04;
   *  H = (A*B - C) - lambda*(A+B)^2;
     *
   * @param width
   * @param height
   */
  private void harrisResponse(int width, int height) {
        int index = 0;
        for(int row=0; row<height; row++) {
          for(int col=0; col<width; col++) {
            index = row * width + col;
            HarrisMatrix hm = harrisMatrixList.get(index);
            double c =  hm.getIxIy() * hm.getIxIy();
            double ab = hm.getXGradient() * hm.getYGradient();
            double aplusb = hm.getXGradient() + hm.getYGradient();
            double response = (ab -c) - lambda * Math.pow(aplusb, 2);
            hm.setR(response);
          }
        }   
  }
 
  private void calculateGaussianBlur(int width, int height) {
        int index = 0;
        int radius = (int)window_radius;
        double[][] gw = get2DKernalData(radius, sigma);
        double sumxx = 0, sumyy = 0, sumxy = 0;
        for(int row=0; row<height; row++) {
          for(int col=0; col<width; col++) {            
            for(int subrow =-radius; subrow<=radius; subrow++)
            {
              for(int subcol=-radius; subcol<=radius; subcol++)
              {
                int nrow = row + subrow;
                int ncol = col + subcol;
                if(nrow >= height || nrow < 0)
                {
                  nrow = 0;
                }
                if(ncol >= width || ncol < 0)
                {
                  ncol = 0;
                }
                int index2 = nrow * width + ncol;
                HarrisMatrix whm = harrisMatrixList.get(index2);
                sumxx += (gw[subrow + radius][subcol + radius] * whm.getXGradient());
                sumyy += (gw[subrow + radius][subcol + radius] * whm.getYGradient());
                sumxy += (gw[subrow + radius][subcol + radius] * whm.getIxIy());
              }
            }
            index = row * width + col;
            HarrisMatrix hm = harrisMatrixList.get(index);
            hm.setXGradient(sumxx);
            hm.setYGradient(sumyy);
            hm.setIxIy(sumxy);
            
            // clean up for next loop
            sumxx = 0;
            sumyy = 0;
            sumxy = 0;
          }
        }   
  }
  
  public double[][] get2DKernalData(int n, double sigma) {
    int size = 2*n +1;
    double sigma22 = 2*sigma*sigma;
    double sigma22PI = Math.PI * sigma22;
    double[][] kernalData = new double[size][size];
    int row = 0;
    for(int i=-n; i<=n; i++) {
      int column = 0;
      for(int j=-n; j<=n; j++) {
        double xDistance = i*i;
        double yDistance = j*j;
        kernalData[row][column] = Math.exp(-(xDistance + yDistance)/sigma22)/sigma22PI;
        column++;
      }
      row++;
    }
    
//    for(int i=0; i<size; i++) {
//      for(int j=0; j<size; j++) {
//        System.out.print("\t" + kernalData[i][j]);
//      }
//      System.out.println();
//      System.out.println("\t ---------------------------");
//    }
    return kernalData;
  }
 
  private void extractPixelData(int[] pixels, int type, int height, int width)
  {
        int index = 0;
        for(int row=0; row<height; row++) {
          int ta = 0, tr = 0, tg = 0, tb = 0;
          for(int col=0; col<width; col++) {
            index = row * width + col;
            ta = (pixels[index] >> 24) & 0xff;
                tr = (pixels[index] >> 16) & 0xff;
                tg = (pixels[index] >> 8) & 0xff;
                tb = pixels[index] & 0xff;
                HarrisMatrix matrix = harrisMatrixList.get(index);
                if(type == GaussianDerivativeFilter.X_DIRECTION)
                {
                  matrix.setXGradient(tr);
                }
                if(type == GaussianDerivativeFilter.Y_DIRECTION)
                {
                  matrix.setYGradient(tr);
                }
          }
        }
  }
  
  private void initSettings(int height, int width)
  {
        int index = 0;
        for(int row=0; row<height; row++) {
          for(int col=0; col<width; col++) {
            index = row * width + col;
            HarrisMatrix matrix = new HarrisMatrix();
                harrisMatrixList.add(index, matrix);
          }
        }
  }
 
}

最后注意:

我是把彩色图像变为灰度图像来计算,这个计算量小点

处理容易点,此外很多图像处理软件都会用标记来显示

关键点像素,我没有,只是将关键点像素改为绿色。

所以可以从这些方面有很大的提高空间。

相关文章
|
3天前
|
监控 算法 自动驾驶
目标检测算法:从理论到实践的深度探索
【7月更文第18天】目标检测,作为计算机视觉领域的核心任务之一,旨在识别图像或视频中特定对象的位置及其类别。这一技术在自动驾驶、视频监控、医疗影像分析等多个领域发挥着至关重要的作用。本文将深入浅出地介绍目标检测的基本概念、主流算法,并通过一个实际的代码示例,带您领略YOLOv5这一高效目标检测模型的魅力。
37 11
|
8天前
|
机器学习/深度学习 人工智能 算法
计算机视觉:目标检测算法综述
【7月更文挑战第13天】目标检测作为计算机视觉领域的重要研究方向,近年来在深度学习技术的推动下取得了显著进展。然而,面对复杂多变的实际应用场景,仍需不断研究和探索更加高效、鲁棒的目标检测算法。随着技术的不断发展和应用场景的不断拓展,相信目标检测算法将在更多领域发挥重要作用。
|
11天前
|
机器学习/深度学习 数据采集 算法
Python实现Prophet时间序列数据建模与异常值检测(Prophet算法)项目实战
Python实现Prophet时间序列数据建模与异常值检测(Prophet算法)项目实战
|
11天前
|
机器学习/深度学习 算法 计算机视觉
|
16天前
|
机器学习/深度学习 人工智能 文字识别
一种基于YOLOv8改进的高精度红外小目标检测算法 (原创自研)
【7月更文挑战第2天】 💡💡💡创新点: 1)SPD-Conv特别是在处理低分辨率图像和小物体等更困难的任务时优势明显; 2)引入Wasserstein Distance Loss提升小目标检测能力; 3)YOLOv8中的Conv用cvpr2024中的DynamicConv代替;
73 3
|
13天前
|
机器学习/深度学习 运维 算法
Python基于局部离群因子LOF算法(LocalOutlierFactor)实现信用卡数据异常值检测项目实战
Python基于局部离群因子LOF算法(LocalOutlierFactor)实现信用卡数据异常值检测项目实战
|
13天前
|
机器学习/深度学习 数据采集 运维
Python基于孤立森林算法(IsolationForest)实现数据异常值检测项目实战
Python基于孤立森林算法(IsolationForest)实现数据异常值检测项目实战
|
4天前
|
传感器 算法
基于无线传感器网络的MCKP-MMF算法matlab仿真
MCKP-MMF算法是一种启发式流量估计方法,用于寻找无线传感器网络的局部最优解。它从最小配置开始,逐步优化部分解,调整访问点的状态。算法处理访问点的动态影响半径,根据带宽需求调整,以避免拥塞。在MATLAB 2022a中进行了仿真,显示了访问点半径请求变化和代价函数随时间的演变。算法分两阶段:慢启动阶段识别瓶颈并重设半径,随后进入周期性调整阶段,追求最大最小公平性。
基于无线传感器网络的MCKP-MMF算法matlab仿真
|
6天前
|
机器学习/深度学习 算法 数据挖掘
基于改进K-means的网络数据聚类算法matlab仿真
**摘要:** K-means聚类算法分析,利用MATLAB2022a进行实现。算法基于最小化误差平方和,优点在于简单快速,适合大数据集,但易受初始值影响。文中探讨了该依赖性并通过实验展示了随机初始值对结果的敏感性。针对传统算法的局限,提出改进版解决孤点影响和K值选择问题。代码中遍历不同K值,计算距离代价,寻找最优聚类数。最终应用改进后的K-means进行聚类分析。
|
8天前
|
算法 数据安全/隐私保护
基于GA遗传优化算法的Okumura-Hata信道参数估计算法matlab仿真
在MATLAB 2022a中应用遗传算法进行无线通信优化,无水印仿真展示了算法性能。遗传算法源于Holland的理论,用于全局优化,常见于参数估计,如Okumura-Hata模型的传播损耗参数。该模型适用于150 MHz至1500 MHz的频段。算法流程包括选择、交叉、变异等步骤。MATLAB代码执行迭代,计算目标值,更新种群,并计算均方根误差(RMSE)以评估拟合质量。最终结果比较了优化前后的RMSE并显示了SNR估计值。
21 7