harries算子特征提取与图像匹配

概述：图像处理过程中经常要遇到图像匹配的问题，通过灰度相关法可以解决比较简单的平移、旋转、亮度差异匹配的问题，当平移角度过大，灰度差异过大、旋转角度较大时，灰度相关法会变得比较麻烦，采用角点匹配的方法能取得不错的成果，高级的有sift和suft，这里我们介绍一种比较简单的角点检测方法harris
，对harris原理上的东西不多讲，主要谈一下图像匹配注意点，以及我自己的处理方式，方式是以前学习时自己想的，没有查阅相关资料，仅供参考。
一、Harris角点检测
注意点1，非极大值抑制，采用阈值分割的方法得到的角点往往都扎堆存在，进行非极大值抑制，得到最合适的角点，如图。

注意点2，边界处理，由于原始图像可能是经过拼接后的图像，形状不规则，其边缘可能不是矩形，在边界上非常容易产生角点，这些角点都要除去，如下图绿色圈内的点

注意点3，但是由于先验知识，我们可以知道图像的相对位置关系，比如一张在左一张在右 ，在选取角点 时，对于左图我们可以只保留靠近右边的角点。右下图为删除靠左角点后的，这两张图是在第12展图匹配拼接后求的Harris角点用来和第三张图匹配

二、Harris角点粗匹配
得到两张图的harris角点特征后，我们便可以开始粗匹配了。
粗匹配时其实还是采用的灰度相关加上双向最大原则的匹配方法，即取每个角点周围的一个邻域，然后求图一中角点和图二中角点的互相关函数。互为最大值则视为两点匹配成功。
这里要特别关心在边界周围的点，由于matlab中第一次拼接后的图像中含有nan（表示没有数据，但是运算时又会看成max）因此在角点检测的时候建议把邻域内含nan的点全部去掉。

三、基于randsac方法的Harris角点精确匹配
在精匹配前我们需要选定变换的模型，这里我们采用仿射变换，关于仿射变换，不多介绍。仿射变换可以用下面这个式子表示
P1=F*P2 F是仿射矩阵，P2是原始坐标，P1是变换后坐标
randsac方法也不多介绍，原理很简单，随机选取满足模型的最少匹配点（仿射模型时为3对），进行逆变换得到仿射矩阵F，然后将所有P2变换得到P1’，比较P1与P1’，记录误差小于阈值的点的个数，点数最多的那组点解出来的F便是最优模型。
得到的精匹配图如下
四、反变换更改图像坐标系
这里其实就是知道模型，进行仿射变换的过程。
注意，要将两张图调整到统一坐标系下，因为仿射变换的过程中可能会产生负的坐标值，因此在一开始处理的时候，我就将基准图，向上平移了，消除坐标中的幅值。
变换后图片如下，可以看出其实前两张图片主要是经过了平移变换

五、图像拼接
图像拼接的方式有很多，可以采用距离重叠中心的方式进行加权。
第一步，将两张图调整为一样的大小，注意调整大小时，要保持坐标不变。
第二步，计算每一行中重叠元素个数，
第三步，计算每一行中第一个重叠的位置
第四步，生成权值矩阵
第五步，加权求和。
拼接结果如图

六、附件：
原始图片三张

结果
程序代码

clc
clear
close all

%已知条件一，由于在统一位置拍摄，因此图像变换后面积变化不大
path1=imgetfile;
path2=imgetfile;
path3=imgetfile;
imrgb1=im2double(imread(path1));
imrgb1(size(imrgb1,1)+50,size(imrgb1,2))=nan;%上拉50个点，在仿射变换中，防止先到原图坐标系有负数
harrispipei(imrgb1,im2double(imread(path2)),5,1);%打包好的函数实现,不足的是点是自己找的
harrispipei(ans,im2double(imread(path3)),4,2);%打包好的函数实现,不足的是点是自己找的


