Add files via upload

Author: Rondapapa

bezier.m

@@ -0,0 +1,28 @@
+function Y=bezier(y)
+%Bezier Curve 
+% figure;
+% plot(1:length(y),y,'color',[0.1 0.2 0.5],'marker','o');
+% hold on;
+x=1:length(y);
+n=length(y);
+t=0:0.001:1;
+A=0;B=0;
+for k=0:n-1
+    tmp=nchoosek(n-1,k)*t.^k.*(1-t).^(n-1-k);
+    A=A+tmp*x(k+1);
+    B=B+tmp*y(k+1);
+end
+
+%figure;
+% plot(A,B,'color',[0.2 0.5 0.1]);
+Y=zeros(1,length(y));
+Y(1)=y(1);
+for i=2:length(y)-1
+    position=round(length(t)/(length(y)-1)*(i-1));
+    Y(i)=B(position);
+end
+Y(end)=y(end);
+
+% plot(1:length(y),Y,'color',[0.7 0.2 0.1],'marker','+');
+% xlabel('x');
+% ylabel('y');

cosSimu.m

@@ -0,0 +1,26 @@
+function result=cosSimu(a,b)
+% Loop Cosine Similarity
+n=fix(size(a,2)/200)+1; % 分成 n 段
+n=1;
+temp_a=cell(1,n);    %  每段成为一个数组
+temp_a_node=zeros(1,n-1);
+for i=1:n-1
+    temp_a_node(i)=fix(length(a)/n)*i;
+end
+temp_a_node=[1 temp_a_node length(a)];
+for i=2:n+1
+    temp_a{i-1}=a(temp_a_node(i-1):temp_a_node(i));
+end
+
+overPoint=length(b)-1;
+b=[b b]; % 双倍，用于滚动搜索
+
+coeff=zeros(n,overPoint);
+for i=1:n
+    for j=1:overPoint
+        temp_b=b(j:j+length(temp_a{i})-1);  
+        coeff(i,j)=dot(temp_a{i},temp_b)/(norm(temp_a{i})*norm(temp_b));
+    end
+end
+coeff=max(coeff');
+result=mean(coeff);

grayGrade.m

@@ -0,0 +1,19 @@
+function [grade,totalNetValue]=grayGrade(image,netWide,netHeight)
+totalNetValue=zeros(fix(size(image,1)/netHeight),fix(size(image,2)/netWide)); % 网格汇总
+for i=1:size(image,1)/netHeight
+    for j=1:size(image,2)/netWide
+        temp_block=image((i-1)*netHeight+1:i*netHeight,(j-1)*netWide+1:j*netWide);
+        [r,~]=find(temp_block<200);
+        totalNetValue(i,j)=length(r);
+    end
+end
+totalNetValue=sort(totalNetValue(:));
+totalNetValue=totalNetValue/(netWide*netHeight);
+
+node=0.15; % 建立节点，将 0 到 1 分成 2 大段，前段的筛孔变化较精细
+grade=zeros(1,round(node/0.001)); % 确定初始的级配筛孔
+for i=1:round(node/0.001)
+    [r,~]=find(totalNetValue<i*0.001);
+    grade(i)=length(r);
+end
+grade=grade/length(totalNetValue);

mainProcess.m

@@ -0,0 +1,101 @@
+clc
+clear;
+warning('off');
+%% 第一步 图片预处理
+I=cell(1,2);
+I{1} = imread('realTraceMap.bmp');
+I{2} = imread('simulatedTraceMap1.bmp');
+%figure;
+figure('Name','Real Trace Map');
+imshow(I{1});
+%figure;
+figure('Name','Simulated Trace Map');
+imshow(I{2});
+totalgray_a=totalgray(I{1},'r');
+totalgray_b=totalgray(I{2},'b');
+yita_totalgray=1-dist(totalgray_a,totalgray_b)/totalgray_a  % 灰度相关的百分比
+
+%% 第三步 计算灰度级配曲线
+netWide=size(I{1},2)/16;
+netHeight=size(I{1},1)/16;
+
+[grade_a,totalNetValue_a]=grayGrade(I{1},netWide,netHeight);
+[grade_b,totalNetValue_b]=grayGrade(I{2},netWide,netHeight);
+
+x=1:length(grade_a);
+
+
+figure('Name','Gray Grade Curve');
+plot(x,grade_a,'r--');
+hold on;
+plot(x,grade_b,'b');
+set(gca,'ylim',[0 1]);
+set(gca,'yticklabel',{'0%','10%', '20%', '30%', '40%' ,'50%','60%','70%','80%','90%','100%'});
+Mhtdist=zeros(1,length(grade_a)); %计算 曼哈顿距离
+for i=1:length(grade_a)
+    Mhtdist(i)=dist(grade_a(i),grade_b(i));
+end
+stem(1:length(grade_a),Mhtdist,'Marker','none');
+xlabel('screen');
+ylabel('proportion');
+legend('gray grade curve of real trace map','gray grade curve of simulation trace map', 'Location','Best');
+%yita_grade=pdist2(grade_a,grade_b,'cityblock')/length(grade_a);
+yita_graygrade=riter(Mhtdist)
+
+%% 第四步 以按比例积分的方式对 Radon变换 改进，搜索特征方向、方差及曲线
+[varians_a,curves_a,fdegrees_a]=radonTrs(I{1});
+[varians_b,curves_b,fdegrees_b]=radonTrs(I{2});
+% 密度曲线方差的比较柱状图
+M=zeros(length(curves_a),141);   % 画拉东变换 的图
+for i=1:length(curves_a)
+    M(i,1:length(curves_a{i}))=curves_a{i};
+end 
+
+figure('Name','Characteristic Directions');
+x=(1:length(varians_a))*1;
+h=bar(x,[varians_a' varians_b']);
+xlabel('\theta');
+ylabel('variance');
+legend('real trace map','simulation trace map','Location','Best');
+set(gca,'xticklabel',{'15','30', '45', '60' ,'75','90','105' ,'120', '135','150','165','180'});
+
+
+yita_direction=yitaDirection(fdegrees_a,fdegrees_b)
+%% 第五步 以 调整的余弦相似度 方法，沿特征方向滚动匹配两组图像
+fIndex=round(fdegrees_a/15); 
+% figure;
+
+for i=1:4
+    result(i)=cosSimu(curves_a{fIndex(i)},curves_b{fIndex(i)});
+    x=1:length(curves_a{fIndex(i)});
+%     plot(x,curves_a{fIndex(i)}+i-0.94,'r');
+%     hold on;
+%     plot(x,curves_b{fIndex(i)}+i-0.94,'b');
+%     hold on
+end
+yita_loopcosSimu=mean(result)
+% 
+% xlabel('Cross section');
+% legend('Gray Density Distribution Curve (Real Trace Map)',...
+%     'Gray Density Distribution Curve (Simulated Trace Map)', 'Location','Best');
+% set(gca,'ylim',[-0.1,4]);
+% 
+% set(gca,'yticklabel',{strcat('\theta=',num2str(fdegrees_a(1)),',0.0'),'0.5',...
+%     strcat('\theta=',num2str(fdegrees_a(2)),',0.0'),'0.5',...
+%     strcat('\theta=',num2str(fdegrees_a(3)),',0.0'),'0.5',...
+%     strcat('\theta=',num2str(fdegrees_a(4)),',0.0'),'0.5',' '});
+% axes;
+% 
+% reresult=0.99-result;
+% bh=barh(0:3,[result',reresult'],'stack','BarWidth',0.05,'LineStyle','none');%'FaceColor',[0.423,0.251,0.392],
+% set(bh(1),'facecolor',[0.635294117647059,0.0784313725490196,0.184313725490196]);
+% set(bh(2),'facecolor',[0.854901960784314,0.701960784313725,1]);
+% set(gca,'xaxislocation','top','color','none');
+% set(gca,'ylim',[-0.1 4],'xlim',[0,1],'XTick',0:0.2:1,'ytick',[]);
+% ylabel('Gray Density');
+% xlabel('Loop Cosine Similarity');
+% legend('Loop Cosine Similarity');
+%% 综合相似度
+comprehensive_Similarity=0.2*(yita_totalgray+yita_graygrade)+0.3*(yita_loopcosSimu+yita_direction)
+
+

radonTrs.m

@@ -0,0 +1,47 @@
+function [varians,curves,fdegrees]=radonTrs(I)
+% density Radon transform
+IBack=zeros(size(I)); % 做背景矩阵,与 I 同时旋转，叠加计算
+limit=min([size(I,2) size(I,1)]);
+
+step=15;
+varians=zeros(1,180/step);
+curves=cell(1,fix(180/step));
+sliceWide=5;   % 切片的宽度， 最少为迹线的宽度
+for i=step:step:180   % 减 1 是因为 180 度和 0 度是一样的
+    tempI=~I; tempIBack=~IBack;
+    tempI=imrotate(tempI,i); tempIBack=imrotate(tempIBack,i);
+    tempI=~tempI;tempIBack=~tempIBack;
+    % 开始切片计算
+    sliceNum=fix(size(tempI,2)/sliceWide);  % 切片的数量
+    sliceDensity=zeros(1,sliceNum);     %  切片的相对灰度密度
+    
+    tmp=zeros(1,sliceNum); % 消除边缘影响用
+    for j=1:sliceNum
+        SliceI=tempI(:,(j-1)*sliceWide+1:j*sliceWide);
+        SliceIBack=tempIBack(:,(j-1)*sliceWide+1:j*sliceWide);
+        [r1,~]=find(SliceI==0);
+        [r2,~]=find(SliceIBack==0);
+        if length(r2)>0 
+         sliceDensity(j)=length(r1)/length(r2);
+        end
+        tmp(j)=length(r2);
+    end
+    % 密度曲线 消除边缘影响
+    leftnode=0;
+    tmpsum=0;
+    for k=1:sliceNum
+        if tmpsum>limit*sliceWide
+            leftnode=k;
+            break;
+        end
+        tmpsum=tmpsum+tmp(k);
+        
+    end
+    curves{i/step}=sliceDensity;
+    varians(i/step)=std(sliceDensity);
+end
+%figure;
+x=[1:length(varians)]*step;
+
+[~, fdegrees]=sort(varians,'descend');
+fdegrees=fdegrees(1:4)*step;%

riter.m

@@ -0,0 +1,18 @@
+function [yitagrade]=riter(Mhtdist)
+xx=0:0.001:1;
+n=zeros(1,length(xx));
+index=0;
+tmp=1;
+for i=1:length(xx)
+    n(i)=1-length(find(Mhtdist>=1-xx(i)))/73;
+    
+    if abs(n(i)-xx(i))<tmp
+        tmp=abs(xx(i)-n(i));
+        index=i;
+    end
+end
+yitagrade=xx(index);
+figure;
+plot(xx,n,'r');
+hold on;
+plot(xx,xx,'b');

totalgray.m

@@ -0,0 +1,31 @@
+function result=totalgray(I,fig)
+% fig=='r',输出红色图像
+% fig=='b'，输出蓝色图像
+% 其它，不输入图像
+
+[r,~]=find(I<200);
+% 获得个数,总的迹线像素点的个数
+result=length(r);
+sliceWide=5;
+sliceNum=fix(size(I,2)/sliceWide);
+sliceDensity=zeros(1,sliceNum);
+for i=1:sliceNum
+    Slice=I(:,(i-1)*sliceWide+1:i*sliceWide);
+    [r,~]=find(Slice==0);
+    sliceDensity(i)=length(r);
+end
+x=1:sliceNum;
+if strcmp('r',fig)
+    figure;
+    bar(x,sliceDensity,'r');
+    set(gca,'ylim',[0 100]);
+    xlabel('x');
+    ylabel('gray statistic');
+end
+if strcmp('b',fig)
+    figure;
+    bar(x,sliceDensity,'b');
+    set(gca,'ylim',[0 100]);
+    xlabel('x');
+    ylabel('gray statistic');
+end

yitaDirection.m

@@ -0,0 +1,27 @@
+function yita_direction=yitaDirection(fdegrees_a,fdegrees_b)
+
+p=[];q=[];
+for i=1:length(fdegrees_a)
+    for j=1:length(fdegrees_a)
+        if fdegrees_a(i)==fdegrees_b(j)
+            p(i,j)=1;
+            q=[q;i j];
+        else
+            p(i,j)=0;
+        end
+    end
+end
+
+count=[];k=1;
+for i=1:length(q)-1
+    for j=i+1:length(q)
+        if (q(j,1)-q(i,1))*(q(j,2)-q(i,2))>0
+            count(k)=1;
+        else
+            count(k)=0;
+            
+        end
+        k=k+1;
+    end
+end
+yita_direction=length(q)/length(fdegrees_a)*0.8+sum(count)/length(count)*0.2;