diff --git a/src/Resources/DigImCorr_resources/icon_16.png b/src/Resources/DigImCorr_resources/icon_16.png new file mode 100644 index 0000000..82511d4 Binary files /dev/null and b/src/Resources/DigImCorr_resources/icon_16.png differ diff --git a/src/Resources/DigImCorr_resources/icon_24.png b/src/Resources/DigImCorr_resources/icon_24.png new file mode 100644 index 0000000..c47babc Binary files /dev/null and b/src/Resources/DigImCorr_resources/icon_24.png differ diff --git a/src/Resources/DigImCorr_resources/icon_48.png b/src/Resources/DigImCorr_resources/icon_48.png new file mode 100644 index 0000000..d9cedd9 Binary files /dev/null and b/src/Resources/DigImCorr_resources/icon_48.png differ diff --git a/src/Resources/DigImCorrtt_resources/icon_16.png b/src/Resources/DigImCorrtt_resources/icon_16.png new file mode 100644 index 0000000..f7121b2 Binary files /dev/null and b/src/Resources/DigImCorrtt_resources/icon_16.png differ diff --git a/src/Resources/DigImCorrtt_resources/icon_24.png b/src/Resources/DigImCorrtt_resources/icon_24.png new file mode 100644 index 0000000..0208c3b Binary files /dev/null and b/src/Resources/DigImCorrtt_resources/icon_24.png differ diff --git a/src/Resources/DigImCorrtt_resources/icon_48.png b/src/Resources/DigImCorrtt_resources/icon_48.png new file mode 100644 index 0000000..7d7f650 Binary files /dev/null and b/src/Resources/DigImCorrtt_resources/icon_48.png differ diff --git a/src/dic/Beginner4.m b/src/dic/Beginner4.m new file mode 100644 index 0000000..309748c --- /dev/null +++ b/src/dic/Beginner4.m @@ -0,0 +1,810 @@ +function varargout = Beginner4(varargin) +% BEGINNER4 MATLAB code for Beginner4.fig +% BEGINNER4, by itself, creates a new BEGINNER4 or raises the existing +% singleton*. +% +% H = BEGINNER4 returns the handle to a new BEGINNER4 or the handle to +% the existing singleton*. +% +% BEGINNER4('CALLBACK',hObject,eventData,handles,...) calls the local +% function named CALLBACK in BEGINNER4.M with the given input arguments. +% +% BEGINNER4('Property','Value',...) creates a new BEGINNER4 or raises +% the existing singleton*. Starting from the left, property value pairs are +% applied to the GUI before Beginner4_OpeningFcn gets called. An +% unrecognized property name or invalid value makes property application +% stop. All inputs are passed to Beginner4_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 Beginner4 + +% Last Modified by GUIDE v2.5 23-May-2018 01:44:50 + +% Begin initialization code - DO NOT EDIT +gui_Singleton = 1; +gui_State = struct('gui_Name', mfilename, ... + 'gui_Singleton', gui_Singleton, ... + 'gui_OpeningFcn', @Beginner4_OpeningFcn, ... + 'gui_OutputFcn', @Beginner4_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 Beginner4 is made visible. +function Beginner4_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 Beginner4 (see VARARGIN) + +% Choose default command line output for Beginner4 +handles.output = hObject; +handles.Load_RADIO ='Images'; +axes(handles.axes7);cla;imshow((imread('lines.png')),[]); +axes(handles.axes8);cla;imshow(imread('Shan_Logo1.jpg'),[]); +axes(handles.axes6);cla;imshow(imread('iisc_logo.png'),[]); + +seting_buttons_state(hObject, eventdata, handles,0) +set(handles.push_start_DIC, 'Enable','off'); +set(handles.Status_check,'String','In pixels'); +% set(handles.txt_ref, 'String','Hi'); + +handles.image_crop =0; + +% show_the_image1() +% handles.aarea = 'Default'; +% handles.filename_RADIO ='DEFORMED'; + +handles.val_overlap= 0.5; +handles.val_windowsize_x = 64; +handles.val_windowsize_y = 64; +handles.val_dicmethod = 'method_2'; +% handles.flg1= 0; +handles.displayData = 'U'; +handles.scal =1; +handles.val_scal=1; +handles.Logfile =1; +handles.select_resolution=0; + +axes(handles.axes3);cla;imshow(imread('p1.png') ,[]); +set(handles.checkbox1_LogFile,'Value',1); +seting_buttons_state(hObject, eventdata, handles,0) + +% set(handles.TEXT2,'String','In pixels'); +% set(handles.Status_check,'String','In pixels'); + + + + +% Update handles structurelogfile +guidata(hObject, handles); + + +% UIWAIT makes Beginner4 wait for user response (see UIRESUME) +% uiwait(handles.figure1); + + +% --- Outputs from this function are returned to the command line. +function varargout = Beginner4_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; + + +% --- Executes on button press in LogFile_check. +function LogFile_check_Callback(hObject, eventdata, handles) +% hObject handle to LogFile_check (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of LogFile_check + + +% --- Executes on button press in push_load_reference. +function push_load_reference_Callback(hObject, eventdata, handles) +% hObject handle to push_load_reference (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + + +% --- Executes on button press in push_load_deformed. +function push_load_deformed_Callback(hObject, eventdata, handles) +% hObject handle to push_load_deformed (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + + +% --- Executes on button press in push_startDIC. +function push_startDIC_Callback(hObject, eventdata, handles) +% hObject handle to push_startDIC (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + + +% --- Executes on button press in pushbutton5. +function pushbutton5_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton5 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + + +% --- Executes on button press in checkbox1_LogFile. +function checkbox1_LogFile_Callback(hObject, eventdata, handles) +% hObject handle to checkbox1_LogFile (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of checkbox1_LogFile +handles.Logfile= get(hObject,'Value'); +guidata(hObject, handles) + + + +% --- Executes on button press in push_start_DIC. +function push_start_DIC_Callback(hObject, eventdata, handles) +% hObject handle to push_start_DIC (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +Files = dir(handles.path); +% tmp=length(Files) +[~,Ref_name,EXT] = fileparts(handles.val_filename_1); +tmpp=0; + + + if iscell(handles.def_name)==1 + tmpp = numel(handles.def_name); + FFiles=handles.def_name; + else + tmpp = 1; + FFiles={handles.def_name}; + end + frame_no=tmpp; + + + + + try + [user_response, destination ] = gui_confirm_DICstart2(tmpp,handles.path); + catch + user_response='NO'; + end + + +switch lower(user_response) +case 'no' + % take no action +case 'yes' +seting_buttons_state11(hObject, eventdata, handles,0) +pause(0.1) +% s1= handles.val_filename_1; +% s=strfind(s1,'.');s1(s(end)+1:length(s1)-1); +% + imr11 = imread(handles.val_filename_1); + IM1 = imr11; + imr1=preset_image(imr11); + flies = {}; + M=[]; + + +hhhh = waitbar(0,'Please wait...'); + + + DATA=[]; + kkk=0; + for k =1:tmpp + + waitbar(k/length(Files),hhhh); + + FileNames = FFiles{k}; + file = strcat(handles.path, FFiles{k}); + DefaultORcrop='Full Frame'; + if strcmp(file, handles.val_filename_1)==0 + + + imr22 = imread(file); + imr2=preset_image(imr22); + DefaultORcrop='Full Frame'; + + if handles.image_crop ==1 + DefaultORcrop='Crop'; + yD= handles.Cordinates; + im1 = imr1( yD(2):yD(2)+yD(4),yD(1):yD(1)+yD(3),:); + im2 = imr2( yD(2):yD(2)+yD(4),yD(1):yD(1)+yD(3),:); + IM = IM1( yD(2):yD(2)+yD(4),yD(1):yD(1)+yD(3),:); + + else handles.image_crop ==0; + im1 = imr1; + im2 = imr2; + IM = IM1; + end + + windowsize_x = handles.val_windowsize_x; + windowsize_y= handles.val_windowsize_y; + overlap= handles.val_overlap; + iu_max =32; + iv_max =32; + dt = 1; + dic_type = handles.val_dicmethod; + i_recur = 1; + i_plot=0; + + + + [MM, AA]=DIC_Data_check_Beginner3(im1, im2, windowsize_x, windowsize_y, overlap, iu_max, iv_max,... + handles.scal,FileNames, IM,destination, handles.val_dicmethod,handles.select_resolution, handles.Cordinates, IM1); + + +% figure,imshow(im2,[]), title(strcat('Vector Field:-:', FileNames)) + else + + MM = [0 0 0 0 0 0 0 0]; + dic_type = 'Not-Procssed, Same Image'; + windowsize_x = 0; + overlap=0; + end + + flies{k} = FileNames; + + M=[M;MM]; + AA.name=FileNames; + DATA{k}=AA; + + end + +% % save strcat(destination, 'DATA') DATA; + + + des=strcat(destination,'\DATA'); + save(des,'DATA') + + + if (handles.Logfile==1), myLogfile(destination, 'Image',frame_no,DefaultORcrop, Ref_name, dic_type, ... + windowsize_x,overlap); end + + M = round(M*1000)/1000; +seting_buttons_state11(hObject, eventdata, handles,1) +try +close(hhhh) +catch +end + +Mean_U = M(:,1); +Mean_V = M(:,2); +MIN_U = M(:,3); +MAX_U = M(:,4); +MIN_V = M(:,5); +MAX_V = M(:,6); +StdDev_U = M(:,7); +StdDev_V = M(:,8); + + +% Tabb11 = table(flies', Mean_U, Mean_V, StdDev_U, StdDev_V ); +Tabb11 = table(flies', Mean_U, Mean_V, MIN_U, MAX_U, MIN_V, MAX_V,StdDev_U, StdDev_V ); +writetable(Tabb11,strcat(destination, '\result11_summary' )) + +end +guidata(hObject, handles) + + + + + + +% --- Executes on button press in pushbutton9. +function pushbutton9_Callback(hObject, eventdata, handles) +% hObject handle to pushbutton9 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +open('Help_DIC_BJ2.pdf'); + + +% --- Executes on button press in push_reference. +function push_reference_Callback(hObject, eventdata, handles) +% hObject handle to push_reference (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +Sw= handles.Load_RADIO; + + + +switch Sw + + case {'Images'} + + files= uigetfile_with_preview2(); + + if ~strcmp(files,'') + + seting_buttons(hObject, eventdata, handles); + seting_buttons_state(hObject, eventdata, handles,0) + set(handles.push_deformed, 'Enable','on'); + + + + +% % % seting_buttons_state(hObject, eventdata, handles,1) + s=files; s=s(1:max(find(s=='\'))); pathname=s; + s=files; s=s(max(find(s=='\'))+1:end); filename=s; + handles.val_filename_1 = strcat(pathname, filename); + + + set(handles.txt_ref, 'String','filename'); + guidata(hObject, handles); +% set(handles.push_reference,'String',filename); + set(handles.txt_ref, 'String', filename ); + + tmp = imread(handles.val_filename_1); + + handles.val_file_1 = (tmp); + [~,~,EXT] = fileparts(handles.val_filename_1); + handles.EXT=EXT; + handles.path= pathname; + handles.Cordinates = [1 1 size(handles.val_file_1,2) size(handles.val_file_1,1)]; + + handles.Scal=1; + set(handles.Status_check ,'String',' '); + handles.image_crop =0; + guidata(hObject, handles); + + + files= dir(pathname); + + + axes(handles.axes3);cla;imshow(handles.val_file_1 ,[]); + + set(handles.number_image,'String',''); + end + +% end + + + + case {'Video'} + + + + otherwise + disp('Unknown method') + +end + set(handles.push_start_DIC, 'Enable','off'); + + + + + + + + + +% --- Executes on button press in push_deformed. +function push_deformed_Callback(hObject, eventdata, handles) +% hObject handle to push_deformed (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% % % Files = dir(handles.path); +% tmp=length(Files) +try +[~,Ref_name,EXT] = fileparts(handles.val_filename_1); +catch + EXT = 0; +end +% % % % % % % % % % % % % tmp=0; + + + [filename_def, pathname_def, filterindex_def] = uigetfile( ... + {strcat('*', EXT), strcat(EXT,'-files') }, ... + 'MultiSelect', 'on', 'Selet the DEFORMED images', handles.path); + +handles.def_name= filename_def; +handles.def_path= pathname_def; +handles.correctSelect= filterindex_def; +% % set(handles.push_deformed,'String',filename_def(end)); +im1=imread(handles.val_filename_1); + +if filterindex_def==1 + seting_buttons_state(hObject, eventdata, handles,1) + if iscell(filename_def)==1 + tmp = strcat(num2str(numel(filename_def)), 'files selected'); +% set(handles.push_deformed,'String',filename_def{end}); + im2 = imread(strcat(pathname_def, filename_def{end})); + axes(handles.axes3);cla;imshow(imread(strcat(pathname_def,filename_def{end})) ,[]); + set(handles.push_start_DIC, 'Enable','on'); + else + tmp = strcat('1', 'files selected'); +% set(handles.push_deformed,'String',filename_def); + im2 = imread(strcat(pathname_def, filename_def)); + axes(handles.axes3);cla;imshow(imread(strcat(pathname_def,filename_def)) ,[]); + set(handles.push_start_DIC, 'Enable','on'); + end + set(handles.number_image,'String',tmp); + + if ~(isequal(size(im1),size(im2))) + gui_DIC_sizeCheck1('Title','ERROR'); + tmp = strcat('0', 'files selected'); + set(handles.number_image,'String',tmp); + set(handles.push_start_DIC, 'Enable','off'); + seting_buttons_state(hObject, eventdata, handles,0); + set(handles.push_deformed, 'Enable','on'); + end + + + + +else + tmp = strcat('0', 'files selected'); + set(handles.number_image,'String',tmp); + set(handles.push_start_DIC, 'Enable','off'); + +end + + +seting_buttons(hObject, eventdata, handles); + + +guidata(hObject, handles); + + +% --- Executes when selected object is changed in uibuttongroup5. +function uibuttongroup5_SelectionChangedFcn(hObject, eventdata, handles) +% hObject handle to the selected object in uibuttongroup5 +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +c= get(hObject,'String'); +%%%%%%%%%%%seting_buttons_state(hObject, eventdata, handles,0); +% set(handles.push_reference,'String','Reference'); +% set(handles.push_deformed,'String','Deformed'); + +seting_buttons_state(hObject, eventdata, handles,0); + + switch c + + case {'Images'} + handles.Load_RADIO ='Images'; +% guidata(hObject, handles); + + case 'Video' + handles.Load_RADIO ='Video'; + + EXT='avi'; + [filename_def, pathname_def, filterindex_def] = uigetfile( ... + {strcat('*', EXT), strcat(EXT,'-files') }); + + + if filterindex_def==1 + % seting_buttons_state(hObject, eventdata, handles,1) + % set(handles.number_image,'String',tmp); + + mkdir(pathname_def,'images101'); + + vdo= VideoReader(strcat(pathname_def, '\' ,filename_def)); + ii = 1; + + while hasFrame(vdo) + img = readFrame(vdo); + filename = [sprintf('%03d',ii) '.jpg']; + fullname = fullfile(pathname_def, 'images101',filename); + imwrite(img,fullname) % Write out to a JPEG file (img1.jpg, img2.jpg, etc.) + ii = ii+1; + end + + handles.path= pathname_def; + + handles.Load_RADIO ='Images'; + guidata(hObject, handles) + + else + handles.Load_RADIO ='Images'; + seting_buttons_state(hObject, eventdata, handles,0) + set(handles.radiobutton14, 'Value',0); + set(handles.radiobutton13, 'Value',1); + + guidata(hObject, handles) +% tmp = strcat('0', 'files selected'); +% set(handles.number_image,'String',tmp); +% set(handles.push_start_DIC, 'Enable','off'); + + end + + +% guidata(hObject, handles); +otherwise + disp('Unknown method') + + + end + + seting_buttons(hObject, eventdata, handles); + + guidata(hObject, handles); + + + + function seting_buttons(hObject, eventdata, handles) + + button_list1 =[ handles.radio_scal_full handles.radio_Scal_inpixel ]; + button_list2 =[ handles.radio_scal_crop handles.radio_Scal_inmicron]; + set(button_list1, 'Value',1); + set(button_list2, 'Value',0); + guidata(hObject, handles); + + + function seting_buttons_state11(hObject, eventdata, handles,flg) + + button_list1 =[ handles.radio_scal_full handles.radio_Scal_inpixel handles.push_start_DIC]; + button_list2 =[ handles.radio_scal_crop handles.radio_Scal_inmicron handles.push_deformed]; + button_list3 =[handles.push_reference]; + + if flg==0 + set(button_list1, 'Enable','off'); + set(button_list2, 'Enable','off'); + set(button_list3, 'Enable','off');end + if flg==1 + set(button_list1, 'Enable','on'); + set(button_list2, 'Enable','on'); + set(button_list3, 'Enable','on'); end + + + guidata(hObject, handles); + + + + + + function seting_buttons_state(hObject, eventdata, handles,flg) + + button_list1 =[ handles.radio_scal_full handles.radio_Scal_inpixel handles.push_start_DIC]; + button_list2 =[ handles.radio_scal_crop handles.radio_Scal_inmicron handles.push_deformed]; + + + if flg==0 + set(button_list1, 'Enable','off'); + set(button_list2, 'Enable','off'); end + if flg==1 + set(button_list1, 'Enable','on'); + set(button_list2, 'Enable','on'); end + + + guidata(hObject, handles); + + + + +% --- Executes on button press in radio_scal_full. +function radio_scal_full_Callback(hObject, eventdata, handles) +% hObject handle to radio_scal_full (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of radio_scal_full + + +% --- Executes on button press in radio_Scal_inmicron. +function radio_Scal_inmicron_Callback(hObject, eventdata, handles) +% hObject handle to radio_Scal_inmicron (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of radio_Scal_inmicron + + +% --- Executes during object creation, after setting all properties. +function number_image_CreateFcn(hObject, eventdata, handles) +% hObject handle to number_image (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + + +% --- Executes when selected object is changed in uibuttongroup6. +function uibuttongroup6_SelectionChangedFcn(hObject, eventdata, handles) +% hObject handle to the selected object in uibuttongroup6 +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +handles.aarea= get(hObject,'String'); +% ea= get(hObject,'value') + + + switch handles.aarea + + case {'Crop'} +% handles.val_filename_Crop = imread(handles.val_filename_2); +% imtool(handles.val_filename_Crop ) + handles.image_crop =1; +% handles.val_filename_Crop = handles.val_file_1 ; handles.val_filename_1 + handles.val_filename_Crop = imread(handles.val_filename_1); + + + [~, yD1 ]= gui_Image_Crop_1(handles.val_filename_Crop); + + if ~isempty(yD1) + + if numel(size(handles.val_filename_Crop))==2 + s = size(handles.val_filename_Crop); s(3)=1; + elseif numel(size(handles.val_filename_Crop))==3 + s = size(handles.val_filename_Crop); + end + + yD = cordinate_filter_mex(yD1, s); + + else yD=[]; + end + + + if ~isempty(yD) + yD= round(yD); + handles.Cordinates=yD; + handles.val_filename_Crop = ... + handles.val_filename_Crop( yD(2):yD(2)+yD(4),yD(1):yD(1)+yD(3),:); + + axes(handles.axes3);cla;imshow(handles.val_filename_Crop) + handles.file_1 = handles.val_file_1; + guidata(hObject, handles); + else + seting_buttons(hObject, eventdata, handles); + handles.image_crop =0; + guidata(hObject, handles); + end + + + case {'Full Frame'} + handles.image_crop =0; + handles.val_file_1 = imread(handles.val_filename_1); + axes(handles.axes3);cla;imshow(handles.val_file_1 ) + handles.file_1 = handles.val_file_1; + handles.Cordinates = [1 1 size(handles.val_file_1,2) size(handles.val_file_1,1)]; + + guidata(hObject, handles); + end + + +% --- Executes when selected object is changed in uibuttongroup7. +function uibuttongroup7_SelectionChangedFcn(hObject, eventdata, handles) +% hObject handle to the selected object in uibuttongroup7 +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +cc= get(hObject,'String'); + + switch cc + + case {'In pixels'} + handles.Scal=1; +% guidata(hObject, handles); + tmp3= 'In Pixels'; + handles.select_resolution=0; + guidata(hObject, handles); +% set(handles.TEXT2,'String',tmp3); + set(handles.Status_check,'String',tmp3); + + case 'In micron' + + try + +% BIKI=gui_set_resolution_3(handles.path,handles.EXT); + BIKI=gui_set_resolution_3(1,handles.path,handles.EXT); + catch + BIKI='aadata'; + end + if ~ strcmp(BIKI,'aadata') + BJ = BIKI; + if isstr(BIKI), BIKI=str2num(BIKI); end + if ~isstr(BJ), BJ=num2str(BJ); end + + + + + tmp3= strcat('1 pixel= ',BJ,' micron'); +% set(handles.TEXT2,'String',tmp3); + set(handles.Status_check ,'String',tmp3); + handles.scal=BIKI; + handles.select_resolution=1; + else +% set(handles.Status_check ,'String','NOT Selected', 'DefaultTextColor','green'); + set(handles.radio_Scal_inmicron, 'Value',0); + set(handles.radio_Scal_inpixel, 'Value',1); + set(handles.Status_check ,'String','In Pixel', 'DefaultTextColor','green'); + handles.scal=1; + handles.select_resolution=0; + + end + +% handles.scal=str2num(BIKI) + + + guidata(hObject, handles); + + + + +otherwise + disp('Unknown method.') + + + end + +guidata(hObject, handles) + + +% --- Executes during object creation, after setting all properties. +function Status_check_CreateFcn(hObject, eventdata, handles) +% hObject handle to Status_check (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + + +% --- Executes on button press in radiobutton14. +function radiobutton14_Callback(hObject, eventdata, handles) +% hObject handle to radiobutton14 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of radiobutton14 + + +% --- Executes on button press in radiobutton13. +function radiobutton13_Callback(hObject, eventdata, handles) +% hObject handle to radiobutton13 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hint: get(hObject,'Value') returns toggle state of radiobutton13 + + +% --- Executes during object creation, after setting all properties. +function txt_ref_CreateFcn(hObject, eventdata, handles) +% hObject handle to txt_ref (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + + +% --- Executes during object creation, after setting all properties. +function text21_CreateFcn(hObject, eventdata, handles) +% hObject handle to text21 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + + +% --- Executes during object creation, after setting all properties. +function axes6_CreateFcn(hObject, eventdata, handles) +% hObject handle to axes6 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: place code in OpeningFcn to populate axes6 + + +% --- Executes during object creation, after setting all properties. +function axes7_CreateFcn(hObject, eventdata, handles) +% hObject handle to axes7 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: place code in OpeningFcn to populate axes7 + + +% --- Executes during object creation, after setting all properties. +function axes8_CreateFcn(hObject, eventdata, handles) +% hObject handle to axes8 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: place code in OpeningFcn to populate axes8 diff --git a/src/dic/DIC_Data_check_Beginner3.m b/src/dic/DIC_Data_check_Beginner3.m new file mode 100644 index 0000000..d7f1998 --- /dev/null +++ b/src/dic/DIC_Data_check_Beginner3.m @@ -0,0 +1,150 @@ +function [M, A]=DIC_Data_check_Beginner3(im1, im2, windowsize_x, windowsize_y, overlap, iu_max, iv_max, inum1, FileNames,IM,Hpath, method, pix_mic, cordinates,IM1) + +flg1= 0; +i_recur = 0; + +if pix_mic==0 + inum =inum1; +else + inum =inum1*1e-6; +end + + +if (any( size(im1)~=size(im2)) ) + gui_DIC_sizeCheck1('Title','ERROR'); + flg1=1; +end + + +if overlap > 0.9 + gui_DIC_overlapCheck('Title','ERROR'); + flg1=1; +end + +if (flg1==0) + + switch method + + case 'method_1' + + [xi,yi,iu,iv,da]=DIC_mqd5s( im1', im2', windowsize_x, windowsize_y, ... + overlap, overlap, iu_max, iv_max, i_recur); + + [iu_ft, iv_ft, iu_ip, iv_ip] = DIC_filter2(iu, iv, 2, 1.5, 2, 0 ); + [iu_s, iv_s] = DIC_smooth(iu_ip, iv_ip, 2); + + iu_s(isnan(iu_s))=0; iv_s(isnan(iv_s))=0; + iu_s = iu_s*inum; iv_s = iv_s*inum; % xi = xi*inum; yi = yi*inum; + [x11 y11]= meshgrid(min(xi):max(xi),min(yi):max(yi)); + [a, b]=meshgrid(yi,xi); + F = TriScatteredInterp(b(:),a(:),iu_s(:),'natural'); + U_inp = F(x11,y11); + F = TriScatteredInterp(b(:),a(:),iv_s(:),'natural'); + V_inp = F(x11,y11); + + M =[mean(iu_s(:)) mean(iv_s(:)) ... + min(iu_s(:)) max(iu_s(:)) min(iv_s(:)) max(iv_s(:))... + std2(iu_s) std2(iv_s)]; + + A.iu_s =iu_s; + A.iv_s =iv_s; + A.U_inp =U_inp; + A.V_inp =V_inp; + A.xi = xi; + A.yi = yi; + + + + + case 'method_2' + + + iu_max=windowsize_x; iv_max=windowsize_y; + [xi,yi,iu,iv]=DIC_cor5( im1', im2', windowsize_x, windowsize_y, ... + overlap, overlap, iu_max, iv_max, i_recur); + iu(1,:)=nan; iu(end,:)=nan;iu(:,1)=nan; iu(:,end)=nan; + iv(1,:)=nan; iv(end,:)=nan;iv(:,1)=nan; iv(:,end)=nan; + + [iu_ft, iv_ft, iu_ip, iv_ip] = DIC_filter2(iu, iv, 2, 2.0, 2, 0 ); + [iu_s, iv_s] = DIC_smooth(iu_ip, iv_ip, 2); + iv_s=-iv_s; + iu_s(isnan(iu_s))=0; iv_s(isnan(iv_s))=0; + iu_s = iu_s*inum; iv_s = iv_s*inum; % xi = xi*inum; yi = yi*inum; + [x11 y11]= meshgrid(min(xi):max(xi),min(yi):max(yi)); + [a, b]=meshgrid(yi,xi); + + F = TriScatteredInterp(b(:),a(:),iu_s(:),'natural'); + U_inp = F(x11,y11); + F = TriScatteredInterp(b(:),a(:),iv_s(:),'natural'); + V_inp = F(x11,y11); + +% M =[mean(iu_s(:)) mean(iv_s(:)) std2(iu_s) std2(iv_s)]; + + M =[mean(iu_s(:)) mean(iv_s(:)) ... + min(iu_s(:)) max(iu_s(:)) min(iv_s(:)) max(iv_s(:))... + std2(iu_s) std2(iv_s)]; + + A.iu_s =iu_s; + A.iv_s =iv_s; + A.U_inp =U_inp; + A.V_inp =V_inp; + A.xi = xi; + A.yi = yi; + + + + end + + [xx, yy]= meshgrid(xi,yi); + fig122=figure; set(fig122, 'Visible', 'off'); + imshow(IM*0.5),colormap(gray),hold on, + quiver (xx,yy, iu_s',-iv_s','g'), axis equal; + T=strcat('Vector Field--', FileNames); + title(T); + TT= strcat(Hpath,'\',T); + saveas(gcf, TT); + close(fig122) + + if pix_mic==0, titl='[Pixel]'; + else titl='[m]'; + end + +%%% if ~max((U_inp(:)))==0; ttm=max(abs(U_inp(:)));else ttm=1e-4; end + ttm= max(U_inp(:)); + fig126=figure; set(fig126, 'Visible', 'off'); + imagesc(U_inp,[min(U_inp(:)), ttm]), + colormap(jet), hcb=colorbar;title(hcb,titl), + axis image %Tt=strcat('U--', FileNames); title(Tt); + Tt=strcat('Displacement along X, File-', FileNames); title(Tt); xlabel('X-direction (pixels)'), ylabel('Y-direction (pixels)') + saveas(gcf, strcat(Hpath,'\',Tt)) + close(fig126) + + + image_background2(IM1 , xi,yi,U_inp,cordinates,titl, Tt,Hpath); + + +%%% if ~max(V_inp(:))==0; ttm=max(abs(V_inp(:)));else ttm=1e-4; end + ttm=max(V_inp(:)); + fig125=figure; set(fig125, 'Visible', 'off'); + imagesc(V_inp,[min(V_inp(:)), ttm]); + colormap(jet), hcb=colorbar;title(hcb,titl); + axis image %Tt=strcat('U--', FileNames); title(Tt); + Tt=strcat('Displacement along Y, File-', FileNames); title(Tt); xlabel('X-direction (pixels)'), ylabel('Y-direction (pixels)') + saveas(gcf, strcat(Hpath,'\',Tt)) + close(fig125) + + + + image_background2(IM1 , xi,yi,V_inp,cordinates,titl, Tt, Hpath); + + Tt=strcat('Displacement Magnitude, File-', FileNames); + image_background2(IM1 , xi,yi,sqrt(U_inp.^2+V_inp.^2),cordinates,titl, Tt, Hpath); + + + + + +end + + + diff --git a/src/dic/DIC_Linear_Interp_vector3.m b/src/dic/DIC_Linear_Interp_vector3.m new file mode 100644 index 0000000..979dbe4 --- /dev/null +++ b/src/dic/DIC_Linear_Interp_vector3.m @@ -0,0 +1,81 @@ + +function[ u_inp1 ] = DIC_Linear_Interp_vector3( ui ) +%======================================================================== + +tmp=ui; +[row,col]=find(isnan(ui(2:end-1,2:end-1))); row=row+1; col=col+1; +try +for ii=1:numel(col) + ro=row(ii);cl=col(ii); +% tmp(ro,cl)= (ui(ro-1,cl)+ui(ro+1,cl)+ui(ro,cl-1)+ui(ro,cl+1))/4; + tmp11=[ui(ro-1,cl) ui(ro+1,cl) ui(ro,cl-1) ui(ro,cl+1)]; + chk_nan= ~isnan(tmp11); + tmp11(isnan(tmp11))=0; + tmp(ro,cl)= sum(tmp11.*chk_nan)/ sum(chk_nan(:)); + +end +catch +end + + + + +[row,col]=find(isnan(ui(2:end-1,1))); row=row+1; +try +for ii=1:numel(col) + ro=row(ii);cl=col(ii); + if ~isnan(ui(ro,2)) + tmp(ro,1)= ui(ro,2); + else tmp(ro,1)= tmp(ro,2); + end +end +catch +end + + +[row,col]=find(isnan(ui(2:end-1,end))); row=row+1; +try +for ii=1:numel(col) + ro=row(ii);cl=col(ii); + if ~isnan(ui(ro,end-1)) + tmp(ro,end)= ui(ro,end-1); + else + tmp(ro,end)= tmp(ro,end-1); + end +end +catch +end + +[row,col]=find(isnan(ui(end,2:end-1))); col=col+1; +try +for ii=1:numel(row) + cl=col(ii); + if ~isnan(ui(end-1,cl)) + tmp(end,cl)= ui(end-1,cl); + else + tmp(end,cl)= tmp(end-1,cl); + end +end +catch +end + +[row,col]=find(isnan(ui(1,2:end-1))); col=col+1; +try +for ii=1:numel(row) + ro=row(ii);cl=col(ii); + if ~isnan(ui(2,cl)) + tmp(1,cl)= ui(2,cl); + else + tmp(1,cl)= tmp(2,cl); + end + +end +catch +end + +if isnan(ui(1,1)), tmp(1,1)=tmp(2,2);end +if isnan(ui(end,end)), tmp(end,end)=tmp(end-1,end-1);end +if isnan(ui(1,end)), tmp(1,end)=tmp(2,end-1);end +if isnan(ui(end,1)), tmp(end,1)=tmp(end-1,2);end + +u_inp1=tmp; \ No newline at end of file diff --git a/src/dic/DIC_check_vector.m b/src/dic/DIC_check_vector.m new file mode 100644 index 0000000..dbaf2fc --- /dev/null +++ b/src/dic/DIC_check_vector.m @@ -0,0 +1,25 @@ +function[ iu_f, iv_f, i_cond ] = DIC_check_vector( iu, iv, vec_std, i_filter ) +%======================================================================== +%%======================================================================== + +mx = size(iu,1); +my = size(iu,2); + + +if i_filter == 2 + +[ iu_f, i_cond ] = DIC_vector_filter_median2( iu, vec_std ); +[ iv_f, i_cond ] = DIC_vector_filter_median2( iv, vec_std ); + + +iu_f = iu_f - iv_f*0; +iv_f = iv_f - iu_f*0; + +else + iu_f=iu; + iv_f=iv; + +end + + + diff --git a/src/dic/DIC_cor5.m b/src/dic/DIC_cor5.m new file mode 100644 index 0000000..099c9ec --- /dev/null +++ b/src/dic/DIC_cor5.m @@ -0,0 +1,90 @@ +function [xi, yi, iu, iv] = DIC_cor5( im1, im2, ... + nx_pixel, ny_pixel, ... + overlap_x, overlap_y, ... + iu_max, iv_max, ... + i_mode ) +%======================================================================== +%====================================================================== + + +% set area of search +p_search = 1/3; % percentage of subwindow + +nx = size(im1,1); +ny = size(im1,2); + +pixel_twx = ceil(p_search*nx_pixel); +pixel_twy = ceil(p_search*ny_pixel); + +% (dx_center and dy_center) +[ xi, yi, nx_start, ny_start, nx_overlap, ny_overlap, dx_center, dy_center ] ... + = diccwindowsize32( 'mqd', ... + nx, ny, nx_pixel, ny_pixel, overlap_x, overlap_y ); + +% total number of vectors +mx = max(size(xi)); +my = max(size(yi)); + +if rem(ny_pixel,2)==0, tem=ny_pixel; +else tem=ceil((ny_pixel+ny_pixel-1)/2); end + +for iy = 1: my + + c_proc = strcat( 'process accomplished : ', ... + num2str( 100*(iy-1)/(my-1),' %03.0f' ), '/100' ); + disp( c_proc ) + + for ix = 1: mx + + ix1 = xi(ix) - dx_center; + ix2 = ix1 + nx_pixel - 1; + iy1 = yi(iy) - dy_center; + iy2 = iy1 + ny_pixel - 1; + + wn=nx_pixel/4; + ix11 = xi(ix) - dx_center -wn; if(ix11<1), ix11=1; end + ix22 = ix1 + nx_pixel - 1+ wn; if(ix22>size(im2,1)), ix22=size(im2,1); end + + iy11 = yi(iy) - dy_center -wn; if(iy11<1), iy11=1; end + iy22 = iy1 + ny_pixel - 1 + wn; if(iy22>size(im2,2)), iy22=size(im2,2);end + + f1 = im1( ix1:ix2, iy1:iy2 ); + f2 = im2( ix11:ix22, iy11:iy22 ); + + C3=dic_xcor_m_mex(f1, f2); + [ ip_x3, ip_y3, MMR] = findpeak32_mex( C3, 1); + + SNR = 3.00; + PPR = 1.10; + + tx = round(size(C3,1)/2); + ix_peak = -( ip_y3-tx); + iy_peak = ( ip_x3-tx); + +% +% --- eliminate odd displacement +% + + if (iu_max == 0) & (iv_max == 0) + u_max_dipl = p_search*nx_pixel; + v_max_dipl = p_search*ny_pixel; + else + u_max_dipl = iu_max; + v_max_dipl = iv_max; + end + + if ( abs(ix_peak) > u_max_dipl ) | ( abs(iy_peak) > v_max_dipl ) + ix_peak = NaN; + iy_peak = NaN; + end + + is_x(ix,iy) = ix_peak; + is_y(ix,iy) = iy_peak; + + end +end + iu = -is_x; + iv = is_y; + + + diff --git a/src/dic/DIC_filter2.m b/src/dic/DIC_filter2.m new file mode 100644 index 0000000..86aeee1 --- /dev/null +++ b/src/dic/DIC_filter2.m @@ -0,0 +1,28 @@ +function [iu_ft,iv_ft,iu_ip,iv_ip]=DIC_filter2( iu,iv, ... + i_filter, vec_std,... + i_interp, i_plot) +%======================================================================== +%%======================================================================== +% +% Local fitering to smooth +% + +if i_filter ~= 0 + + [ iu_ft, iv_ft, i_cond ] = DIC_check_vector( iu, iv, vec_std, i_filter ); + +end + + +if i_interp ~= 0 + + [ iu_ip ] = DIC_interp_vector( iu_ft, i_interp ); + [ iv_ip ] = DIC_interp_vector( iv_ft, i_interp ); + +else + + iu_ip = iu_ft; + iv_ip = iv_ft; + +end + diff --git a/src/dic/DIC_interp_vector.m b/src/dic/DIC_interp_vector.m new file mode 100644 index 0000000..c5da4eb --- /dev/null +++ b/src/dic/DIC_interp_vector.m @@ -0,0 +1,11 @@ +function[ uo ] = DIC_interp_vector( ui, ii ) +%======================================================================== + +uo = ui; + +% interpolation + +[ uo ] = DIC_spline_Interp_vector2( ui ); + + + diff --git a/src/dic/DIC_mqd5s.m b/src/dic/DIC_mqd5s.m new file mode 100644 index 0000000..9c50d68 --- /dev/null +++ b/src/dic/DIC_mqd5s.m @@ -0,0 +1,96 @@ +function [xi, yi, iu, iv, D] = DIC_mqd5s( im1, im2, ... + windo_x, windo_y, ... + X_overlap, Y_overlap, ... + u_max, v_max, ... + i_mode ) +%======================================================================== +%% + + +% set filter to eliminate stray vectors: 1-std, 2-median +i_filter = 2; + +% Set threshold value (times of standard deviation) +% Vectors will be eliminated if exceeded +vec_std = 1.5; + +% set interpolation of missing vector : 1-linear, 2-spline, 3-kriging +i_interp = 2; + +% SNR: signal to noise ratio to find peak +r_SNR = 3.0; + +% r_peak: ratio of maximum and mean +r_MMR = 1.10; + +% r_peak: ratio of 1st peak and 2nd peak +r_PPR = 1.1; + +% set area of search +p_search = 1/2; % percentage of subwindow + +% min and max values for MQD +d_min = 10^(-5); +d_max = Inf; + +%======================================================================== + + +% --- initialization +% + +im1 = double(im1); +im2 = double(im2); + +nx = size(im1,1); +ny = size(im1,2); + +% Get the seacrh area +if ( u_max <= 0 ) | ( v_max <= 0 ) + x_maxDis = ceil( p_search*windo_x ); + y_maxDis = ceil( p_search*windo_y ); +else + + x_maxDis = floor( u_max ); + y_maxDis = floor( v_max ); + if x_maxDis >= windo_x + x_maxDis = windo_x; + end + if y_maxDis >= windo_y + y_maxDis = windo_y; + end +end + +X_search = 2*x_maxDis + 1; +Y_search = 2*y_maxDis + 1; + +% to obtain the center locations of all the subwindows +% (dx_center and dy_center) +[ xi, yi, nx_start, ny_start, nx_overlap, ny_overlap, dx_center, dy_center ] ... + = diccwindowsize2( 'mqd', nx, ny, windo_x, windo_y, X_overlap, Y_overlap ); + +dx_center + +dy_center + + +bx=xi-dx_center+windo_x/2; +by=yi-dy_center+windo_y/2; + + +% number of vectors in x and y +mx = max(size(xi)); +my = max(size(yi)); + + + [xi,yi,iu,iv]=DIC_cor5( im1, im2, ... + windo_x, windo_y, ... + X_overlap, Y_overlap, ... + u_max, v_max, ... + i_mode ) +D=[]; + + + + + diff --git a/src/dic/DIC_smooth.m b/src/dic/DIC_smooth.m new file mode 100644 index 0000000..f1c9df0 --- /dev/null +++ b/src/dic/DIC_smooth.m @@ -0,0 +1,9 @@ +function [iu_s,iv_s] = DIC_smooth( iu, iv, i_plot) +%======================================================================== +% +%======================================================================== + +iu_s = DIC_smooth_data(iu); +iv_s = DIC_smooth_data(iv); + +i_plot =i_plot *1; \ No newline at end of file diff --git a/src/dic/DIC_smooth_data.m b/src/dic/DIC_smooth_data.m new file mode 100644 index 0000000..b0917b1 --- /dev/null +++ b/src/dic/DIC_smooth_data.m @@ -0,0 +1,8 @@ +function [y] = DIC_smooth_data(f) +%======================================================================== +%%======================================================================== + +H= [1 2 1; 2 4 2; 1 2 1]; +H= H / sum(H(:)); + +y=imfilter(f,H,'replicate'); \ No newline at end of file diff --git a/src/dic/DIC_spline_Interp_vector2.m b/src/dic/DIC_spline_Interp_vector2.m new file mode 100644 index 0000000..d78beb4 --- /dev/null +++ b/src/dic/DIC_spline_Interp_vector2.m @@ -0,0 +1,19 @@ +function[ tmp1] = DIC_spline_Interp_vector2( ui ) +%======================================================================== +%======================================================================== + +[ u_inp1 ] = DIC_Linear_Interp_vector3( ui ); +tmp=u_inp1; + +[row,col]=find(isnan(ui)); %row=row1+2; col=col1+2; +[X, Y]= meshgrid(1:size(tmp,2),1:size(tmp,1)); +Vq = interp2(X,Y,tmp,row,col); + +tmp1=tmp; + +for ii=1:numel(row) + tmp1(row(ii),col(ii))=Vq(ii); +end + +if any(isnan(tmp1)), tmp1=u_inp1; end + diff --git a/src/dic/DIC_vector_filter_median2.m b/src/dic/DIC_vector_filter_median2.m new file mode 100644 index 0000000..b88d4ec --- /dev/null +++ b/src/dic/DIC_vector_filter_median2.m @@ -0,0 +1,38 @@ +function[ uo, ii ] = DIC_vector_filter_median2( ui, vec_std ) +%======================================================================== +dia=4; +err_int = 0.25; +u_new=zeros(size(ui,1)+dia, size(ui,2)+dia)+nan; +u_new(3:end-2, 3:end-2)=ui; +uo= zeros(size(ui)); +[row, col]=size(ui); +ii=0; +for ri=1:row + for ci=1:col + tmp1 = u_new(ri:ri+4, ci:ci+4); + + + tmp = find(~isnan(tmp1)); + f = tmp1(tmp); + + [f_mean f_std f_med] = dic_filterHist3(f); + + if abs( ui(ri,ci) - f_med ) > (max(err_int, vec_std*f_std)) + uo(ri,ci)=nan; ii=1; + else + uo(ri,ci)=ui(ri,ci); + end + + end +end + + + + + + + + + + + diff --git a/src/dic/cordinate_filter.m b/src/dic/cordinate_filter.m new file mode 100644 index 0000000..1cde78d --- /dev/null +++ b/src/dic/cordinate_filter.m @@ -0,0 +1,14 @@ +function yD = cordinate_filter(yD1, s) +%#codegen + +if ~isempty(yD1) + + yD= round(yD1); + +else + yD=[]; +end + + + + diff --git a/src/dic/dic_filterHist3.m b/src/dic/dic_filterHist3.m new file mode 100644 index 0000000..ced278e --- /dev/null +++ b/src/dic/dic_filterHist3.m @@ -0,0 +1,11 @@ +function[ u_mean, u_std, u_med] = dic_filterHist3( u ) + +limt = 2.0; + +f = u(~isnan(u)); +tm= f.*(abs(f-mean(f))<(limt*std(f))); +tm=tm(find(tm)); + +u_mean=mean(tm); +u_std= std(tm); +u_med= median(tm); diff --git a/src/dic/dic_xcor_m.m b/src/dic/dic_xcor_m.m new file mode 100644 index 0000000..fba35b7 --- /dev/null +++ b/src/dic/dic_xcor_m.m @@ -0,0 +1,217 @@ +function C = dic_xcor_m(varargin) + +%#codegen + +[T, A] = ParseInputs(varargin{:}); + +% Haralick and Shapiro, Volume II (p. 317), +% "Fast Normalized Cross-Correlation", by J. P. Lewis, Industrial Light & Magic. +% + +xcorr_TA = xcorr2_fast(T,A); + +[m, n] = size(T); +mn = m*n; + +local_sum_A = local_sum(A,m,n); +local_sum_A2 = local_sum(A.*A,m,n); + +% Note: diff_local_sums should be nonnegative, but may have negative +% values due to round off errors. Below, we use max to ensure the +% radicand is nonnegative. +diff_local_sums = ( local_sum_A2 - (local_sum_A.^2)/mn ); +denom_A = sqrt( max(diff_local_sums,0) ); + +denom_T = sqrt(mn-1)*std(T(:)); +denom = denom_T*denom_A; +numerator = (xcorr_TA - local_sum_A*sum(T(:))/mn ); + +% We know denom_T~=0 from input parsing; +% so denom is only zero where denom_A is zero, and in +% these locations, C is also zero. +C = zeros(size(numerator)); +tol = sqrt( eps( max(abs(denom(:)))) ); +i_nonzero = find(denom > tol); +C(i_nonzero) = numerator(i_nonzero) ./ denom(i_nonzero); + +% Another numerics backstop. If any of the coefficients are outside the +% range [-1 1], the numerics are unstable to small variance in A or T. In +% these cases, set C to zero to reflect undefined 0/0 condition. +C( ( abs(C) - 1 ) > sqrt(eps(1)) ) = 0; + +%------------------------------- +% Function local_sum +% +function local_sum_A = local_sum(A,m,n) + +% We thank Eli Horn for providing this code, used with his permission, +% to speed up the calculation of local sums. The algorithm depends on +% precomputing running sums as described in "Fast Normalized +% Cross-Correlation", by J. P. Lewis, Industrial Light & Magic. + +B = padarray(A,[m n]); +s = cumsum(B,1); +c = s(1+m:end-1,:)-s(1:end-m-1,:); +s = cumsum(c,2); +local_sum_A = s(:,1+n:end-1)-s(:,1:end-n-1); + +%------------------------------- +% Function xcorr2_fast +% +function cross_corr = xcorr2_fast(T,A) + +T_size = size(T); +A_size = size(A); +outsize = A_size + T_size - 1; + +% figure out when to use spatial domain vs. freq domain +conv_time = time_conv2(T_size,A_size); % 1 conv2 +fft_time = 1*time_fft2(outsize); % 2 fft2 + 1 ifft2 + +cross_corr = conv2(rot90(T,2),A); + + +%------------------------------- +% Function freqxcorr +% +function xcorr_ab = freqxcorr(a,b,outsize) + +% calculate correlation in frequency domain +Fa = fft2(rot90(a,2),outsize(1),outsize(2)); +Fb = fft2(b,outsize(1),outsize(2)); +% xcorr_ab = ifft2(Fa .* Fb,'symmetric'); +xcorr_ab = real(ifft2(Fa .* Fb)); + +%------------------------------- +% Function time_conv2 +% +function time = time_conv2(obssize,refsize) + +% time a spatial domain convolution for 10-by-10 x 20-by-20 matrices + +% a = ones(10); +% b = ones(20); +% mintime = 0.1; + +% t1 = cputime; +% t2 = t1; +% k = 0; +% while (t2-t1)