function imo=harrispipei(imrgb1,imrgb2,r,cnt)
    [imrgb1,imrgb2]=im2samesizeup(imrgb1,imrgb2,nan);
    img1=(rgb2gray(imrgb1));
    img2=(rgb2gray(imrgb2));
    img1(find(img1==0))=nan;
    img2(find(img2==0))=nan;
    xyh1=HarrisPoint(img1,0.8,r,1,cnt-1);%  
    xyh2=HarrisPoint(img2,0.8,r,2,2*cnt-2);
    [xya1,xya2]=HarrisPointMosaic(img1,img2,xyh1,xyh2,r);%得到粗匹配点
    [A,xy1,xy2]=randsac([xya1,xya2],img1,img2)
    
%     if cnt==1
%         xy1=[147,158;139 193;216 296;187 207];
%         xy2=[253,163;246 199;320 296;289 212]-[224 0]+[0 50];
%     end
%     if cnt==2
%         xy1=[273 133;320 190;273 235];
%         xy2=[399 183;434 239;392 285]-[374 0];
%     end
%     if cnt==3
%         xy1=[469 245;356 422;370 490];
%         xy2=[683 220;576 397;589 471]-[500,0];
%     end
    A=rfangshe(xy2,xy1);
    [im2to1 dx1 dy1]=fangshe(img2,A);
    HarrisPoint(im2to1,0.8,r,3,0);
    plotcpl(img1,img2,xya1,xya2);
    plotcpl(img1,img2,xy1,xy2);
    %拼接
    imrgb21=fangshe(imrgb2,A);
    [imp1,imp2]=im2samesizeup(imrgb1,imrgb21,nan);%调整到同一个大小，这个是在图像上方补nan
    imo=pj12(imp1,imp2);%拼接在一起
    figure; imshow(imo);
end
%辅助函数
function plot_rxm_wh(img,pos,fig)%画红叉,pos为[w,h],支持多个图画到一个里面
    if fig==1
        figure
    end
    subplot(1,3,fig);
    imshow(img);hold on;
    for k=1:size(pos,1)
        plot(pos(k,2),pos(k,1),'r+');
    end
    hold off;
end
function [x y]=exchange(x,y)
    temp=x;x=y;y=temp;
end
function hwo=dropoutedge(hwi,mat2,r);%去掉靠近边界的点
    [hm,wm]=size(mat2);
    h=find(hwi(:,1)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,1)<hm-r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)<wm-r); 
    hwo=hwi(h,:);%删掉不满足的点
end
function [imo1,imo2]=im2samesizeup(im1,im2,num)%将两张图调整到一样大小
    [h1,w1,d1]=size(im1),[h2,w2,d2]=size(im2);
    hm=max(h1,h2);
    wm=max(w1,w2);
    if h1~=hm||w1~=wm
        temp=ones(hm,wm,d1)*num;
        im1=flipud(im1);
        temp(1:h1,1:w1,1:d1)=im1;
        im1=flipud(temp);
    end
    if h2~=hm||w2~=wm
        temp=ones(hm,wm,d2)*num;
        im2=flipud(im2);
        temp(1:h2,1:w2,1:d2)=im2;
        im2=flipud(temp);
    end
    imo1=im1;
    imo2=im2;
end
function [imo1,imo2]=im2samesizedown(im1,im2,num)%将两张图调整到一样大小
    [h1,w1,~]=size(im1),[h2,w2,~]=size(im2);
    hm=max(h1,h2);
    wm=max(w1,w2);
    if h1~=hm||w1~=wm
        im1(hm,wm,:)=num;
    end
    if h2~=hm||w2~=wm
        im2(hm,wm,:)=num;
    end
    imo1=im1;
    imo2=im2;
end
function imo=im2nan(im,dh)%将图片周围填充nan
    [h w d]=size(im);
    imo=ones(h+dh,w,d)*nan;
    imo(1:h,1:w,1:d)=im;
end
function plotcpl(im1,im2,xy1,xy2)%画控制点的连线
    figure;
    [im1,im2]=im2samesizeup(im1,im2,nan);
    [y,x]=size(im1);
    temp=[im1,im2];
    imshow(temp);
    hold on;
    xy2=xy2+[x,0];
    xx=[xy1(:,1),xy2(:,1)]';
    yy=[xy1(:,2),xy2(:,2)]';
    plot(xx,yy,'r');
    hold off;
