【数学建模】CUMCM-2014A 嫦娥三号软着陆过程 避障阶段图像处理

实现

图像高斯滤波→求梯度图→二次高斯滤波→二值化→开运算、膨胀→面积带通滤波→求霍夫圆进行优化→搜索定位。

1.图像模糊

也就是对图像进行高斯滤波。在我们得到高程图后，它还存在着不少孤立噪声点，因此我们需要对其进行滤波。我们选用的是高斯滤波。
高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到。具体操作是：用一个掩模扫描图像中的每一个像素，以邻域内像素的加权平均灰度值去替代模板中心像素点的值。

2.求梯度

由于需要寻找平坦的合适降落的区域，因此我们对2300m*2300m的高程图求梯度，我们认为梯度越大越不平坦，梯度越小越平坦。我们使用sobel算子求取梯度。
Sobel算子的原理是，使用两个3*3矩阵（如下图）和原始图片进行卷积，分别得到横向和纵向的梯度值，将其取绝对值并相加就能得到梯度图。

得到的梯度图如下所示：

3.再次图像模糊

由于仅仅依靠梯度图得到的效果不是特别好，我们还需要做其他处理。因此我们再做一次图像模糊。

4.图像二值化

原始图像经过模糊后，我们需要确定一个阈值并对图像进行二值化，将高于阈值的像点变成1，低于阈值的像素点变成0，从而图像变为黑白的二值图。关于阈值的选取有不少方法，比如全局平均值，即sum(所有点)/2300^2。尽管阈值的选取对结果有一定的影响，但不会很大，因此暂不考虑复杂的情况。我们一开始使用的是matlab自带的根据灰度图选取阈值的函数，但是觉得效果不佳，因此稍微做了人工调整，定下了阈值。黑白二值图如下图所示：

5.开闭运算、膨胀、腐蚀

尽管之前进行了滤波，不过仍旧存在着两个问题，其一是还是有许多零散的点分布在图中，需要进一步降噪；其二是由于图像本身的像素不够高，导致前两步处理之后有些陨石坑的边缘没了，导致它不够完整，航天器如果误判可能会飞到陨石坑上，因此会有潜在的危险。综合这两个问题，我们需要对图像进行开运算和闭运算。
其中开运算是先腐蚀后膨胀，它完全删除了不能包含结构元素的对象区域，平滑了对象的轮廓，断开了狭窄的连接，去掉了细小的突出部分。因此使用开运算能一定程度上解决第一个问题。闭运算是先膨胀后腐蚀，它会将狭窄的缺口连接起来形成长的弯口，并填充比结构元素小的洞，即将区域联通，因此它能一定程度上解决第二个问题，不过鉴于膨胀的效果更为直接，因此我们用膨胀来联通区域。
至于腐蚀和膨胀的次数问题，需要根据实际图像来确定，我们选取的原则类似Neyman-Pearson原则，在控制第二类问题带来的误差小于某个值的情况下尽可能减小第一类问题带来的误差。即优先保证安全性(毕竟零散的野点其实影响没有很大)。
我们对前面得到的图先进行开运算，然后膨胀，得到如下的图：

6.联通区域面积的带通滤波

相当于利用面积进行区域的筛选。我们求取联通区域的面积，将过大的和过小的都滤掉，留下来的是主要的陨石坑区域。其中过大区域是陨石坑外围区域，包括平坦区域，过小区域是前面操作后仍旧遗留下来的零散的点域，经过这样的带通滤波，我们可以得到可降落区域，如下图白色区域所示：

7.霍夫圆优化

由于较深的陨石坑能够在梯度图中得到显著的表达，但是较浅的陨石坑就不够明显了，因此需要在前面的基础上进行一些纠正。首先我们考虑到陨石坑的二值化图像形态是一圈一个白色的圆形（凸起的山）中间有一个黑色的小圆（凹进去的坑），即便是较浅的陨石坑也有这样的特征。因此我们用求取中间黑色区域的霍夫圆，并将其膨胀，使得其能涵盖它所在的不平坦区域，从而将梯度图无法涵盖到的不平坦区域涵盖了，提高了安全性。可以从下图中看到，它明显更加符合月球的地貌高程图：

代码http://cn.mathworks.com/help/images/examples/identifying-round-objects.html?s_eid=PSM_10906

8.搜索算法

我们搜索的原则是（从上到下按重要度降序排列）：
1） 保证安全，尽量不落入不平坦区域
2） 保证降落区域尽可能大，不过不必特别大，能够保证降落精度范围内即可
3） 尽可能靠近图像中心，因为图像中心是飞行器此时所在位置的投影点，如果降落在边上，需要进行额外的姿态调整，增加了风险和燃料损耗。
那么，基于上面得到的图像，我们考虑了两种搜索的方法，其一是将图像按照飞行器的尺寸为基元进行矩形离散化(分割成多个小矩形)，将其以放射型扩大，拼接起来，选择靠近中心且足够大的矩形，作为降落点；其二是以安全半径螺旋式搜索，该算法时间上有优势，但编程起来有点麻烦。

