www.gusucode.com > Matlab动力系统和时间序列分析工具箱 > Matlab动力系统和时间序列分析工具箱/lab432/toolbox/recurrence_plot.m
function varargout = recurrence_plot(varargin)
% RECURRENCE_PLOT M-file for recurrence_plot.fig
%
% cross-recurrence:
% RP=recurrence_plot(x,y,gui_mode,norm,treshold,embedding,delay);
% recurrence:
% RP=recurrence_plot(x,gui_mode,norm,treshold,embedding,delay);
% RP=recurrence_plot(x);
% Parameters:
% x,y - data vectors
% gui_mode - 'gui' or 'silent'
% norm - 1- Maximum norm
% 2- Euclidean norm
% 3- Minimum norm
% 4- Fixed amount of nearest neighbours
% 5- Distance coded matrix
% 6- Interdependent neighbours
%
% Example:
% RP=recurrence_plot(x,y,'gui',1,0.1,3,1);
%
% based on N. Marwan's CRP Toolbox
% Last Modified by GUIDE v2.5 19-Jul-2004 20:50:39
% last modified 10.01.05
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @recurrence_plot_OpeningFcn, ...
'gui_OutputFcn', @recurrence_plot_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
% End initialization code - DO NOT EDIT
% --- Executes just before recurrence_plot is made visible.
function recurrence_plot_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to recurrence_plot (see VARARGIN)
% Choose default command line output for recurrence_plot
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes recurrence_plot wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% RP=recurrence_plot(x,y,gui_mode,norm,treshold,embedding,delay);
global rpGLOBALS
set(handles.norm,'string',{'Maximum norm' 'Euclidean norm' 'Minimum norm' 'Fixed amount of nearest neighbours' 'Distance coded matrix' 'Interdependent neighbours'});
L=length(varargin);
rpGLOBALS.auto=0;
switch L
case 0
error('Not enougth input parameters');
case 1
rpGLOBALS.x=varargin{1};
rpGLOBALS.y=varargin{1};
set(hObject,'name','Recurrence plot');
rpGLOBALS.gui_mode='gui';
if length(rpGLOBALS.y)>1
set([handles.text_embedding, handles.frame_embedding handles.text_dimension...
handles.text_delay handles.edit_delay handles.edit_dimension],'enable','off');
set(handles.ax_graph,'visible','off');
else
axes(handles.ax_graph);
plot(rpGLOBALS.x);
end
case {2; 3; 4; 5; 6; 7}
rpGLOBALS.x=varargin{1};
if isa(varargin{2},'double')
rpGLOBALS.y=varargin{2};
set(hObject,'name','Cross recurrence plot');
rpGLOBALS.gui_mode='gui';
if length(rpGLOBALS.y)~=length(rpGLOBALS.x)
error('Time series must have equal dimensions');
end
axes(handles.ax_graph);
plot(rpGLOBALS.x);
if length(rpGLOBALS.y)>1
set([handles.text_embedding, handles.frame_embedding handles.text_dimension...
handles.text_delay handles.edit_delay handles.edit_dimension],'enable','off');
hold on
axes(handles.ax_graph);
plot(rpGLOBALS.y,'r');
% xlim([]);
hold off
end
if L>2
rpGLOBALS.gui_mode=varargin{3};
end
if L>3
set(handles.norm,'value',varargin{4});
end
if L>4
set(handles.edit_treshold,'string',num2str(varargin{5}));
rpGLOBALS.auto=1;
end
if L>5
set(handles.edit_dimension,'value',varargin{6});
rpGLOBALS.auto=1;
end
if L>6
set(handles.edit_delay,'string',num2str(varargin{7}));
rpGLOBALS.auto=1;
end
else
if strcmp(varargin{2},'gui')||strcmp(varargin{2},'silent')
rpGLOBALS.gui_mode=varargin{2};
else
warning('Unrecognized parameter ''gui_mode''. ''gui_mode'' set to ''gui''');
rpGLOBALS.gui_mode='gui';
end
rpGLOBALS.x=varargin{1};
rpGLOBALS.y=varargin{1};
set(hObject,'name','Recurrence plot');
if length(rpGLOBALS.y)>1
set([handles.text_embedding, handles.frame_embedding handles.text_dimension...
handles.text_delay handles.edit_delay handles.edit_dimension],'enable','off');
set(handles.ax_graph,'visible','off');
else
axes(handles.ax_graph);
plot(rpGLOBALS.x);
end
if L>2
set(handles.norm,'value',varargin{3});
end
if L>3
set(handles.edit_treshold,'string',num2str(varargin{4}));
rpGLOBALS.auto=1;
end
if L>4
set(handles.edit_dimension,'value',varargin{5});
rpGLOBALS.auto=1;
end
if L>5
set(handles.edit_delay,'string',num2str(varargin{6}));
rpGLOBALS.auto=1;
end
end
end
% --- Outputs from this function are returned to the command line.
function varargout = recurrence_plot_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
global rpGLOBALS
if rpGLOBALS.auto
compute_Callback(handles.compute, [], handles);
varargout{1} =rpGLOBALS.RP;
if strcmp(rpGLOBALS.gui_mode,'silent')
try
close(handles.figure1);
catch
end
end
else
uiwait;
try
varargout{1} =rpGLOBALS.RP;
catch
end
end
% try
% close(handles.figure1);
% catch
% end
function compute_Callback(hObject, eventdata, handles)
% hObject handle to compute (see GCBO)
% handles structure with handles and user data (see GUIDATA)
global rpGLOBALS
X=rpGLOBALS.x;
Y=rpGLOBALS.y;
E=str2num(get(handles.edit_treshold,'string'));
Norm=get(handles.norm,'value');
switch Norm
case 1
n='maxnorm';
case 2
n='euclidean';
case 3
n='minnorm';
case 4
n='fan';
case 6
n='inter';
case 5
n='distance';
end
M=get(handles.edit_dimension,'value');
T=str2num(get(handles.edit_delay,'string'));
RP=crp_big(X,Y,M,T,E,n,'sil');
% RP=crp_big(X,Y,1,1,E,n,'sil');
[s1 s2]=size(RP);
RP=rot90(double(RP));
if ~strcmp(n,'distance')
RP=-RP+max(max(RP));
end
RP=ceil(RP./max(max(RP)).*s1);
RP(find(RP==0))=1;
if strcmp(rpGLOBALS.gui_mode,'gui')
axes(handles.ax_plot);
image(RP);
colormap gray;
end
rpGLOBALS.RP=RP;
if ~rpGLOBALS.auto
uiresume;
end
% --- Executes during object creation, after setting all properties.
function norm_CreateFcn(hObject, eventdata, handles)
% hObject handle to norm (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: popupmenu 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
set(hObject,'string','Maximum norm|Euclidean norm|Minimum norm|Fixed amount of nearest neighbours|Interdependent neighbours|Distance coded matrix');
% --- Executes on selection change in norm.
function norm_Callback(hObject, eventdata, handles)
% hObject handle to norm (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: contents = get(hObject,'String') returns norm contents as cell array
% contents{get(hObject,'Value')} returns selected item from norm
if get(hObject,'Value')==5
set(handles.edit_treshold,'enable','off');
else
set(handles.edit_treshold,'enable','on');
end
% --- Executes during object creation, after setting all properties.
function edit_treshold_CreateFcn(hObject, eventdata, handles)
% hObject handle to edit_treshold (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: edit 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
function edit_treshold_Callback(hObject, eventdata, handles)
% hObject handle to edit_treshold (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of edit_treshold as text
% str2double(get(hObject,'String')) returns contents of edit_treshold as a double
function cancel_Callback(hObject, eventdata, handles)
% hObject handle to compute (see GCBO)
% handles structure with handles and user data (see GUIDATA)
global GSD_GLOBALS
GSD_GLOBALS.recurrence_plot=[];
close(gcf);
function edit_delay_CreateFcn(hObject, eventdata, handles)
% hObject handle to edit_delay (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: edit 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
function edit_delay_Callback(hObject, eventdata, handles)
% hObject handle to edit_delay (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of edit_delay as text
% str2double(get(hObject,'String')) returns contents of edit_delay as a double
% --- Executes during object creation, after setting all properties.
function edit_dimension_CreateFcn(hObject, eventdata, handles)
% hObject handle to edit_dimension (see GCBO)
% handles empty - handles not created until after all CreateFcns called
% Hint: edit 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
function edit_dimension_Callback(hObject, eventdata, handles)
% hObject handle to edit_dimension (see GCBO)
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of edit_dimension as text
% str2double(get(hObject,'String')) returns contents of edit_dimension as a double
V=str2double(get(hObject,'String'));
if V(get(hObject,'value'))>1
set(handles.edit_delay,'enable','on');
else
set(handles.edit_delay,'enable','off');
end
function xout=crp_big(varargin)
%CRP_BIG Creates a cross recurrence plot/ recurrence plot.
% CRP_BIG(X [,Y] [,param1,param2,...]) creates a cross
% recurrence plot/ recurrence plot. Results can be
% stored into the workspace. In contrast to CRP, long
% data series can be used. Results can be stored into
% the workspace.
%
% R=CRP_BIG(X,M,T,E) uses the dimension M, delay T
% and the size of neighbourhood E and creates a recurrence
% plot of X.
%
% R=CRP_BIG(X,Y,'distance','nonormalize') creates a
% distance coded matrix plot without normalization
% of the data.
%
% The source-data X and test-data Y can be one- or
% a two-coloumn vectors (then, in the first column
% have to be the time); if the test-data Y is not
% specified, a simple recurrence plot is created.
%
% Parameters: dimension M, delay T and the size of
% neighbourhood E are the first three numbers after
% the data series; further parameters can be used
% to switch between various methods of finding the
% neighbours of the phasespace trajectory, to suppress
% the normalization of the data and to suppress the
% GUI (useful in order to use this programme by other
% programmes).
%
% Methods of finding the neighbours.
% maxnorm - Maximum norm.
% euclidean - Euclidean norm.
% minnorm - Minimum norm.
% nrmnorm - Euclidean norm between normalized vectors
% (all vectors have the length one).
% fan - Fixed amount of nearest neighbours.
% inter - Interdependent neighbours.
% distance - Distance coded matrix (global CRP, Euclidean norm).
%
% Normalization of the data series.
% normalize - Normalization of the data.
% nonormalize - No normalization of the data.
%
% Parameters not needed to be specified.
%
% Examples: a=sqrt(100^2-(-71:71).^2); b=1:100;
% b(101:100+length(a))=-(a)+170;
% b(end+1:end+100)=100:-1:1;
% crp_big(b,1,1,.1,'max')
%
% See also CRP, CRP2 and CRQA.
% Copyright (c) 1998-2003 by AMRON
% Norbert Marwan, Potsdam University, Germany
% http://www.agnld.uni-potsdam.de
%
% Last Revision: 2004-01-23
% Version: 4.5.8
%
% This program is part of the new generation XXII series.
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or any later version.
%
% last modified 12.09.04 by Max Logunov
warning off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties
nogui=[];obj=[];mflag=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check and read the input
error(nargchk(1,9,nargin));
if isnumeric(varargin{1})==1 % read commandline input
varargin{9}=[];
i_double=find(cellfun('isclass',varargin,'double'));
i_char=find(cellfun('isclass',varargin,'char'));
% check the text input parameters for method, gui and normalization
check_meth={'ma','eu','mi','nr','fa','in','di'}; % maxnorm, euclidean, nrmnorm, fan, distance
check_norm={'non','nor'}; % nonormalize, normalize
check_gui={'gui','nog','sil'}; % gui, nogui, silent
temp_meth=0;
temp_norm=0;
temp_gui=0;
if ~isempty(i_char)
for i=1:length(i_char),
varargin{i_char(i)}(4)='0';
temp_meth=temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth');
temp_norm=temp_norm+strcmpi(varargin{i_char(i)}(1:3),check_norm');
temp_gui=temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui');
end
method=min(find(temp_meth));
nonorm=min(find(temp_norm))-1;
nogui=min(find(temp_gui))-1;
if isempty(method), method=1; end
if isempty(nonorm), nonorm=1; end
if isempty(nogui), nogui=0; end
if method>7, method=7; end
if nonorm>1, nonorm=1; end
if nogui>2, nogui=2; end
else
method=1; nonorm=1; nogui=0;
end
if nogui==0 & nargout>0, nogui=1; end
% get the parameters for creating RP
if max(size(varargin{1}))<=3
error('To less values in data X.')
end
x=double(varargin{1});
% if isempty(varargin{2}) | ~isnumeric(varargin{2}), y=x; end
% if ~isempty(varargin{2}) & isnumeric(varargin{2}), y=double(varargin{2}); end
if isempty(varargin{2}) | ~isnumeric(varargin{2}), y=x;
else, y=double(varargin{2}); end
if (isnumeric(varargin{2}) & max(size(varargin{2}))==1) | ~isnumeric(varargin{2})
y=x;
if ~isempty(varargin{i_double(2)}), m=varargin{i_double(2)}(1); else m=1; end
if ~isempty(varargin{i_double(3)}), t=varargin{i_double(3)}(1); else t=1; end
if ~isempty(varargin{i_double(4)}), e=varargin{i_double(4)}(1); else e=.1; end
else
if ~isempty(varargin{i_double(3)}), m=varargin{i_double(3)}(1); else m=1; end
if ~isempty(varargin{i_double(4)}), t=varargin{i_double(4)}(1); else t=1; end
if ~isempty(varargin{i_double(5)}), e=varargin{i_double(5)}(1); else e=.1; end
end
if e<0, e=1; disp('Warning: The threshold size E cannot be negative and is now set to 1.'), end
if t<1, t=1; disp('Warning: The delay T cannot be smaller than one and is now set to 1.'), end
t=round(t); m=round(m); mflag=method;
action='init';
Nx=length(x); Ny=length(y);
NX=Nx-t*(m-1);NY=Ny-t*(m-1);
if (NX<1 | NY<1) & nogui==0;
errordlg('The embedding vectors cannot be created. Dimension M and/ or delay T are to big. Please use smaller values.','Dimension/ delay to big')
waitforbuttonpress
end
if size(x,1)<size(x,2), x=x'; end
if size(y,1)<size(y,2), y=y'; end
if ~isempty(find(isnan(x)))
disp('Warning: NaN detected (in first variable) - will be cleared.')
for k=1:size(x,2), x(find(isnan(x(:,k))),:)=[]; end
end
if ~isempty(find(isnan(y)))
disp('Warning: NaN detected (in second variable) - will be cleared.')
for k=1:size(y,2), y(find(isnan(y(:,k))),:)=[]; end
end
if size(x,2)>=2
xscale=x(:,1);
if ~isempty(find(diff(xscale)<0)), error('First column of the first vector must be monotonically non-decreasing.'),end
if nonorm==1, x=(x(:,2)-mean(x(:,2)))/std(x(:,2)); else x=x(:,2); end
else
if nonorm==1, x=(x-mean(x))/std(x); end
xscale=(1:length(x))';
end
if size(y,2)>=2
yscale=y(:,1);
if ~isempty(find(diff(yscale)<0)), error('First column of the second vector must be monotonically non-decreasing.'),end
if nonorm==1, y=(y(:,2)-mean(y(:,2)))/std(y(:,2)); else y=y(:,2); end
else
if nonorm==1, y=(y-mean(y))/std(y); end
yscale=(1:length(y))';
end
ds=eye(m);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nogui
if nogui>0
hCRP=9999;
if nogui~=2
tx(1)={'Maximum norm'};
tx(2)={'Euclidean norm'};
tx(3)={'Minimum norm'};
tx(4)={'Euclidean norm of normalized distance'};
tx(5)={'fixed amount of nearest neighbours'};
tx(6)={'interdependent neighbours'};
tx(7)={'distance plot'};
disp(['use method: ', char(tx(method))]);
if nonorm==1, disp('normalize data'); else disp('do not normalize data'); end
end
action='compute';
if (NX<1 | NY<1)
disp('Warning: The embedding vectors cannot be created.')
disp('Dimension M and/ or delay T are to big. Please use smaller values.')
action='';
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch routines
switch(action)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% initialization
case 'init'
errcode=1;
xshuttle(:,1)=xscale;
yshuttle(:,1)=yscale;
xshuttle(:,2)=x;
yshuttle(:,2)=y;
ds=eye(m);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% computation
case 'compute'
errcode=11;
txt_cross='Cross ';
if size(x)==size(y)
if x==y, txt_cross=''; end
end
maxN=sqrt(m)*max(NX,NY);
if maxN>800
n=round(log(maxN/350)/log(2)); % orig war 150
maxN=maxN/(2^n);
maxN=round(maxN);NX0=NX;NY0=NY;smax=sqrt(m)*(max(abs(x))+max(abs(y)));
else
maxN=round(sqrt(m)*800);NX0=NX;NY0=NY;smax=sqrt(m)*(max(abs(x))+max(abs(y)));
end
s_norm=sqrt(m)*abs(max(x(:))-min(x(:)));
X=uint8(zeros(NY+(m-1)*t,NX+(m-1)*t));
% set(findobj('Tag','Status','Parent',findobj('Parent',hCRP,'Tag','CRPPlot')),'String','Build Embedding Vectors'),drawnow
if m*Nx>maxN
ax=ceil(m*Nx/maxN);
Nx2=ceil(Nx/ax);
x(Nx+1:Nx2*ax+t*(m-1))=x(Nx);
else
ax=1; Nx2=Nx;
end
x(Nx+1:Nx2*ax+t*(m-1))=x(Nx);
x0=x;
if m*Ny>maxN
ay=ceil(m*Ny/maxN);
Ny2=ceil(Ny/ay);
else
ay=1; Ny2=Ny;
end
y(Ny+1:Ny2*ay+t*(m-1))=y(Ny);
y0=y;
% set(findobj('Tag','Status','Parent',findobj('Parent',hCRP,'Tag','CRPPlot')),'String','Compute CRP'),drawnow
k=0;
for i=1:ax,
x=x0(1+Nx2*(i-1):Nx2+Nx2*(i-1)+t*(m-1));
for j=1:ay,
y=y0(1+Ny2*(j-1):Ny2+Ny2*(j-1)+t*(m-1));
k=k+1;
Nx=length(x); Ny=length(y);
NX=Nx-t*(m-1);NY=Ny-t*(m-1);
%%
ii=(1:NX)';jj=0:t:0+(m-1)*t;
ii=reshape(ii(:,ones(1,m))+jj(ones(NX,1),:),m*NX,1);
x1=x(ii);
x2=reshape(x1,NX,m);
ii=(1:NY)';jj=0:t:0+(m-1)*t;
ii=reshape(ii(:,ones(1,m))+jj(ones(NY,1),:),m*NY,1);
y1=y(ii);
y2=reshape(y1,NY,m);
y2=y2*ds; % switch vectors
[NX, mx] = size(x2);
[NY, my] = size(y2);
clear jx jy
switch(mflag)
%%%%%%%%%%%%%%%%% local CRP, fixed distance
case {1,2,3}
errcode=111;
px = permute(x2, [ 1 3 2 ]);
py = permute(y2, [ 3 1 2 ]);
s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
switch mflag
case 1
%%%%%%%%%%%%%%%%% maximum norm
s = max(abs(s1),[],3);
matext=[num2str(round(100*e)/100) '\sigma (fixed distance maximum norm)'];
case 2
%%%%%%%%%%%%%%%%% euclidean norm
s = sqrt(sum(s1.^2, 3));
matext=[num2str(round(100*e)/100) '\sigma (fixed distance euclidean norm)'];
case 3
%%%%%%%%%%%%%%%%% minimum norm
s = sum(abs(s1), 3);
matext=[num2str(round(100*e)/100) '\sigma (fixed distance minimum norm)'];
end
X1=255*s/max(s(:))<(255*e/max(s(:)));
X0=(uint8(X1))'; clear s s1 x1 y1 px py X1
X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
X(NY0+1:end,:)=[];
X(:,NX0+1:end)=[];
NX=NX0; NY=NY0;
% matext=[num2str(round(100*e)/100) '\sigma (fixed distance maximum norm)'];
%%%%%%%%%%%%%%%%% local CRP, normalized distance euclidean norm
case 4
% set(findobj('Tag','Status','Parent',findobj('Parent',hCRP,'Tag','CRPPlot')),'String',['Normalize Embedding Vectors ',num2str(ax*(i-1)+j)]),drawnow
Dx=sqrt(sum(((x2(:,:)).^2)'))';
Dy=sqrt(sum(((y2(:,:)).^2)'))';
x1=x2./repmat(Dx,1,m);x2=x1;
y1=y2./repmat(Dy,1,m);y2=y1; clear Dx Dy y1 x1
px = permute(x2, [ 1 3 2 ]);
py = permute(y2, [ 3 1 2 ]);
s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
s = sqrt(sum(s1.^2, 3));
X0=(uint8(255*s/max(s(:)))<(255*e/max(s(:))))'; clear s s1 x1 y1 px py
X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
X(NY0+1:end,:)=[];
X(:,NX0+1:end)=[];
NX=NX0; NY=NY0;
matext=[num2str(round(100*e)/100) '\sigma (normalized distance euclidean norm)'];
%%%%%%%%%%%%%%%%% local CRP, fixed neigbours amount
case 5
if e>=1
e=round(e)/100;
txt=['The value for fixed neigbours amount has to be smaller '...
'than one. Continue the computation with a value of ' ...
num2str(e)];
if nogui==0
warndlg(txt,'Threshold value mismatch');
drawnow
waitforbuttonpress
set(findobj('Tag','Size','Parent',gcf),'String',num2str(e))
end
end
px = permute(x2, [ 1 3 2 ]);
py = permute(y2, [ 3 1 2 ]);
s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
s = sqrt(sum(s1.^2, 3));
mine=round(NY*e);
[SS, JJ]=sort(s'); JJ=JJ';
X1(NX*NY)=uint8(0); X1(JJ(:,1:mine)+repmat([0:NY:NX*NY-1]',1,mine))=uint8(1);
X0=reshape(X1,NY,NX); clear X1 SS JJ s px py
X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
X(NY0+1:end,:)=[];
X(:,NX0+1:end)=[];
matext=[num2str(round(1000*e)/10) '% (fixed neighbours amount)'];
NX=NX0; NY=NY0;
%%%%%%%%%%%%%%%%% local CRP, interdependent neigbours
case 6
if e>=1
e=round(e)/100;
txt=['The value for fixed neigbours amount has to be smaller '...
'than one. Continue the computation with a value of ' ...
num2str(e)];
end
px = permute(x2, [ 1 3 2 ]);
py = permute(y2, [ 1 3 2 ]);
px2 = permute(x2, [ 3 1 2 ]);
py2 = permute(y2, [ 3 1 2 ]);
s1 = px(:,ones(1,NX),:) - px2(ones(1,NX),:,:);
sx = sqrt(sum(s1.^2, 3));
s1 = py(:,ones(1,NY),:) - py2(ones(1,NY),:,:);
sy = sqrt(sum(s1.^2, 3));
mine=round(min(NX,NY)*e);
[SSx, JJx]=sort(sx);%SSx(1,:)=[]; JJx(1,:)=[];
[SSy, JJy]=sort(sy);%SSy(1,:)=[]; JJy(1,:)=[];
ee=mean(SSy(mine:mine+1,:));
X0=uint8(zeros(NY,NX));
for ii=1:min(NX,NY),
JJx((JJx(1:mine,ii)>min(NX,NY)),ii)=1;
X0(ii,JJx(1:mine,ii))=sy(ii,JJx(1:mine,ii))<=ee(ii);
end
X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
X(NY0+1:end,:)=[];
X(:,NX0+1:end)=[];
NX=NX0; NY=NY0;
clear X1 SS* JJ* s sx sy s1 px py
matext=[num2str(round(1000*e)/10) '% (interdependent neighbours)'];
%%%%%%%%%%%%%%%%% global CRP
case 7
px = permute(x2, [ 1 3 2 ]);
py = permute(y2, [ 3 1 2 ]);
s1 = px(:,ones(1,NY),:) - py(ones(1,NX),:,:);
s = sqrt(sum(s1.^2, 3))';
X0=uint8(255*s/s_norm);
X(1+Ny2*(j-1):Ny2+Ny2*(j-1),1+Nx2*(i-1):Nx2+Nx2*(i-1))=X0;
X(NY0+1:end,:)=[];
X(:,NX0+1:end)=[];
NX=NX0; NY=NY0;
matext='';
end
end, end
%%%%%%%%%%%%%%%%% show CRP
xout=X;
end
warning on