end
% harris算法
function xy=HarrisPoint(img,sigma,r,fig,flag)%找到harris角点，输出为坐标[w,h]与xy对应
% img是输入的图像
% sigam是高斯滤波器参数
% T是阈值参数
    [X,Y]=size(img);
    gausFilter = fspecial('gaussian',r,sigma);%高斯滤波器
    
    fx=[-1,0,1]; Ix=filter2(fx,img);% 对图像求一阶梯度
    fy=[-1;0;1]; Iy=filter2(fy,img);
    Ix2=Ix.^2;Iy2=Iy.^2;Ixy=Ix.*Iy;
    
    gx=filter2(gausFilter,Ix2);
    gy=filter2(gausFilter,Iy2);
    gxy=filter2(gausFilter,Ixy);
    % 遍历所有像元，计算响应因子R 
    
    k = 0.04; 
    corrmat=gx.*gy-gxy.*gxy-k*(gx+gy).^2;%自相关矩阵
    corrmat(find(isnan(corrmat)==1))=0;
    sortedarray= sort(corrmat(:),'descend');%降序排序 
    corrmat(find(corrmat==0))=nan;
    count1 = size(find(sortedarray>0)); 
    T = mean(sortedarray(1:count1));
    
    [h w]=find(corrmat>T);%取较大值
    hw=deleteedge([h w],corrmat,r);%去掉靠近边界的点
    hw=deletenan(hw,img,r);%去掉邻域中含有nan的点
    hw=judgelinyumax(hw,corrmat,r);%找邻域最大
    xy=fliplr(hw);
    %显示用
    if fig~=0
        plot_rxm_wh(img,hw,fig);
    end
    if flag==1%删去图像左边的点,%基于先验知识的特征点排除
        xy=xy(find(xy(:,1)>X/2),:);
    end
    if flag==2%删去图像右边的点
        xy=xy(find(xy(:,1)<X/4),:);
    end
    %显示用
    if fig~=0
        plot_rxm_wh(img,fliplr(xy),fig);
    end
end
function [xyo1,xyo2]=HarrisPointMosaic(img1,img2,xy1,xy2,r)%对两幅图中的harris角点匹配,
    flag=0;%数据格式要求，前面那个是点少的
    hw1=fliplr(xy1);
    hw2=fliplr(xy2);
    nump1=size(hw1,1);%p1特征点数目
    nump2=size(hw2,1);%p2特征点数目
  
    if nump1>nump2%单纯为了减少计算量
        temp=nump1;nump1=nump2;nump2=temp;%交换
        temp=img1;img1=img2;img2=temp;%交换
        temp=hw1;hw1=hw2;hw2=temp;%交换
        flag=1;%表示前面那个是点少的
    end
    %计算误差矩阵
    P12=ones(nump1,nump2)*Inf;%保存差值
%     r=3;%周围邻域大小7,计算邻域灰度用 
    for k=1:nump1;
        ip1=img1(hw1(k,1)-r:hw1(k,1)+r,hw1(k,2)-r:hw1(k,2)+r);%用角点周围的一小个邻域
        for kk=1:nump2;
            ip2=img2(hw2(kk,1)-r:hw2(kk,1)+r,hw2(kk,2)-r:hw2(kk,2)+r);
            temp=(ip1-ip2).^2;P12(k,kk)=sum(temp(:));%求误差
        end
    end
    %找到匹配的点
    cnt=0;
    for k=1:nump1;%找最小值
        [maxw]=find(P12(k,:)==min(P12(k,:)));%找到 第h1行 最小值行所在位置
        if  P12(k,maxw)<mean(P12(:))
            if min(P12(:,maxw))==P12(k,maxw) %比较是不是 列 最小值
                cnt=cnt+1;%是 加一
                if flag==0
                    hwo1(cnt,:)=hw1(k,:);%输出点中对应起来
                    hwo2(cnt,:)=hw2(maxw,:);%输出点中对应起来
                else
                    hwo1(cnt,:)=hw2(maxw,:);%输出点中对应起来
                    hwo2(cnt,:)=hw1(k,:);%输出点中对应起来
                end
            end
        end
    end
    xyo1=fliplr(hwo1);
    xyo2=fliplr(hwo2);
