基于内容的遥感影像场景检索
基于内容的遥感影像场景检索（一）
（一）基础练习
1、从遥感影像文件夹里批量读取遥感影像，并计算其灰度直方图特征
①首先学习matlab批量读取文件
（由于之前此部分没认真学习，所以试图理解每一个函数并进行简单整理）
参考：https://blog.csdn.net/u012273127/article/details/76150861
filePath1 = 'D:\data\RS_CBIR\RS_Dataset\'; %存储图像的路径
fileExt = '*.tif'; %待读取图像的后缀名
%dir:读取一个文件夹下的所有的文件，并存储到一个结构体中。
%在结构体中存储了该文件夹下所有的文件的名字以及文件的创建日期。
files = dir(fullfile(filePath1,fileExt));
%size（data, n）返回多维矩阵的第n维的维度，
len1 = size(files,1);%本来就只有一维，其实也就是通过size返回了文件中图像的个数
%遍历路径下每一幅图像
for i=1:len1
fileName = strcat(filePath1,files(i).name); %strcat：字符串附加
image = imread(fileName);%读取图像数据
imshow(image);%一张张显示图片
end
【补充1】：字符串与字符数组
不常用：char 表示字符数组，将短文本片段存储为字符向量，如 c = 'Hello World'。
常用：string表示字符串，使用双引号创建字符串，例如 str = "Greetings friend"
参考：https://ww2.mathworks.cn/help/matlab/characters-and-strings.html
②计算灰度直方图特征：
由于暂时不清楚灰度直方图特征有哪些，所以待补充
补充：了解到暂时灰度直方图特征在此仅计算均值以及方差即可，开始查资料
参考：https://blog.csdn.net/yffly0406/article/details/80228103?biz_id=102&utm_term=matlab%E8%AE%A1%E7%AE%97%E7%81%B0%E5%BA%A6%E7%9B%B4%E6%96%B9%E5%9B%BE%E5%9D%87%E5%80%BC&utm_medium=distribute.pc_search_result.none-task-blog-2allsobaiduweb~default-4-80228103&spm=1018.2118.3001.4449
计算此两项还是较为简单的，直接把更改后的for循环列如下
for i=1:len1
fileName = strcat(filePath1,files(i).name); %strcat：字符串附加
image = imread(fileName);%读取图像数据
%imshow(image);一张张显示图片
GrayImage=rgb2gray(image);%将彩色图转为灰度图
GrayImage=double(GrayImage); %将uint8型转换为double型，否则不能计算统计量
[m,n]=size(GrayImage);%得到矩阵的行数以及列数
avg=mean2(i); %求图像均值
a1=var(GrayImage(:));%求图像方差
end
2、随机选一张影像作为检索参考影像，从遥感影像文件夹里
（例如，通过灰度直方图特征的相似性对文件夹里遥感影像进行相似性排序，其中灰度直方图特征的相似性可以考虑使用特征向量差的范数来计算）
①首先不着急编代码和面向百度编程，先想一想有没有思路
去看遥感文件夹，发现分为六种图像的变换，所以我们要做的
首先是随机选取一张影像作为参考。
随后根据提示，计算图片灰度直方图的特征向量差的范数，
然后把所有图片的范数进行排序，
最后取合适的区间，若在区间内则算是视觉内容相似的遥感影像。
②把思路理清晰之后，我们可以开始尝试一步步完成。首先完成第一步，取随机图像，我们想想只需要在第一问中变量len1表示的图片个数范围内取随机数即可
这一步还是非常简单的，经过网上搜索，取1，n之间的整数随机数即用 randi 函数
参考：https://ww2.mathworks.cn/help/matlab/random-number-generation.html
x=randi(len1,1,1);
%以下是测试程序
fileName = strcat(filePath1,files(x).name);
image = imread(fileName);
imshow(image);
③随后就到了最关键的一步：计算图片灰度直方图的特征向量差的范数。由于啥都不懂，直接开始查资料，据参考，此处求范数和求余弦相似度都可以，找到了余弦相似度有关的博客，试着进行修改尝试
参考：https://blog.csdn.net/sinat_28826891/article/details/95310592?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160627066619726890300947%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=160627066619726890300947&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allfirst_rank_v2~rank_v28-7-95310592.pc_first_rank_v2_rank_v28&utm_term=matlab%E8%BF%9B%E8%A1%8C%E4%BD%99%E5%BC%A6%E7%9B%B8%E4%BC%BC%E5%BA%A6&spm=1018.2118.3001.4449
%求随机图像并显示
x=randi(len1,1,1);
fileName = strcat(filePath1,files(x).name); %strcat：字符串附加
imageR =imread(fileName);%读取图像数据
imageRR=imageR;%这一步的意义：在于如果不经过这样的转换，后面显示图像的时候会突然转成灰度图
figure;
imshow(imageRR);%显示所抽取的原始图像
imageR=rgb2gray(imageR);
vC=zeros(1,len1);%将算得的余弦距离放入其中
for k=1:len1
fileName2 = strcat(filePath1,files(k).name);
imageC=imread(fileName2);
imageC=rgb2gray(imageC);
t1=imageR;
[a1,b1]=size(t1);
t2=imageC;
t2=imresize(t2,[a1 b1],'bicubic');%将两张图片缩放为一致大小
t1=round(t1);%取四舍五入
t2=round(t2);
e1=zeros(1,256);%生成零矩阵
e2=zeros(1,256);
%获取直方图分布/
for i=1:a1
for j=1:b1
m1=t1(i,j)+1;
m2=t2(i,j)+1;
if(m2==257)
m2=256;
end
e1(m1)=e1(m1)+1;%读入灰度值
e2(m2)=e2(m2)+1;
end
end
%将直方图分为64个区
m1=zeros(1,64);
m2=zeros(1,64);
for i=0:63
m1(1,i+1)=e1(4*i+1)+e1(4*i+2)+e1(4*i+3)+e1(4*i+4);
m2(1,i+1)=e2(4*i+1)+e2(4*i+2)+e2(4*i+3)+e2(4*i+4);
end
%计算余弦相似度
A=sqrt(sum(sum(m1.^2)));
B=sqrt(sum(sum(m2.^2)));
C=sum(sum(m1.*m2));
cos1=C/(A*B);%计算余弦值
cos2=acos(cos1);%弧度
v=cos2*180/pi;%换算成角度
vC(1,k)=v;
end
④把所有图片的余弦进行排序，取最小的五个（现有min函数是取最小值，则每取到一个最小值就接着把其换成999，下一步取则取下一个值，依次推）
for i=1:5 %这里取的是排序前五的最小
ind = find(vC==min(min(vC)));%取一维矩阵中最小值
fileName = strcat(filePath1,files(ind).name); %strcat：字符串附加
imageMin =imread(fileName);%读取图像数据
figure;
imshow(imageMin);
vC(1,ind)=999;
end
3、通过多种Matlab的多种可视化手段来展示算法的中间结果
（例如，可视化展示图象的直方图特征，展示检索图象特征与文件夹中图象特征之间相似性数值结果）
先理一下思路：其实就是对2进行拓展，展示图像直方图特征暂时可以理解为显示图像直方图，然后最好的办法就是两个图片放一块（原图+直方图）。显示特征相似性数值计算结果，这个可以暂时理解成直接输出cos值，即在后续图像上加小标题。于是根据本思路开始实现。
①显示两个图片：
首先经过了解可以用Hold on;来将两个图片放到一起，但是我在尝试之后的结果竟然变成了两个图片重叠在一起。后经了解可以用subplot同时显示多个图像。
参考：https://blog.csdn.net/u013925378/article/details/82878407?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160634980719724848153040%2522%252C%2522scm%2522%253A%252220140713.130102334…%2522%257D&request_id=160634980719724848153040&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2alltop_click~default-1-82878407.pc_first_rank_v2_rank_v28&utm_term=matlab%E6%98%BE%E7%A4%BA%E5%A4%9A%E4%B8%AA%E5%9B%BE%E5%83%8F&spm=1018.2118.3001.4449
https://baike.baidu.com/item/subplot/7801858?fr=aladdin
②加小标题&cos值
for i=1:5
ind = find(vC==min(min(vC)));
fileName = strcat(filePath1,files(ind).name); %strcat：字符串附加
imageMin =imread(fileName);%读取图像数据
imageMinn=imageMin;
figure;
subplot(1,2,1);
imshow(imageMin);
j=i;
title(['近似图',num2str(j)]);
hold on;
imageMin=rgb2gray(imageMin);
subplot(1,2,2);
imhist(uint8(imageMin));
title(['余弦夹角为',num2str(vC(1,ind))]);
vC(1,ind)=999;
end
4、将检索参考影像和与其相似性排在前五的遥感影像通过可视化的手段显示出来**
（例如，可以考虑使用subplot函数功能）；
在②中已经实现，故在此不再赘述。
（二）高级练习
1，遥感影像的多特征表达方式（例如，尝试利用灰度直方图、灰度共生矩阵、梯度方向直方图、深度特征等方式来表达遥感影像的视觉内容）；
①灰度直方图：在前面已经实现，故不再进行阐述
②灰度共生矩阵：灰度共生矩阵实现相当简单，从网上了解后发现已经有函数可以实现
参考：https://blog.csdn.net/lang_yubo/article/details/79151132?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160646788719215668896710%2522%252C%2522scm%2522%253A%252220140713.130102334…%2522%257D&request_id=160646788719215668896710&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2allbaidu_landing_v2~default-3-79151132.pc_first_rank_v2_rank_v28&utm_term=matlab%E6%98%BE%E7%A4%BA%E7%81%B0%E5%BA%A6%E5%85%B1%E7%94%9F%E7%9F%A9%E9%98%B5&spm=1018.2118.3001.4449
glcm = graycomatrix(imageR,'Offset',[0 2]);
set(0,'defaultFigurePosition',[100,100,1000,500]);
set(0,'defaultFigureColor',[1 1 1]);
③梯度方向直方图：
这个比较复杂，我其实概念也没有理解的太清楚，还需要再理解一下。经过计算并实现3D视图可视化之后，代码如下：
参考：https://blog.csdn.net/zhengtu009/article/details/46652275?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522160635396919721942218409%2522%252C%2522scm%2522%253A%252220140713.130102334…%2522%257D&request_id=160635396919721942218409&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2alltop_click~default-2-46652275.pc_first_rank_v2_rank_v28&utm_term=matlab%E6%A2%AF%E5%BA%A6%E6%96%B9%E5%90%91%E7%9B%B4%E6%96%B9%E5%9B%BE&spm=1018.2118.3001.4449
imageRRR=sqrt(imageRRR); %伽马校正
%下面是求边缘
fy=[-1 0 1]; %定义竖直模板
fx=fy'; %定义水平模板
Iy=imfilter(imageRRR,fy,'replicate'); %竖直边缘
Ix=imfilter(imageRRR,fx,'replicate'); %水平边缘
Ied=sqrt(Ix.^2+Iy.^2); %边缘强度
Iphase=Iy./Ix; %边缘斜率，有些为inf,-inf,nan，其中nan需要再处理一下
%下面是求cell
step=16; %step*step个像素作为一个单元
orient=9; %方向直方图的方向个数
jiao=360/orient; %每个方向包含的角度数
Cell=cell(1,1); %所有的角度直方图,cell是可以动态增加的，所以先设了一个
ii=1;
jj=1;
for i=1:step:m %如果处理的m/step不是整数，最好是i=1:step:m-step
ii=1;
for j=1:step:n %注释同上
tmpx=Ix(i:i+step-1,j:j+step-1);
tmped=Ied(i:i+step-1,j:j+step-1);
tmped=tmped/sum(sum(tmped)); %局部边缘强度归一化
tmpphase=Iphase(i:i+step-1,j:j+step-1);
Hist=zeros(1,orient); %当前step*step像素块统计角度直方图,就是cell
for p=1:step
for q=1:step
if isnan(tmpphase(p,q))==1 %0/0会得到nan，如果像素是nan，重设为0
tmpphase(p,q)=0;
end
ang=atan(tmpphase(p,q)); %atan求的是[-90 90]度之间
ang=mod(ang*180/pi,360); %全部变正，-90变270
if tmpx(p,q)<0 %根据x方向确定真正的角度
if ang<90 %如果是第一象限
ang=ang+180; %移到第三象限
end
if ang>270 %如果是第四象限
ang=ang-180; %移到第二象限
end
end
ang=ang+0.0000001; %防止ang为0
Hist(ceil(ang/jiao))=Hist(ceil(ang/jiao))+tmped(p,q); %ceil向上取整，使用边缘强度加权
end
end
Hist=Hist/sum(Hist); %方向直方图归一化
Cell{ii,jj}=Hist; %放入Cell中
ii=ii+1; %针对Cell的y坐标循环变量
end
jj=jj+1; %针对Cell的x坐标循环变量
end
%下面是求feature,2*2个cell合成一个block,没有显式的求block
[m n]=size(Cell);
feature=cell(1,(m-1)*(n-1));
for i=1:m-1
for j=1:n-1
f=[];
f=[f Cell{i,j}(:)' Cell{i,j+1}(:)' Cell{i+1,j}(:)' Cell{i+1,j+1}(:)'];
feature{(i-1)*(n-1)+j}=f;
end
end
l=length(feature);
f=[];
for i=1:l
f=[f;feature{i}(:)'];
end
mesh(f);
title("梯度方向直方图");
④深度特征：暂时没有搞清楚需要做什么，故咨询以后再进行相关补充
成果如下：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-npoyrWDw-1606698041727)(C:\Users\admin\AppData\Roaming\Typora\typora-user-images\image-20201128111245280.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-aTjzr7db-1606698041730)(C:\Users\admin\AppData\Roaming\Typora\typora-user-images\image-20201128111448531.png)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0iFWs97H-1606698041732)(C:\Users\admin\AppData\Roaming\Typora\typora-user-images\image-20201128111614034.png)]
1)+j}=f;
end
end
l=length(feature);
f=[];
for i=1:l
f=[f;feature{i}(????’];
end
mesh(f);
title(“梯度方向直方图”);
④深度特征：暂时没有搞清楚需要做什么，故咨询以后再进行相关补充
成果如下：
[外链图片转存中...(img-npoyrWDw-1606698041727)]
[外链图片转存中...(img-aTjzr7db-1606698041730)]
[外链图片转存中...(img-0iFWs97H-1606698041732)]
上一篇： Halcon 轮廓合并算子