图像处理-天体中心确定及边缘拟合

标签：gray 高斯 img 天体边缘圆心图像处理拟合半径

概述

确定一张天体照片的中心点和边缘拟合对于天文学、空间科学以及相关领域的研究和应用具有重要意义。尤其在进行多幅天体图像的对比或组合时，精确的中心点和边缘能够帮助进行图像配准，确保图像正确对齐。

本文基于一张太阳图片，介绍了定位太阳圆心及拟合边缘的过程及实现，该过程来自于论文：

Monitoring Image Alignments and Flat Fields for AIA/SDO Data Images（Richard A. Shine, C. Jacob Wolfson. et al. 2011）

处理过程中涉及的主要技术手段包括：边缘检测、霍夫变换。这里我们将太阳视为一个规则的圆形，因此中心点确定及边缘拟合问题转换成了确定一个圆的圆心及半径问题。

处理流程

依照论文中的描述，其处理流程如下，包括高斯平滑（消噪）、边缘检测、径向梯度边缘处理、霍夫变换。

高斯平滑

高斯平滑又称高斯模糊，是常见的图像去噪的一种方法。它是利用二维高斯函数作为卷积核与原始图像进行卷积运算：

$G(x,y)=\frac{1}{2\pi\sigma ^{2}}e^{-\frac{x^{2}+y^{2}}{2\sigma ^{2}}}$

通俗的讲，就是可以生成一个n*n的卷积核，这个2维卷积核里面的数据是服从高斯分布的，中心点的值最大，越靠近边缘越小，而后卷积核覆盖到图像的每一个像素上，并将核内的值与对应像素的值相乘，求和后得到中心点新的像素值。如下是二维高斯函数以及生成的3*3的高斯核。

公式中， $\sigma$ 是标准差，控制着曲线的宽度。在图像处理中， $\sigma$ 值越大，平滑效果越明显，噪声抑制越强，但同时也会损失更多的图像细节。

（后文代码中引入了一个参数FWHM来控制标准差 $\sigma$ ，其全称为“Full Width at Half Maximum”，即半高全宽，是描述波峰或者波形宽度的一个参数。在高斯函数中，它指的是高斯曲线达到最大高度一半时的宽度。FWHM 与高斯函数的标准差 $\sigma$ 存在确定的数学关系: $FWHM=2\cdot \sqrt{2ln2}\cdot \sigma$ ）

边缘检测：

本次处理在完成高斯平滑后，使用了Sobel算子进行边缘检测，有关边缘检测算子的叙述，可参考：图像处理-边缘检测算法的原理和实现_图像边缘检测卷积核-CSDN博客

图1

由于圆形轮廓内部存在着很多明显的暗斑，因此内部也存在着很多被认为是边缘的信号，于是我们

采用一个圆环掩模来覆盖不适宜的点，只留下大致边缘处的像素点。当然圆环掩模的产生依赖于我们设定的经验圆心和半径。

此外在边缘附近仍然会产生很多不适当的信号，视觉上看，我们的边缘乱且厚。于是使用径向梯度作为掩模进一步筛选边缘点。

径向梯度处理

Sobel滤波的原理是找到图像中颜色（或灰度值）快速变化的点，这些点会被视为边缘点，它使用梯度来量化这种变化率。通过计算每个像素点的梯度值，将那些大于所设阈值的像素点标记为边缘点。

梯度作为一个矢量，不只具有幅值，还具有方向，Sobel滤波的处理只考虑了像素点的梯度值，并未考虑该点的梯度方向，但在一个标准的圆形中，边缘点应是那些沿半径方向灰度值变化最快的点（即梯度最大的点），尤其在天体中。因此我们计算了图像中每个像素点到圆心的角度，并利用Sobel卷积核计算了每个像素点的梯度方向。设定一个较小的度数阈值，认为那些角度与梯度方向之差小于阈值的点为潜在边缘点，由此生成一个掩模，并应用于经Sobel边缘检测后的图像，进一步删选边缘点。得到图2中的第2幅图，第1幅图是我们生成的径向掩模。

图2

迭代霍夫变换

霍夫变化可以检测并识别一个标准圆形，关于霍夫变换的叙述可以参考：图像处理-霍夫变换的原理及实现-CSDN博客

因此我们使用霍夫变换来拟合边缘并得到一个圆心和半径。为了使圆心和半径更加准确，我们采用了多次迭代的方式。具体的：

根据经验提供一个初始圆心和半径，再探查该初始圆边缘一定范围内（一个圆环）的边缘点，利用这些边缘点拟合出一个新的圆，这会提供一个新的圆心和半径，然后利用新圆心和半径继续探查一定范围内（新圆环）的边缘点拟合新圆。每次迭代我们可以减少圆环的厚度，最终，圆心和半径都会趋向收敛，设定一个收敛阈值，认为迭代结束，这通常可以收敛到0.1像素。

Matlab实现

FWHM = 5;
angleThreshold = 20;
delta = FWHM/(2 * sqrt(2 * log(2)));
img_orignal = imread('sun.png'); 
img_gray = rgb2gray(img_orignal);
% 定义经验圆心,半径
centerX = size(img_gray, 1) / 2;
centerY = size(img_gray, 2) / 2;
radius = 350;
%% 高斯平滑
gaussianFilter = fspecial('gaussian', [5, 5], delta); 
img_smooth = imfilter(img_gray, gaussianFilter, 'same');