end
function hw=deleteedge(hwi,mat2,r);%去掉角点靠近边界的点
    [hm,wm]=size(mat2);
    h=find(hwi(:,1)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,1)<hm-r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)<wm-r);    hw=hwi(h,:);%删掉不满足的点
end
function hw=deletenan(hw,img,r)%删除角点邻域含有nan的点
    for k=1:size(hw,1)
        linyu=img(hw(k,1)-r:hw(k,1)+r,hw(k,2)-r:hw(k,2)+r);
        if isempty(find(isnan(linyu)==1))==0%邻域存在nan
            hw(k)=0;
        end
    end
    hw=hw(find(hw(:,1)~=0),:);
end
function hwo=judgelinyumax(hwi,mat2,r)%输出邻域最大值的hw
    [hm,wm]=size(mat2);
    h=find(hwi(:,1)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,1)<hm-r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)>r);    hwi=hwi(h,:);%删掉不满足的点
    h=find(hwi(:,2)<wm-r);    hwi=hwi(h,:);%删掉不满足的点
    %舍弃边界周围的点 
    h=size(h,1);
    hwo(h,2)=0;%构造输出数据
    for h=1:h;
        linyu=mat2(hwi(h,1)-r:hwi(h,1)+r,hwi(h,2)-r:hwi(h,2)+r);
%         if isempty(find(isnan(linyu)==1))~=0
%             break;
%         end
        if max(linyu(:))==mat2(hwi(h,1),hwi(h,2));
            hwo(h)=hwi(h);
        end
    end
    hwo=hwi(find(hwo(:,1)~=0),:);
end
function [F xyo1 xyo2]=randsac(xy12,img1,img2)%采用随机方法得到最适合的模型,
%   img1是大图、基准图，img2是小图待匹配图,
%   输出A是从img2变到img1的变换矩阵
%   hwo1是img1中的控制点
%   hwo2是img2中的控制点
%         randp%随机产生的点
%         F%变换矩阵
%          temp=[0,0];变换后的坐标
    pointnum=size(xy12,1);%点的数目
    jilu(pointnum*30,pointnum+5)=0;%记录选的点与得到的均方误差
    ss=0;%误差累积
    x1=xy12(:,1);%转化到xy坐标
    y1=xy12(:,2);%转化到xy坐标
    x2=xy12(:,3);%转化到xy坐标
    y2=xy12(:,4);%转化到xy坐标
    xy1=[x1,y1];
    xy2=[x2,y2];
    k=100;
    cnt=1;
    while 1
        randp=randperm(pointnum,3);%产生3个不同随机数
        xy1t=xy1(randp,:);
        xy2t=xy2(randp,:);
        F=rfangshe(xy2t,xy1t);%%得到2变到1的变换矩阵,输入格式是xy坐标
        cdxs=F(1,1);
        cyds=F(1,2);
        saxc=F(2,1);
        syac=F(2,2);
        if det(F)>1.5 || det(F)<0.5%基于面积相同的改进策略
            continue;
        end
        if 1%基于小角度旋转的改进策略 
        end
        xy2t1=F*xy2xy1(xy2);%将img2点的坐标变到imga下,坐标变换公式
        %得到的temp是如下格式    %矩阵 [ x x x x ...
                                %       y y y y ...
                                %       1 1 1 1 ...]
        errors=xy2t1(1:2,:)-xy1';
        jilu(cnt,1:pointnum)=errors(1,:).^2+errors(2,:).^2; %得到误差矩阵
        jilu(cnt,pointnum+1:pointnum+3)=randp;              %得到误差矩阵
        jilu(cnt,pointnum+4)=sum(jilu(cnt,1:pointnum));     %临时存放总误差
        
        cnt= cnt+1;%自动选择迭代参数
        if cnt>=pointnum*k
            break;
        end
    end
