前言

程序流程图如下：

主程序

• 这是一个MatlabGUI界面，
• 包含了读入参考图像，浮动图像，
• 显示配准前后参考图像与浮动图像的融合效果图
• 调用GLPF.m完成高斯低通预处理
• 调用函数Powll.m，实现图像的配准
• 显示配准结果

function varargout = ImageRegistration(varargin)
%%%%%%% %%% %%% %%%%%%% %%%%% %%%%
gui_Singleton= 1;
gui_State = struct( 'gui_Name',      mfilename, ...
                    'gui_Singleton', gui_Singleton, ...
                    'gui_OpeningFcn',@ImageRegistration_OpeningFcn, ...
                   'gui_OutputFcn',  @ImageRegistration_OutputFcn, ...
                    'gui_LayoutFcn', [],...
                    'gui_Callback',  []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
[varargout{1:nargout}]=gui_mainfcn(gui_State,varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
addpath(pwd);
function ImageRegistration_OpeningFcn(hObject,eventdata,handles,varargin)
handles.output = hObject;
% handles.OAname='PSO';
% handles.MIname='MI';
% handles.NumOfVar=0;
% set(handles.edit1,'visible','off');
% set(handles.edit2,'visible','off');
% set(handles.text3,'visible','off');
% set(handles.text4,'visible','off');
guidata(hObject, handles);

function varargout = ImageRegistration_OutputFcn(hObject,eventdata, handles)
varargout{1} = handles.output;
function pushbutton1_Callback(hObject, eventdata, handles)
% clc
%%%%%%%%调用OpenImage.m读入参考图像并获取
%%%%文件名、图像大小
Image_I.FileInformation.IsImage=0;
while Image_I.FileInformation.IsImage==0
Image_I=OpenImage;
end
delete(Image_I.figurel);
handles.Imsizel=Image_I.FileInformation.imsize;
handles.filenameI=Image_I.FileInformation.filename;
handles.names_dispI=Image_I.FileInformation.names_disp;
set(handles.text2,'String' ,handles.names_dispI);
guidata(hObject, handles);
%%%%%%%% 显示参考图像
axes(handles.axes1 )
I=imread(handles.filenameI);
imshow(I)

function pushbutton2_Callback(hObject, eventdata, handles)
clc
%%%%%%%%调用OpenImage.m读入浮动图像并获取
% 文件名、图像大小
Image_J.FileInformation.IsImage=0;
while Image_J.FileInformation.IsImage==0
Image_J=OpenImage;
end
delete(Image_J.figure1);
handles.ImsizeJ=Image_J.FileInformation.imsize;
handles.filenameJ=Image_J.FileInformation.filename;
handles.names_dispJ=Image_J.FileInformation.names_disp;
set(handles.text3,'String' ,handles.names_dispJ);
guidata(hObject, handles);
%%%%%%%%显示浮动图像%%%%%%%%%%
axes(handles.axes2)
J=imread(handles.filenameJ);
imshow(J)

function pushbutton3_Callback(hObject, eventdata, handles)
clc;
%%%%%%%检测是否已输入参考图像与浮动图像%%
axesIbox=get(handles.axes1,'box');
axesJbox=get(handles.axes2,'box');
if strcmp(axeslbox,off) | strcmp(axesJbox,'off')
errordlg('Please select Image for Registration',Error')
error('NO Image !')
end

%%%%%%检测参考图像与浮动图像大小是否相同%%
handles.isSameSizelJ=strcmp(handles.ImsizeI,handles.ImsizeJ);
if handles.isSameSizeIJ~= 1
errordlg('Please Select the Same Size Image',Error')
error('Image Size doesn"t match')
end 
%%%%%%%%读入并复制图像，一幅用于配准过程，
% 另一-幅用于配准后输出
I1=imread(handles.filenameD);
J1=imread(handles.filenameJ);
handles.Old_I=I1;
handles.Old_J=J1;
I=ImageTransfer_add(I1);
J=ImageTransfer_add(J1);
handles.I=I;
handles.J=J;
guidata(hObject,handles);
%%%%%%%%显示配准前参考图像与浮动图像的
% “融合”效果图%%%%%%%
axes(handles.axes3)
imshow(uint8(I+J))
%%%%%%%%调用函数GLPF.m,完成高斯低通预
%处理%%%%%%%%%%%%%%%
[I,J]=GLPF(handles);
handles.I=I;
handles.J=J;
guidata(hObject, handles);
%%%%%%调用函数Powell.m,实现图像配准%%%%%
tic
RegistrationParameters=Powell(handles);
toc
ElapsedTime=toc;

%%%%%%%%显示配准结果%%%%%%%%%%%%%
handles.RegistrationParameters=RegistrationParameters;
x=RegistrationParameters(1);
y=RegistrationParameters(2);
ang=RegistrationParameters(3);
MI_Value=RegistrationParameters(4);
RegistrationResult=sprintf('X,Y,Angle=[%.5f][%.5f][%.5f]',x,y,ang)
MI_Value=sprintf('MI_Value=[%.4f]',MI_Value)
ElapsedTime=sprintf('Elapsed Time=[%.3f]',ElapsedTime)
axes(handles.axes4)
[FusionImage,RegistrationImage]=Fusion(handles);
imshow(FusionImage)
axes(handles.axes5)
imshow(RegistrationImage)
set(handles.text8,'string',RegistrationResult);
set(handles.text9,'string' ,MI_Value);
set(handles.text10,'string' ,ElapsedTime);

Powell（鲍威尔优化法）

鲍威尔法又称方向加速法，它由Powell于1964年提出，是利用共轭方向可以加快收敛速度的性质形成的一种搜索方法。该方法不需要对目标函数进行求导，当目标函数的导数不连续的时候也能应用，因此，鲍威尔算法是一种十分有效的直接搜索法。改进后的Powell算法步骤为： (转载于：原文链接)

• Step1：选取初始数据。选取初始点x0,n个线性无关的初始搜索方向d0,d1,...,dn-1;给定允许误差Err>0,令k=0;

• Step2：进行基本搜索。令y0=x(k),依次沿d0,d1,...,dn-1进行一维搜索。对一切j=1,2,...,n，记f(y(j-1)+λ(j-1)d(j-1)) = minf(y(j-1)+λd(j-1))， y(j) = f(y(j-1)+λ(j-1)d(j-1)) ;

• Step3：检查是否满足终止条件。取加速度方向d(n)=y(n)-y(0);若||d(n)||<Err,迭代终止，得yn为问题的最优解；否则转Step4。

• Step4：确定搜索方向，按照上式公式确定m，若验证式子成立，转Step5；否则，转Step6。

• Step5：调整搜索方向。从点yn出发沿方向dn进行一维搜索，求出λn,使得f(yn+λn$\times$dn)=minf(yn+λ$\times$dn);令x(k+1)=yn+λn$\times$dn。再令d(j)=d(j+1)，j=m,m+1,...,n-1.k=k+1,返回Step2。

• Step6：不调整搜索方向，令x(k+1)=yn; k=k+1,返回Step2。

  4.3 典型例题
  用修正的Powell算法求f(x)=x1*x1+2*x2*x2-4*x1-2*x1*x2的最优解。X0=[1,1]';收敛精度Err=0.001。
  解：X0(1)=X0=[1,1]'；f1=f(x0(1))=-3;
  第一轮进行循环：
  沿坐标轴方向向S1=[1,0]方向进行一维搜索:
  X1(1) = X0(1)+a(1)*S1(1) = [1,1]'+a1(1)[1,0]'=[1+a1(1),1]'
  将X1(1)带入上式，求导可以确定a1(1)=2时，得到极小值。此时，X1(1)=[3,1]; f(X1(1))=-7。
  沿坐标轴方向向S2=[0,1]方向进行一维搜索:
  X2(1) = X1(1)+a2(1)*S2(1) = [3,1]'+[0,a2(1)]'=[3,1+a2(1)]'。
  将X2(1)代入上式，求导可以确定a2(1)=1/2,得到极小值。此时，X2(1)=[3,1.5]; f(X2(1))=-7.5。
  检验是否满足终止迭代条件||X2(1)-X0(1)|| = 2.06>Err.
  计算各个方向上的函数下降量：
  △1=f(X0(1))-f(X1(1))=4;
  △2=f(X1(1))-f(X2(1))=0.5;
  △m=min{4,0.5}=4;
  映射点：
  X(1) = 2X2(1)-X0(1)=2*[3,1.5]'-[1,1]=[5,2];
  条件验证：f3=f(X(1))=-7, f3<f1
  (f1-2f2+f3)(f1-f2-△m)^2 = 1.25 < △m/2*(f1-f3)^2 =32 满足。
  因此，得到新的搜索方向：S(1)=X2(1)-X0(1)=[3,1.5]-[1,1]=[2,1.5]';
  
  沿S(1)方向再次做一维搜索：
  X3(1)=X2(1)+a3(1)*S(1)=[3,1.5]+a3(1)[2,1.5] = [3+2*a3(1),1.5+1.5a3(1) ];
       当a3(1)=0.4时，取得极小值。此时，X3(1)=[19/5,17/10]; f(X3(1)) = -7.9。
  
  以X0(2) = X3(1)作为新的起点，沿(e2,S(1))方向进行一维搜索，开始进入第二个循环：
  X2(0)=[19/5,17/10]; f(X2(0))=-7.9;
  沿e2=[0,1]方向进行一维搜索：
  X1(2) = [19/5,19/10]; f(X1(2))=-7.98 
  沿S(1)=[2,1.5]方向进行一维搜索:
  X2(2) = [99/25,97/50];  f(X2(2)) =-7.996
  检验：||X2(2)-X2(0)|| = 0.288>Err
  继续进行迭代：△1=0.08;
   △2=0.016;
   △m=0.08;
  映射点：X(2)=2*X2(2)-X0(2)=[103/25,109/50]; f(X(2))=-7.964>f1 条件不满足。
  新的方向：S(2)=[99/25,97/50]-[19/5,17/10]=[4/25,12/50];
  
  进行继续迭代时，取X3(0)=X2(2)，沿着(e2,S(1))方向一维搜索进行第三轮循环：
  X3(1)=[99/25,99/50] ; f1=-7.9992;
  X3(2)=[3.9992,1.988]; f2=-7.99984;
  检验：||X2(3)-X0(3)||=0.0577 >Err;
  X(3)=[4.024,2.036]'; f3=-7.99856;  f3>f1.
  由于S(1)即S(2)为共轭方向，目标函数是二次函数，若沿S(2)方向进行一维搜索得到
  X(2)=[4,2]，F(X(2))=-8，此为目标函数的最优解。
  

OpenImage

查看图片

function varargout = testListBoxlj(varargin)
% TESTLISTBOXLJ M-file for testListBoxlj.fig
%      TESTLISTBOXLJ, by itself, creates a new TESTLISTBOXLJ or raises the existing
%      singleton*.
%
%      H = TESTLISTBOXLJ returns the handle to a new TESTLISTBOXLJ or the handle to
%      the existing singleton*.
%
%      TESTLISTBOXLJ('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in TESTLISTBOXLJ.M with the given input arguments.
%
%      TESTLISTBOXLJ('Property','Value',...) creates a new TESTLISTBOXLJ or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before testListBoxlj_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to testListBoxlj_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help testListBoxlj

% Last Modified by GUIDE v2.5 19-Sep-2005 10:36:50

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @testListBoxlj_OpeningFcn, ...
                   'gui_OutputFcn',  @testListBoxlj_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before testListBoxlj is made visible.
function testListBoxlj_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to testListBoxlj (see VARARGIN)

% Choose default command line output for testListBoxlj
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

initial_dir = pwd;
ss = load_listbox(initial_dir,handles);
handles.ss = ss;
% UIWAIT makes testListBoxlj wait for user response (see UIRESUME)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %     add by lijing
% handles.parhandles=varargin{1};

guidata(handles.figure1, handles);

% while handles.ss.IsImage==0
uiwait(handles.figure1); %%%%%%%只有在这里加才不会出错
% handles.figure1
% handles.ss.IsImage
% end

% --- Outputs from this function are returned to the command line.
function varargout = testListBoxlj_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

% varargout{1} = handles.output;
varargout{1} = handles;
% delete(handles.figure1)%%%%%%%%%%

% --- Executes during object creation, after setting all properties.
function listbox1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

% --- Executes on selection change in listbox1.
function listbox1_Callback(hObject, eventdata, handles)
% hObject    handle to listbox1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% get(handles.figure1,'SelectionType');
clc
index_selected = get(handles.listbox1,'Value');
% file_list = get(handles.listbox1,'String');	
%%%%%%%%%%%%大问题％％％％％％％％％％％％％％％％％％％％％％％
filename = handles.ss.sorted_names{index_selected}; 
handles.ss.filename=[pwd,'\',filename];
handles.ss.names_disp=handles.ss.sorted_names_disp{index_selected};
handles.ss.imsize=handles.ss.imsize_disp{index_selected}; % （返回）
% handles.ss.f=handles.figure1


    if index_selected<=handles.ss.cnPiont    
        
        cd(handles.ss.filename)
        ss=load_listbox(pwd,handles);
        handles.ss=ss;
        handles.ss.IsImage=0; % 返回一个‘文件’或‘图片’的记录
% filename = handles.ss.sorted_names{index_selected}; 
% handles.ss.filename=[pwd,'\',filename];
% handles.ss.names_disp=handles.ss.sorted_names_disp{index_selected};
% handles.ss.imsize=handles.ss.imsize_disp{index_selected}; % （返回）
    else
        handles.ss.IsImage=1;  % （返回）
        
end
%%%%%%%%%%%%
guidata(handles.figure1, handles);
uiresume(handles.figure1);%%%%%%%只有在这里加才不会出错

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ss = load_listbox(dir_path,handles)
cd (dir_path)
dir_struct = dir(dir_path);

[sorted_names1,sorted_index1] = sortrows({dir_struct.name}');

k = max(sorted_index1);
if k==2
    sorted_names{1}='.';
    sorted_names{2}='..';
    sorted_names_disp{1}='.';
    sorted_names_disp{2}='..';
    imsize_disp{1}='null';
    imsize_disp{2}='null';
    cnPiont=2;
else
    sorted_names{1}='.';
    sorted_names{2}='..';
    sorted_names_disp{1}='.';
    sorted_names_disp{2}='..';
    imsize_disp{1}='null';
    imsize_disp{2}='null';
    cn=2;
%%%%%% 找文件夹
for i=1:k
%     [path,name,ext,ver] = fileparts(sorted_names1{i});
      [path,name,ext] = fileparts(sorted_names1{i});
     switch ext
         case''
            cn=cn+1;
            cnarray(cn)=cn;
            sorted_names{cn}=name;
            sorted_names_disp{cn}=name;
            %name
            imsize_disp{cn}='null';
    end
end 
cnPiont=cn;
% guidata(handles.figure1,handles)
%%%%%% 找图片
for i=1:k
    [path,name,ext] = fileparts(sorted_names1{i});
     switch ext
        case {'.bmp','.jpg','.jpeg','.tif','.png'}
            cn=cn+1;
            cnarray(cn)=cn;
            name=[name,ext];
            sorted_names{cn}=name;
            
            temp=imread(name);
            [m,n]=size(temp);
            m=num2str(m);
            n=num2str(n);
            imsize=['   ',m,'*',n,'(size)'];  %%%  要返回
            name_disp=[name,imsize];
            sorted_names_disp{cn}=name_disp;
            imsize_disp{cn}=imsize;
    end
end 
end

ss.sorted_names=sorted_names; %%%  要返回
ss.sorted_names_disp=sorted_names_disp; %%%  要返回
ss.imsize_disp=imsize_disp;
ss.cnPiont=cnPiont; %%%  要返回

% % handles.file_names = sorted_names;
% handles.file_names =sorted_names_disp;
% handles.is_dir = [dir_struct.isdir];
% %handles.sorted_index = [sorted_index];
% handles.sorted_index = [cnarray];

set(handles.listbox1,'String',sorted_names_disp,...
    'Value',1)
set(handles.text1,'String',pwd)  

Fusion

融合图像

function [FusionImage,RegistrationImage]=Fusion(handles)
%%%%%%获取原图像及配准参数%%%%%%%%%%
I=handes.Old_I;
J=handles.Old_ J;
x=handles.RegistrationParameters(1);
y=handles.RegistrationParameters(2);
ang=-handles.RegistrationParameters(3);
%%%%%%%%对浮动图像空间变换及插值，以便得到
% 配准后的输出图像%%%%%%%
[nrows,ncols]=size(J);
width = nrows;
height = ncols;
new_J = uint8(zeros(width,heigh));
a=(width-1)/2;
c-a;
b=(height-1)/2;
d=b;
rad=pi/180* ang;
t1=[1 0 0;0 1 0;x y 1];
t2=[1 0 0;0 1 0;-a-b 1];
t3=[cos(rad)-sin(rad) O;sin(rad) cos(rad) 0;0 0 1];
t4=[1 0 0;0 1 0;c d 1];
T=t2*t3*t4*tl;
tform=maketform('affine',T);
tx=zeros(width,height);
ty=zeros(width,height);
for i=1:width
for j=1:height;
tx(ij)=i;
end
end
for i= l:width
for j=l:height;
ty(ij)=j;
end
end
[w z]=forminv(tform,tx,ty); 
for i= 1:width
for j= 1:height
source_x= w(ij);
source_y = z(ij);
if (source_x>= width-1 || source_y >=height-1 || double(uint16(source_x)) <=0 || double(uint16(source_y))<=0)
new_J(i,j)=J(1,1);
else
if (source_x/double(uint16(source_x))== 1.0) & (source_y/double(uint16(source_y))= 1.0)
new_J(i,j)= J(int16(source_x),int16(source_y));
else
%a = double(round(source_ x));
%b = double(round(source_ y));
%new_ J(des_ x,des_ y)= J(a,b);
a = double(uint16(source_x));
b = double(uint16(source_y));
x11 = double( J(a,b));
x12 = double( J(a,b+1));
x21 = double( J(a+1,b));
x22 = double( J(a+1,b+1));
new_J(i.j) = uint8( (b+l-source_ y)*((source_ x-a)*x21 + (a+l-source_ x)*x11) +(source_ y-b) * ((source_ x-a)*x22 +(a+l-source_ _x)* x12) );
end
end 
end
end
J=new_J;
I=uint8(I);
J=uint8(J);
RegistrationImage=uint8(J);
%%%%%%产生“融合”后的效果图%%%%%%%%%
I=double(I)/255;
J=double(J)/255;
J = double(zeros(width,height));
for m= l:width
for n=1 :height
if I(m,n)>0.999 || J(m,n)>0.999
IJ(m,n)=0.8;
elseif I(m,n)==0 || J(m,n)==0
IJ(m,n)=0.01;
else
IJ(m,n)=(I(m,n)*0.3+J(m,n)*0.7);
end
end
end
IJ=IJ*255;
IJ=uint8(IJ);
FusionImage=IJ;

PV

几何变换及互信息计算

function [mi]=PV(x,y,ang,handles)
a=handles.I;
a=double(a);
b=handles.J;
b=double(b);
[M,N] = size(a);
hab = zeros(256,256);
ha = zeros(1,256);
hb = zeros(1,256);
if max(max(a))~=min(min(a))
a = (a-min(min(a)))/(max(max(a))-min(min(a)));
else
a = zeros(M,N);
end
if max(max(b))-min(min(b))
b = (b-min(min(b)))/(max(max(b))-min(min(b)));
else
b = zeros(M,N);
end
a = double(int16(a*255))+1;
b = double(int16(b*255))+1;
[width,height]=size(b);
u=(width-1)/2;
v=(height-1)/2;
rad=pi/180*ang;
t1=[1 0 0;0 1 0;x y 1];
t2=[1 0 0;0 1 0;-u -v 1];
t3=[cos(rad) -sin(rad) 0;sin(rad) cos(rad) 0;00 1];
t4=[1 0 0;0 1 0;u v 1];
T=t2*t3*t4*t1;
tform=maketform( affine',T);
coordinate_x=zeros(width,height);
coordinate_y=zeros(width,height);
for i=1:width
for j=1:height;
coordinate_x(i,j)=i;
end
end

for i=1:width
for j=1:height;
coordinate_x(i,j)=j;
end
end
[w z]=tforminv(tform,coordinate_x,coordinate_y);
for i= 1:width
for j= 1:height
source_X = w(ij);
source_y = z(ij);
if (source_X > width-1 || source_y >height-1 ||double(uint16(source_x))<=1 ||double(uint16(source_y))<=1)
hab(a(1,1),a(1,1))=hab(a(1,1),a(1,1))+1;
else
m=fix(source_x);
n=fix(source_y);
index_b =b(i,j);
index_a0=a(m,n);
index_a1=a(m+1,n);
index_a2- a(m,n+1);
index_a3=a(m+1,n+1); 
dx=source_x-m;
dy=source_y-n;
hab(index_a0,index_b)= hab(index_a0,index_b)+(1-dx)*(1-dy);
hab(index_al,index_b)= hab(index_al,index_b)+dx*(1-dy);
hab(index_a2,index_b)= hab(index_a2,index_b)+(1-dx)*dy;
hab(index_a3,index_b)= hab(index_a3,index_b)+dx*dy;
end
end
end
habsum = sum(sum(hab));
index = find(hab~=0);
pab = hab/habsum;
Hab = sum(sum(-pab(index).*log2(pab(index))));
pa=sum(pab');
index=find(pa~=0);
Ha= sum(sum(pa(index).*log2(pa(index))));
pb=sum(pab);
index=find(pb~=0);
Hb= sum(sum( pb(index).*log2(pb(index))));
mi = Ha+Hb-Hab;

GLPF

低通滤波

function [I,J]=GLPF(handles)
I=handles.I;
J=handles.J;
I=double(I);
J=double(J);
PQ=2*size(I); 
M=PQ(1);
N=PQ(2);
u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D0=0.05*PQ(2);
H=exp(-(U.^2+V.^2)/(2*(D0^2)));
FI=fft2(I,size(H,1),size(H,2));
FJ=fft2(J,size(H,1),size(H,2));
gi=real(ifft2(H.*FD)); 
gj=real(ifft2(H.*FJ));
gi=gi(1:size(I,1),1:size(I,2));
gj=gj(1:size(L,1),1:size(I,2));
I=uint8(gi);
J=uint8(gj);