%% 使用Sobel算子进行边缘检测
sobel_imgs = edge(img_smooth, 'sobel');

%% 生成经验圆环，用作淹模
r1 = 380;  % 大半径
r2 = 330;  % 小半径
[rows, cols] = meshgrid(1:size(img_gray, 2), 1:size(img_gray, 1)); 
distances = sqrt((rows - centerX).^2 + (cols - centerY).^2);
% 生成圆环
ring = (distances <= r1) & (distances >= r2);
sobel_imgsm = sobel_imgs & ring;

%% 径向梯度处理边缘
% 创建网格来计算每个像素的径向梯度
[columns, rows] = meshgrid(1:size(img_gray, 2), 1:size(img_gray, 1));
dx = columns - centerX;
dy = rows - centerY;
% 计算每个像素点到圆心的角度
angles = atan2(dy, dx);
% 计算图像的梯度方向
[gradientX, gradientY] = imgradientxy(img_smooth, 'sobel');
gradientDirection = atan2(gradientY, gradientX);
% 计算梯度方向和点的角度之间的差异
angleDifference = abs(angles - gradientDirection);
angleDifference = min(angleDifference, 2 * pi - angleDifference);
angleThreshold = angleThreshold * (pi / 180); 
% 应用掩模
radial_mask = angleDifference < angleThreshold;
filter_edges = sobel_imgsm & radial_mask;
filter_edges = double(filter_edges);

%% 迭代霍夫
% 参数初始化
initialCenter = [centerX, centerY];
initialRadius = radius; 
deltaR = 10;
deltaC = 10;
radiusRange = 1;

% 初始化最佳估计
bestCenter = initialCenter;
bestRadius = initialRadius;

lastCenter = [0, 0];
lastRadius = 0;

maxIterations = 20;
convergenceThreshold = 0.1;

for iter = 1:maxIterations
    rMin = bestRadius - deltaR;
    rMax = bestRadius + deltaR;
    radii = rMin:radiusRange:rMax;
    
    xMin = bestCenter(1) - deltaC;
    xMax = bestCenter(1) + deltaC;
    yMin = bestCenter(2) - deltaC;
    yMax = bestCenter(2) + deltaC;
    
    [newCenter, newRadius, ~] = houghTransform(filter_edges, radii, xMin:xMax, yMin:yMax);
    
    % 检测收敛
    if max(abs(newCenter - lastCenter)) < convergenceThreshold && abs(newRadius - lastRadius) < convergenceThreshold
        fprintf('Convergence reached after %d iterations.\n', iter);
        break;
    end 
    lastCenter = newCenter;
    lastRadius = newRadius;
    % 更新最佳估计
    bestCenter = newCenter;
    bestRadius = newRadius;   
    % 动态调整搜索范围
    deltaR = deltaR * 0.9;
    deltaC = deltaC * 0.9;
end
center = bestCenter;
radius = bestRadius;
if iter == maxIterations
    fprintf('Maximum iterations reached without convergence.\n');
end

%% 
figure;
subplot(2, 2, 1); imshow(img_gray); title('Gray Image');
subplot(2, 2, 2); imshow(sobel_imgs); title('Sobel Filter');
subplot(2, 2, 3); imshow(radial_mask); title('Radial Mask');
subplot(2, 2, 4); imshow(filter_edges); title('Filter Edges');
figure;
imshow(img_gray);
hold on;
viscircles(center, radius, 'edgecolor', 'r','linewidth',0.2);
hold on;
plot(center(1),center(2), '+', 'color', 'red'); 
hold off;

houghTransform()函数：

function [newCenter, newRadius, hough_space] = houghTransform(edgeImage, radii, xRange, yRange)
    % 初始化霍夫空间
    hough_space = zeros(length(yRange), length(xRange), length(radii));
    % 对边缘图像中的每个点进行处理
    [edgeRows, edgeCols] = find(edgeImage);
    for idx = 1:length(edgeRows)
        edgeX = edgeCols(idx);
        edgeY = edgeRows(idx);
        
        for rIdx = 1:length(radii)
            r = radii(rIdx);
            for xIdx = 1:length(xRange)
                x0 = xRange(xIdx);
                for yIdx = 1:length(yRange)
                    y0 = yRange(yIdx);
                    % 计算该点到理论圆心的距离
                    squareDistance = (edgeX - x0)^2 + (edgeY - y0)^2;
                    % 当点位于圆的边缘上(考虑一定误差)
                    if abs(squareDistance - r^2) <= (r * 0.1)^2
                        hough_space(yIdx, xIdx, rIdx) = hough_space(yIdx, xIdx, rIdx) + 1;
                    end
                end
            end
        end
    end

    % 找到霍夫空间中的最大值
    [~, maxIndex] = max(hough_space(:));
    [yIndex, xIndex, rIndex] = ind2sub(size(hough_space), maxIndex);
    newCenter = [xRange(xIndex), yRange(yIndex)];
    newRadius = radii(rIndex);
end

标签：gray,高斯,img,天体,边缘,圆心,图像处理,拟合,半径
From： https://blog.csdn.net/qq_39631801/article/details/142304598