%     es=min(min(jilu(:,1:pointnum)))*1.5;
    es=sum(sum(jilu(:,1:pointnum)))/pointnum/k/pointnum;
    for cnt=1:pointnum*30
        jilu(cnt,pointnum+5)=size(find(jilu(cnt,1:pointnum)<es),2);
    end
    pos=find(jilu(:,pointnum+5)==max(jilu(:,pointnum+5)));%满足的点最多
    xyo1=[xy12(jilu(pos(1),pointnum+1),1:2);xy12(jilu(pos(1),pointnum+2),1:2);xy12(jilu(pos(1),pointnum+3),1:2);];%hw1hw2(jilu(pos(1),pointnum+4),1:2);];
    xyo2=[xy12(jilu(pos(1),pointnum+1),3:4);xy12(jilu(pos(1),pointnum+2),3:4);xy12(jilu(pos(1),pointnum+3),3:4);];%hw1hw2(jilu(pos(1),pointnum+4),3:4);];
    F=rfangshe(xyo2,xyo1);%得到img2toimg1变换矩阵,hw to xy
    
    xy2t1=F*xy2xy1(xy2);%再次将2变到1下，找出所有合适点
    errors=xy2t1(1:2,:)'-xy1;
    errors=errors(:,1).^2+errors(:,2).^2;%得到误差矩阵
    sortedarray= sort(errors);%升序排序 
    num=4%返回点的个数
    xyo1=xy1(find(errors<sortedarray(num+1)),:);
    xyo2=xy2(find(errors<sortedarray(num+1)),:);%最终输出
    F=rfangshe(xyo2,xyo1);%%得到2变到1的变换矩阵,输入格式是xy坐标
end
%匹配算法
function [img1,img2]=imaxmin(img1,img2)%图像按大dao小排序
    [h1,~]=size(img1);
    [h2,~]=size(img2);
    if h2>h1%根据大小自动选择带匹配图像
        temp=img1;
        img1=img2;
        img2=temp;
    end
end    1 1 1 1
%    x y;]
    out=ones(3,size(in,1));
    out(1:2,:)=in';
end
function F=rfangshe(xy1,xy2)%求投影变换参数,输出为从xy1变到xy2的仿射矩阵 F=Y/X
%求ty变换的参数,x,y为对应控制点，输入为列向量
%     [x]=[a1 a2 a3][x]
%      y   b1 b2 b3  y
%      1   c1 c2 c3  1
    %按书上的公式来算
    %矩阵 [ w w w w ...
    %       h h h h ...
    %       1 1 1 1 ...]
    x=xy2xy1(xy1);
    y=xy2xy1(xy2);
    %得到修改了格式的输入量
    F=y*x'/(x*x');%求投影变换的参数
end
function [Ip dx dy]=fangshe(I,Fxy)%进行投影变换 Y=FX,dx,dy是进行平移的值,支持灰度和RGB图像
%投影变换函数 T是投影变换的模型参数，I是输入的图片
%输入：一张图片，投影变换矩阵[A B C]
%                           D E F
%                           G H I
%输出：一张图片
%映射到变换后图像 公式1
%     [x]=[a1 a2 a3][x]
%      y   b1 b2 b3  y
%      1   c1 c2 c3  1
%由变换后图像对应到原图像 公式2
% X =inv[B C]*[x-A]
% Y      E F   y-D
[h,w,d]=size(I);
cp = Fxy*[ 1 1 w w ; 
         1 h 1 h ;
         1 1 1 1];%求范围 执行公式1
Xpr = min( [ cp(1,:) 0 ] ) : max( [cp(1,:) w] ); % min x : max x,如果最下方点大于0，下方的点需要补0
Ypr = min( [ cp(2,:) 0 ] ) : max( [cp(2,:) h] ); % min y : max y


[Xp,Yp] = ndgrid(Xpr,Ypr);%产生对应的网状坐标
[wp,hp] = size(Xp); % = size(Yp)

X=Fxy\[Xp(:)';Yp(:)';ones([1,wp*hp])];%求对应的原始坐标，执行公式2
clear Ip;
xI = reshape( X(1,:),wp,hp)';
yI = reshape( X(2,:),wp,hp)';
Ip(:,:,1) = interp2(I(:,:,1), xI, yI, '*cubic'); % red
    if d==3;
        Ip(:,:,2) = interp2(I(:,:,2), xI, yI, '*cubic'); % green
        Ip(:,:,3) = interp2(I(:,:,3), xI, yI, '*cubic'); % blue
    end
    