然后精确避障阶段和此类似。

反思

可能也被一篇参考文章影响的，感觉自己的思路比较局限，其实总感觉（取了梯度就差不多了），还有就是，作图方式比较局限= =只有平面图(这么一看还真是……就我图是黑白的 哇的一声)。

像优秀论文的评价就是，
粗避障阶段考虑了初始的运动状态（而我们在做第一题的时候由于太垃圾求不出解= =）
然后两个阶段都是直接分块然后上评估函数。
完全是不一样的思路了= =

像在这个过程中我想的主要是（要怎么把月球的特征用起来）然后想到了霍夫圆，并且不知道为什么（觉得不应该停在坑里面）然后就把里面都去掉了=并且我们的话是从中间开始向外搜索

（不过我做这题的时候队友一直表示反正CUMCM里图像处理题不太有= =）

代码

clc;clear;
img = imread('1.tif');

gausFilter = fspecial('gaussian',[25,25],100);
img1=imfilter(img,gausFilter,'replicate');

%%
%求梯度
H1=[-1 -2 -1; 0 0 0; 1 2 1];
H2 = H1';
img21 = filter2(H1,img1);
img22 = filter2(H2,img1);
img2 = (abs(img21)+abs(img22));
[sizex,sizey] = size(img2);
img2 = img2(2:sizex - 1,2:sizey - 1);
img2 = img2 / max(max(img2));
%%
%降噪
gausFilter = fspecial('gaussian',[50,50],100);
img2=imfilter(img2,gausFilter,'replicate');
imshow(img2);

%%
%二值化


thresh = graythresh(img2);     %自动确定二值化阈值
img3 = im2bw(img2,0.17);       %对图像二值化
imshow(img3);

 B=[0 1 0; 1 1 1;  0 1 0];


se = strel('disk',4);
%img3 = bwareaopen(img3,100);
img3 = imopen(img3,se);

 for k = 1:20
     img3=imdilate(img3,B);%图像A1被结构元素B膨胀
 end
 %img3 = imclose(img3,se);
imshow(img3);


%%
%填充联通区域面积过大|过小的
img4 = img3==0;  % 反向
L = bwlabeln(img4);
S = regionprops(L, 'Area');
bw2 = ismember(L, find([S.Area] >= 1000)) & ismember(L, find([S.Area] <= 100000));
%bw2 = ismember(L, find([S.Area] <= 100000));

imshow(bw2);

%%
%求霍夫圆
[B,L] = bwboundaries(bw2,'noholes');

% Display the label matrix and draw each boundary
imshow(label2rgb(L, @jet, [.5 .5 .5]))
hold on
for k = 1:length(B)
  boundary = B{k};
  plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)
end

stats = regionprops(L,'Area','Centroid');

threshold = 0.5;
    r1 = [] ;
    D = [];
% loop over the boundaries
for k = 1:length(B)

  % obtain (X,Y) boundary coordinates corresponding to label 'k'
  boundary = B{k};

  % compute a simple estimate of the object's perimeter
  delta_sq = diff(boundary).^2;    
  perimeter = sum(sqrt(sum(delta_sq,2)));

  % obtain the area calculation corresponding to label 'k'
  area = stats(k).Area;

  % compute the roundness metric
  metric = 4*pi*area/perimeter^2;
  r = perimeter /(2*pi);
  r1 = [r1 , r];

  % display the results
  metric_string = sprintf('%2.2f',metric);

  % mark objects above the threshold with a black circle
  if metric > threshold
    centroid = stats(k).Centroid;
    plot(centroid(1),centroid(2),'ko');
    D = [D,k];
    rectangle('Position',[centroid(1)-r,centroid(2)-r,2*r,2*r],'Curvature',[1,1]),axis equal
  end

  text(boundary(1,2)-35,boundary(1,1)+13,metric_string,'Color','y',...
       'FontSize',14,'FontWeight','bold');

end

title(['Metrics closer to 1 indicate that ',...
       'the object is approximately round']);

 %%
 %显示区域
  L = bwlabeln(img4);
S = regionprops(L, 'Area');
img5 = ismember(L, find([S.Area] >= 100000)) ;
imshow(img5);
hold on
alpha=0:pi/20:2*pi;    %角度[0,2*pi] 
for m = 1:length(D)
    k = D(m);
centroid = stats(k).Centroid;
r = r1(k)*5;
%rectangle('Position',[centroid(1)-r,centroid(2)-r,2*r,2*r],'Curvature',[1,1]),axis equal

 R=r;                   %半径 
 i = floor(R*cos(alpha)+centroid(1)); 
 j = floor(R*sin(alpha)+centroid(2)); 
 for n=1:length(i)
     x = i(n);
     y = j(n);
     if(x>0 && y>0 && (x < sizex - 2) && (y < sizey - 2))
         img5(x,y)=0;  
     end
end
end
imshow(img5);