使用R语言imager包实现Canny边缘检测时非极大值抑制异常求助

我正尝试使用R语言的imager包，通过Canny边缘检测从图像中提取污垢圆圈。我的样本图像中有3个圆圈，已对红色通道进行阈值处理，因为它的区分度最清晰。

注意：imager包的cimg对象为4维结构，若不指定所有维度（如image[i,j,1,1]），R会报错。

我的问题出在*非极大值抑制（non-maximum suppression）*阶段：本想细化检测到的边缘，得到图像中所有边缘的1像素宽骨架（后续还会进行强弱边缘清理），但实际得到的却是梯度幅值图像边缘的椒盐效果。

我用3x3 Sobel核计算得到：

• 边缘梯度幅值图像（代码中变量G）
• 边缘梯度方向图像（代码中变量theta，取值范围-90到90，不易解读）

我的代码如下：

#packages
library(imager)

#define function
non_max <- function(image, edge_direction){
  for (i in 2:(nrow(image)-1)){
    for (j in 2:(ncol(image)-1)){
      #angle 0 (-22.5 to 22.5 degree range)
      if (edge_direction[i,j,1,1] >= -22.5 & 
          edge_direction[i,j,1,1] < 22.5 &
          image[i,j,1,1] < image[i+1,j,1,1] | 
          image[i,j,1,1] < image[i-1,j,1,1]){
        image[i,j,1,1] <- 0
      }
      #angle 45 (22.5 to 67.5)
      else if (edge_direction[i,j,1,1] >= 22.5 & 
            edge_direction[i,j,1,1] < 67.5 &
            image[i,j,1,1] < image[i-1,j-1,1,1] | 
            image[i,j,1,1] < image[i+1,j+1,1,1]){
          image[i,j,1,1] <- 0
      }
      #angle 90 (67.5 to 90)
      else if (edge_direction[i,j,1,1] >= 67.5 & 
               edge_direction[i,j,1,1] < 90 &
               image[i,j,1,1] < image[i,j-1,1,1] | 
               image[i,j,1,1] < image[i,j+1,1,1]){
        image[i,j,1,1] <- 0
      }
      #angle -45 (-22.5 to -67.5)
      else if (edge_direction[i,j,1,1] >= -67.5 & 
               edge_direction[i,j,1,1] < -22.5 &
               image[i,j,1,1] < image[i+1,j-1,1,1] | 
               image[i,j,1,1] < image[i-1,j+1,1,1]){
        image[i,j,1,1] <- 0
      }
      #angle 90 (67.5 to 90)
      else if (edge_direction[i,j,1,1] >= -90 & 
               edge_direction[i,j,1,1] < -67.5 &
               image[i,j,1,1] < image[i,j-1,1,1] | 
               image[i,j,1,1] < image[i,j+1,1,1]){
        image[i,j,1,1] <- 0
      }
    }
  }
  plot(image)
}

#execute function
non_max(G,theta)

我发现如果在每个角度条件中添加额外的IF语句，椒盐效果会减少，但会生成相互平行的额外边缘。示例代码如下，我不确定为何会出现这种情况（除了我对&和|的复杂语句语法理解不足）：

#angle 0 (-22.5 to 22.5 degree range)
if (edge_direction[i,j,1,1] >= -22.5 & 
    edge_direction[i,j,1,1] < 22.5){
  if (image[i,j,1,1] < image[i+1,j,1,1] | 
      image[i,j,1,1] < image[i-1,j,1,1]){
    image[i,j,1,1] <- 0
  }
}

我怀疑边缘方向文件有问题，不过我认为计算是正确的：对y梯度一阶导数（Iy）除以x梯度一阶导数（Ix）取反正切，再乘以180/pi转换为角度：

theta <- atan(Iy/Ix) * (180/pi)

有没有人知道我哪里出错了？我查看了一些Python实现，它们的边缘处理和我的相反，但我确信自己的处理是正确的。

这是我第一次提问，若需其他信息请告知，谢谢！

内容的提问来源于stack exchange，提问作者Michael Sramek