dx=min([floor(min( [ cp(1,:) 0] )) 0]);%表示新的图像左下方点与原图像中的偏差, xn-xo
dy=min([floor(min( [ cp(2,:) 0] )) 0]);% ynew-yo
%     if dx<0%将图像调整到原图像坐标系
%         Ip=Ip(:,-dx:end);
%     end
%     if dy<0
%         Ip=Ip(1:(size(Ip,1)+dy),:);
%     end
end
function imo=pinjie(imall,hwi)%输入为数据立方体，输出拼接后的图像
%imall为数据立方体 hwi是立方体每一张图左上角对应的坐标
    [h,w,d]=size(imall);
    imall(find(isnan(imall)))=0;;
    pos1=find(imall(:,:,1)~=0);
    pos2=find(imall(:,:,2)~=0);
    im1=imall(:,:,1);
    im2=imall(:,:,2);
    im1(find(im1>0))=1;
    im2(find(im2>0))=1;
    im3=im1+im2;
    
    [h w]=find(im3==2);
    imall(h,w,2)=0;
    imo=imall(:,:,1)+imall(:,:,2);
end
function imo=pjgray(im1,im2);%输入的两张图片要求已经到统一坐标系，大小一致
    im1v=1-isnan(im1);%记录有数据的点
    im2v=1-isnan(im2);%记录有数据的点
    im12v=im1v+im2v;im12v(find(im12v==1))=0;im12v=im12v/2;%记录12都有数据的位置
    % w方向权重
    im12vnumw=sum(im12v,2);%对行求和,得到公共元素的个数
    im12vfirstw(size(im2,1),1)=0;%第一个出现的位置
    for k=1:size(im2,1)%
        temp=find(im12v(k,:),1,'first');
        if isempty(temp)==0
            im12vfirstw(k,1)=temp;%每行中第一个0出现的位置
        end
    end
    im1qw=im1v;%im1行加权系数 w方向上的
    im2qw=im2v;%im2行权系数
    for k=1:size(im2,1)
        if im12vnumw(k,1)~=0;
            im1qw(k,im12vfirstw(k,1):im12vfirstw(k,1)+im12vnumw(k,1)-1)=(im12vnumw(k,1):(-1):1)/im12vnumw(k,1);
            im2qw(k,im12vfirstw(k,1):im12vfirstw(k,1)+im12vnumw(k,1)-1)=(1:im12vnumw(k,1))/im12vnumw(k,1);
        end
    end
%     % h方向权重
%     im12vnumh=sum(im12v);%对列求和,得到公共元素的个数
%     im12vfirsth(size(im2,2),1)=0;%第一个出现的位置
%     for k=1:size(im2,2)%
%         temp=find(im12v(:,k),1,'first')
%         if isempty(temp)==0
%             im12vfirsth(k,1)=temp;%每行中第一个0出现的位置
%         end
%     end
%     im1qh=im1v;%im1列加权系数 h方向上的
%     im2qh=im2v;%im2列权系数
%     for k=1:size(im2,2)
%         if im12vnumh(1,k)~=0;
%             im1qh(im12vfirsth(k,1):im12vfirsth(k,1)+im12vfirsth(k,1)-1,k)=(im12vfirsth(k,1):(-1):1)/im12vfirsth(k,1)/2;
%             im2qh(im12vfirsth(k,1):im12vfirsth(k,1)+im12vfirsth(k,1)-1,k)=(1:im12vfirsth(k,1))/im12vfirsth(k,1)/2;
%         end
%     end
    im1(find(isnan(im1)==1))=0;
    im2(find(isnan(im2)==1))=0;
    imo=im1.*im1qw+im2.*im2qw;
    imo(find(imo==0))=nan;
end
function imo=pj12(im1,im2)%支持rgb和gray图像
    [h w d]=size(im1);
    if d==3
        imo(:,:,1)=pjgray(im1(:,:,1),im2(:,:,1));
        imo(:,:,2)=pjgray(im1(:,:,2),im2(:,:,2));
        imo(:,:,3)=pjgray(im1(:,:,3),im2(:,:,3));
    end
end
function hw=x2hw(x,base)%转成坐标值
    hw(:,1)=floor((x-1)./base)+1;
    hw(:,2)=mod(x-1,base)+1;
end

%%