Code covered by the BSD License  

Highlights from
Muscle fascicle tracking - Ultrasound

image thumbnail

Muscle fascicle tracking - Ultrasound

by

 

02 Sep 2011 (Updated )

Implementation of an optical flow algorithm to track muscle length changes imaged with ultrasound.

muscle_track_v1(varargin)
function varargout = muscle_track_v1(varargin)
% MUSCLE_TRACK_V1 M-file for muscle_track_v1.fig
%      MUSCLE_TRACK_V1, by itself, creates a new MUSCLE_TRACK_V1 or raises the existing
%      singleton*.
%
%      H = MUSCLE_TRACK_V1 returns the handle to a new MUSCLE_TRACK_V1 or the handle to
%      the existing singleton*.
%
%      MUSCLE_TRACK_V1('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MUSCLE_TRACK_V1.M with the given input arguments.
%
%      MUSCLE_TRACK_V1('Property','Value',...) creates a new MUSCLE_TRACK_V1 or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before muscle_track_v1_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to muscle_track_v1_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 muscle_track_v1
%
% Last Modified by GUIDE v2.5 26-Aug-2011 15:11:06
%
% Written by Dr Glen Lichtwark (The Univesity of Queensland)
% Copyright Dr Glen Lichtwark 2011.
% 
% Please cite the following manuscripts in any academic publications arising
% from this contribution -
% 
% Jarred G. Gillett, Rod S. Barrett, Glen A. Lichtwark 
% Reliability and accuracy of an automatic tracking algorithm to measure 
% passive and active muscle fascicle length changes from ultrasound.  
% 
% Cronin, NJ, Carty, CP, Barrett, RS and Lichtwark, G. (In Press). 
% Automatic tracking of medial gastrocnemius fascicle length during human 
% locomotion. Journal of Applied Physiology.
% doi:10.?1152/?japplphysiol.?00530.?2011

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @muscle_track_v1_OpeningFcn, ...
                   'gui_OutputFcn',  @muscle_track_v1_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 muscle_track_v1 is made visible.
function muscle_track_v1_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 muscle_track_v1 (see VARARGIN)

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

set(0,'RecursionLimit',3000)
% set some standard values for processing

% default depth setting - please change if you have a standard 
% depth you want to use all of the time
handles.ID = 60; 
% default settings for affine flow tracking
handles.SIGMA = 3;
handles.S_STEP = 3;
% try prevent flickering in figure
set(gcf,'DoubleBuffer','on');

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes muscle_track_v1 wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = muscle_track_v1_OutputFcn(~, 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;

% --------------------------------------------------------------------
function File_Callback(hObject, eventdata, handles)
% hObject    handle to File (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function Image_Callback(hObject, eventdata, handles)
% hObject    handle to Image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --------------------------------------------------------------------
function Tracking_Callback(hObject, eventdata, handles)
% hObject    handle to Tracking (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Fascicle_Callback(hObject, eventdata, handles)
% hObject    handle to Fascicle (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 cut_frames_before.
function cut_frames_before_Callback(hObject, eventdata, handles)
% hObject    handle to cut_frames_before (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % this function effectively trims the video sequence in case there are
    % frames you don't want to analyse in the video - trims frames before
    % the current frame
    frame_no = round(get(handles.frame_slider,'Value'));
    
    handles.start_frame = frame_no + handles.start_frame;
       
    handles.NumFrames = handles.movObj.NumberOfFrames-handles.start_frame+1;
    
    % reset the slider bar
    set(handles.frame_slider,'Min',1);
    set(handles.frame_slider,'Max',handles.NumFrames);
    set(handles.frame_slider,'Value',1);
    set(handles.frame_slider,'SliderStep',[1/handles.NumFrames 5/handles.NumFrames]);

    handles.fas_x(1:handles.start_frame-1)=[];
    handles.fas_y(1:handles.start_frame-1)=[];
    
    % set the string in the frame_number to 1
    set(handles.frame_number,'String',1);

    % update the image axes using show_image function (bottom)
    show_image(hObject,handles)
end

% --- Executes on button press in cut_frames_after.
function cut_frames_after_Callback(hObject, eventdata, handles)
% hObject    handle to cut_frames_after (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % this function effectively trims the video sequence in case there are
    % frames you don't want to analyse in the video - trims frams after
    % current frame
    frame_no = round(get(handles.frame_slider,'Value'));
    
    handles.NumFrames = frame_no;   
    
    % reset the slider bar
    set(handles.frame_slider,'Min',1);
    set(handles.frame_slider,'Max',handles.NumFrames);
    set(handles.frame_slider,'Value',handles.NumFrames);
    set(handles.frame_slider,'SliderStep',[1/handles.NumFrames 5/handles.NumFrames]);
    
    handles.fas_x(frame_no+1:end)=[];
    handles.fas_y(frame_no+1:end)=[];
    
    % set the string in the frame_number to 1
    set(handles.frame_number,'String',frame_no);

    % update the image axes using show_image function (bottom)
    show_image(hObject,handles)
end

% --- Executes on slider movement.
function frame_slider_Callback(hObject, eventdata, handles)
% hObject    handle to frame_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

% get the current value from the slider (round to ensure it is integer)
frame_no = round(get(handles.frame_slider,'Value'));

% set the string in the frame_number box to the current frame value
set(handles.frame_number,'String',num2str(frame_no));

% update the image axes using show_image function (bottom)
show_image(hObject,handles)

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

% Hint: slider controls usually have a light gray background.


if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end


function frame_number_Callback(hObject, eventdata, handles)
% hObject    handle to frame_number (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of frame_number as text
%        str2double(get(hObject,'String')) returns contents of frame_number as a double

% goto frame number
frame_no = str2num(get(handles.frame_number,'String'));
set(handles.frame_slider,'Value',round(frame_no));

% update the image axes using show_image function (bottom)
show_image(hObject,handles)

% --- Executes during object creation, after setting all properties.
function frame_number_CreateFcn(hObject, eventdata, handles)
% hObject    handle to frame_number (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% 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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% -----------------------------------------------
% -----------------------------------------------
function goto_next_frame(hObject, handles)

if get(handles.autorun_but,'BackgroundColor') == [0 1 0]

% get the current value from the slider (round to ensure it is integer)
frame_no = round(get(handles.frame_slider,'Value'))+1;

if frame_no < get(handles.frame_slider,'Max')
    % set the slider
    set(handles.frame_slider,'Value',frame_no);
    % set the string in the frame_number box to the current frame value
    set(handles.frame_number,'String',num2str(frame_no));
else  handles.run = 0;
    set(handles.autorun_but,'BackgroundColor',[1 0 0]);
end

% update the image axes using show_image function (bottom)
show_image(hObject,handles)

end


% --- Executes on button press in define_roi.
function define_roi_Callback(hObject, eventdata, handles)
% hObject    handle to define_roi (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % define a region of interest in the current image using the roipoly
    % function from the image processing toolbox
    axes(handles.axes1)    
    
    [handles.ROI,handles.ROIx, handles.ROIy] = roipoly;
    
    handles.CIm = handles.Im;
    
    % show the ROI on the image
    show_image(hObject,handles)
    
end

% --- Executes on button press in define_fascicle.
function define_fascicle_Callback(hObject, eventdata, handles)
% hObject    handle to define_fascicle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % find current frame number from slilder
    frame_no = round(get(handles.frame_slider,'Value'));
    
    if ~isempty(handles.fas_x{frame_no})   
        handles.fas_x{frame_no} = [];
        handles.fas_y{frame_no} = [];
    end
    
    % make the image axes the current axes
    axes(handles.axes1)
    imshow(handles.Im)
      
    % select the two points that define the fascicle line of action -
    % image this with the rbline function
    [p1,p2]=rbline(handles.axes1);
    
    % store these points in specific locations
    handles.fas_x{frame_no}(1) = p1(1); 
    handles.fas_x{frame_no}(2) = p2(1);
    handles.fas_y{frame_no}(1) = p1(2); 
    handles.fas_y{frame_no}(2) = p2(2);

    % correct which is the upper and with is the lower points relative to
    % the image - to standardise
    if handles.fas_y{frame_no}(1) > handles.fas_y{frame_no}(2)
        handles.fas_y{frame_no} = fliplr(handles.fas_y{frame_no});
        handles.fas_x{frame_no} = fliplr(handles.fas_x{frame_no});
    end
    
    % make the  fascicle end points the current points for subsequent
    % tracking
    handles.current_xy = [handles.fas_x{frame_no};handles.fas_y{frame_no}]';
    
    % calculate the length and pennation for the current frame
    handles.fas_pen(frame_no) = atan2(abs(diff(handles.fas_y{frame_no})),abs(diff(handles.fas_x{frame_no})));
    scalar = handles.ID/handles.movObj.Height;
    handles.fas_length(frame_no) = scalar*sqrt(diff(handles.fas_y{frame_no}).^2 + diff(handles.fas_x{frame_no}).^2);
    
    handles.CIm = handles.Im;
    
    % show the fascicle on the image
    show_image(hObject,handles)
    
end

% --- Executes on button press in clear_fascicle.
function clear_fascicle_Callback(hObject, eventdata, handles)
% hObject    handle to clear_fascicle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % clear the current fascicle points to start again
    if isfield(handles,'fas_length')
        handles = rmfield(handles,'fas_length');
        handles = rmfield(handles,'fas_pen');
        handles = rmfield(handles,'current_xy');
        for i = 1:handles.movObj.NumberOfFrames
            handles.fas_x{i}=[];
            handles.fas_y{i}=[];
        end
    end
    
    if isfield(handles,'ROI')
        handles = rmfield(handles,'ROI');
    end
               
    % set current frame to 1
    set(handles.frame_slider,'Value',1);
    set(handles.frame_number,'String',1);
    
    show_image(hObject,handles)
end


% --- Executes on button press in process_all.
function process_all_Callback(hObject, eventdata, handles)
% hObject    handle to process_all (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj') && isfield(handles,'current_xy') && isfield(handles,'ROI')
    
    % perform the affine flow tracking on the muscle ultrasound for all
    % frames from current frame
    
    % find current frame number from slilder
    frame_no = round(get(handles.frame_slider,'Value'));
    
    if isempty(handles.fas_length(frame_no))
        f = warndlg('First frame must have fascicle and ROI marked', 'Warning');
    else
        % create image from movie frame
        Im = read(handles.movObj, handles.start_frame+frame_no-1);

        % crop image with current cropping value
        Im = imcrop(Im,handles.crop_rect);
        
        % set image as current image
        handles.CIm = Im;    

        % define end frame to get o
        end_frame = get(handles.frame_slider,'Max');
        
        h = waitbar(0,'Please wait while processing...');  
        
        image_depth = handles.ID;
        
        for i = frame_no+1:end_frame
            
            % create image from movie frame
            Im = read(handles.movObj, handles.start_frame+i-1);

            % crop image with current cropping value
            Im = imcrop(Im,handles.crop_rect);

            handles.Im = Im;    
            
            % set new image
            handles.NIm = handles.Im;

            % convert current and new images to greyscale
            im1 = double(rgb2gray(handles.CIm))/256;
            im2 = double(rgb2gray(handles.NIm))/256;             

            sigmaXY = handles.SIGMA;
            sampleStep = handles.S_STEP;
            
            % determine affine flow on region of interest (roi) and convert
            % this into a flow warp to be applied to fascicle end points
            af = affine_flow('image1', im1, 'image2', im2, 'sigmaXY',sigmaXY,'sampleStep',sampleStep,'regionOfInterest',handles.ROI);
            af = af.findFlow;
            flow = af.flowStruct;
            w = affine_flow.warp(flow);
            
            % apply flow warp (affine transformation) to the fascicle end
            % points of current frame
            pos = [handles.fas_x{i-1}(:) handles.fas_y{i-1}(:) ones(length(handles.current_xy),1)] * w;
            
            % save these as the fascicle end points of new frame
            handles.fas_x{i}(1:2) = pos(:,1);
            handles.fas_y{i}(1:2) = pos(:,2);   
                                                   
            % calculate the length and pennation for the new frame
            handles.fas_pen(i) = atan2(abs(diff(handles.fas_y{i})),abs(diff(handles.fas_x{i})));
            scalar = image_depth/handles.movObj.Height;
            handles.fas_length(i) = scalar*sqrt(diff(handles.fas_y{i}).^2 + diff(handles.fas_x{i}).^2);
            
            % make the new frame the current frame for the next iteration
            handles.CIm = handles.NIm;
            
            % update wait bar
            waitbar(i/end_frame,h)
            
        end
        %close waitbar
        close(h)
        
        set(handles.frame_slider,'Value',frame_no);
        set(handles.frame_number,'String',num2str(frame_no));
        
        % update the image axes using show_image function (bottom)
        show_image(hObject,handles)
        
      end
end

% Update handles structure
guidata(hObject, handles);


% --------------------------------------------------------------------
function menu_affine_settings_Callback(hObject, eventdata, handles)
% hObject    handle to menu_affine_settings (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    
    % change the two important parameters used in the affine flow algorithm
    if ~isfield(handles, 'SIGMA')
        handles.SIGMA = 3;
    end
    
    if ~isfield(handles, 'S_STEP')
        handles.SIGMA = 3;
    end
    
    % use a dialog box to get the two variables
    A = inputdlg({'Sigma','Sample Step'}, 'Affine Flow Settings', 1,...
        {num2str(handles.SIGMA),num2str(handles.S_STEP)});
    handles.SIGMA = str2double(A{1});
    handles.S_STEP = str2double(A{2});
    
    % Update handles structure
    guidata(hObject, handles);

% --------------------------------------------------------------------
function menu_crop_image_Callback(hObject, eventdata, handles)
% hObject    handle to menu_crop_image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % define current axes
    axes(handles.axes1)

    % use imcrop tool to determine croppable area
    [~,handles.crop_rect] = imcrop;
    
    handles.crop_rect = round(handles.crop_rect);
       
    % update the image axes using show_image function (bottom)
    show_image(hObject,handles)
    
end

% --------------------------------------------------------------------
function menu_reset_image_Callback(hObject, eventdata, handles)
% hObject    handle to menu_reset_image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % set the image croppable area to the maximum area
    handles.crop_rect = [1 1 handles.movObj.Width handles.movObj.Height];
    
    % update the image axes using show_image function (bottom)
    show_image(hObject,handles)
    
end


% --------------------------------------------------------------------
function menu_set_depth_Callback(hObject, eventdata, handles)
% hObject    handle to menu_set_depth (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    if ~isfield(handles, 'ID')
        handles.ID = 60.4;
    end
    
    A = inputdlg('Enter image depth', 'Image Depth', 1, {num2str(handles.ID)});
    handles.ID = str2double(A{1});
    
    % Update handles structure
    guidata(hObject, handles);
    
end


% --------------------------------------------------------------------
function menu_save_tracking_Callback(hObject, eventdata, handles)
% hObject    handle to menu_save_tracking (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj') && isfield(handles,'fas_length')
    
    dt = 1/handles.FrameRate;
    time = (dt:dt:length(handles.fas_length)*dt)+((handles.start_frame-1)*dt);
        
    %determine any non-zero entries in fascicle length array 
    nz = logical(handles.fas_length ~= 0);
    
    T = time(nz)';
    FL = handles.fas_length(nz);
    PEN = handles.fas_pen(nz);   
    
    [filename, pathname] = uiputfile('*.txt', 'Save Fasicle Tracking as');

    % open a new text file for writing
    fid = fopen([pathname filename],'w');

    fprintf(fid,'TIME\tFASCICLE LENGTH\tPENNATION ANLE\n');
    fprintf(fid,'%3.4f\t%3.6f\t%3.6f\n',[T FL' PEN']');

    fclose(fid);
    
end

% --------------------------------------------------------------------
function menu_load_fascicle_Callback(hObject, eventdata, handles)
% hObject    handle to menu_load_fascicle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
            
    % load a fascicle from a previously saved MAT file with the fascicle
    % end points - from menu_save_fascicle function
    
    % find current frame number from slilder
    frame_no = round(get(handles.frame_slider,'Value'));
    
    if isfield(handles,'fas_x')
        [fname, pname] = uigetfile('*.mat','Load tracking MAT file');
        a = load([pname fname]);
                
        handles.fas_x{frame_no} = a.fas_x;
        handles.fas_y{frame_no} = a.fas_y;
        
        handles.current_xy(1,1) = handles.fas_x{frame_no}(1);
        handles.current_xy(1,2) = handles.fas_y{frame_no}(1);
        handles.current_xy(2,1) = handles.fas_x{frame_no}(2);
        handles.current_xy(2,2) = handles.fas_y{frame_no}(2);
              
        % update the image axes using show_image function (bottom)
        show_image(hObject,handles)
    end       
    
end

% --------------------------------------------------------------------
function menu_save_fascicle_Callback(hObject, eventdata, handles)
% hObject    handle to menu_save_fascicle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % save the current fascicle end points in a MAT file so that the
    % fascicle can be re-loaded at any time
    
    % find current frame number from slilder
    frame_no = round(get(handles.frame_slider,'Value'));
        
    if isfield(handles,'fas_x_corr')
        
        fas_x = handles.fas_x_corr{frame_no};
        fas_y = handles.fas_y_corr{frame_no};
        
        [fname, pname] = uiputfile('*.mat', 'Save tracking MAT file');
        
        save([pname fname],'fas_x','fas_y');
        
    elseif isfield(handles,'fas_x')
        
        fas_x = handles.fas_x{frame_no};
        fas_y = handles.fas_y{frame_no};
        
        [fname, pname] = uiputfile('*.mat', 'Save tracking MAT file');
        
        save([pname fname],'fas_x','fas_y');
        
    end
    
end


% --------------------------------------------------------------------
function menu_save_image_Callback(hObject, eventdata, handles)
% hObject    handle to menu_save_image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if isfield(handles,'movObj')
    
    % save the current frame as an image including any fascicles and/or ROI
    % drawings
    [fileout, pathout, FI] = uiputfile('*.tif', 'Save video as');
    
    if FI > 0
    IM_out = getframe(handles.axes1);
    imwrite(IM_out.cdata,[pathout fileout], 'TIFF');
    end
    
end

% --------------------------------------------------------------------
function menu_load_file_Callback(hObject, eventdata, handles)
% hObject    handle to menu_load_file (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% load a new video file 

% if a video is already opened then delete current video and any current
% fascicles etc
if isfield(handles,'movObj')
    handles = rmfield(handles,'mov');
    handles = rmfield(handles,'movObj');
end

if isfield(handles,'fas_length')
    handles = rmfield(handles,'fas_length');
    handles = rmfield(handles,'fas_pen');
    handles = rmfield(handles,'fas_x');
    handles = rmfield(handles,'fas_y');
end

version = ver;
% do some checking to make sure the image processing toolbox is available
if isempty(strfind([version.Name],'Image Processing Toolbox'))
   error('This application requires the Matlab Image Processing Toolbox to be installed') 
end

% determine the release date and use the appropriate file loading function 
% newer versions > 2010b use VideoReader, older versions use mmreader
if datenum(version(1).Date) >= 734353
    reader = 'VideoReader';
else reader = 'mmreader';
end

% if post 2010a version of Matlab then look at the available video formats
% for the list
if datenum(version(1).Date) >= 734202 
    %determine the available video formats 
    file_formats = eval([reader '.getFileFormats']);
    format_list{1,1} = [];
    format_list{1,2} = 'All Video Files (';
    for i = 1:length(file_formats)
        format_list{1,1} = [format_list{1,1} '*.' file_formats(i).get.Extension ';'];
        format_list{1,2} = [format_list{1,2} file_formats(i).get.Description ', '];
        format_list{i+1,1} = ['*.' file_formats(i).get.Extension];
        format_list{i+1,2} = ['*.' file_formats(i).get.Extension ' - ' file_formats(i).get.Description];
    end

    format_list{1,2} = [format_list{1,2}(1:end-2) ')'];

    %load the avi file
    [handles.fname, handles.pname] = uigetfile(format_list, 'Pick a movie file');

else % pre 2010a mmreader function cannot use the getFileFormats method so only allow AVI files to be selected
    [handles.fname, handles.pname] = uigetfile('*.avi', 'Pick a movie file');
end

if isequal(handles.fname,0) || isequal(handles.pname,0)
       return;
end

cd (handles.pname)    

handles.movObj = eval([reader '([handles.pname handles.fname])']);

handles.NumFrames = handles.movObj.NumberOfFrames;
vidHeight = handles.movObj.Height;
vidWidth = handles.movObj.Width;

handles.FrameRate = handles.movObj.NumberOfFrames/handles.movObj.Duration;

% display the path and name of the file in the filename text box
set(handles.filename,'String',[handles.pname handles.fname])

% set the string in the frame_number box to the current frame value (1)
set(handles.frame_number,'String',num2str(1))

handles.start_frame = 1;

% set the limits on the slider - use number of frames to set maximum (min =
% 1)
set(handles.frame_slider,'Min',1);
set(handles.frame_slider,'Max',handles.NumFrames);
set(handles.frame_slider,'Value',1);
set(handles.frame_slider,'SliderStep',[1/handles.NumFrames 10/handles.NumFrames]);

% set the image croppable area to the maximum area
if ~isfield(handles,'crop_rect')
    handles.crop_rect = [1 1 vidWidth vidHeight];
end

for i = 1:handles.NumFrames
    handles.fas_x{i}=[];
    handles.fas_y{i}=[];
end

cd(handles.pname)

% make a timeline which corresponds to the ultrasound frame data
handles.Time = (1/handles.FrameRate:1/handles.FrameRate:handles.NumFrames/handles.FrameRate)';
    
% update the image axes using show_image function (bottom)
show_image(hObject,handles)


%---------------------------------------------------------
% Function to show image with appropriate image processing
%---------------------------------------------------------
function show_image(hObject, handles)

if isfield(handles,'movObj')
                        
    % find current frame number from slilder
    frame_no = round(get(handles.frame_slider,'Value'));
       
    image_depth = handles.ID;
    
    handles.mov.cdata = read(handles.movObj, frame_no+handles.start_frame-1);

    % create image from movie frame
    Im = handles.mov.cdata;
    
    % crop image with current cropping value
    Im = imcrop(Im,handles.crop_rect);
    
    % make the image axes the current axes
    axes(handles.axes1)
    cla
    
    % show the image
    imshow(Im)
    
    handles.Im = Im;    
    
    handles.NIm = handles.Im;
    
    % plot the ROI if available
    if isfield(handles,'ROI')
        hold on
        plot(handles.ROIx,handles.ROIy,'r:','LineWidth',2)
        hold off
    end
    
    if isfield(handles,'current_xy') 
     
        % if the new frame hasn't been tracked yet and a fascicle for the
        % previous viewed frame is determined along with a ROI defined then
        % perform the tracking to the new frame
        if isempty(handles.fas_x{frame_no}) && isfield(handles,'ROI')
        
            % convert current and new images to greyscale
            im1 = double(rgb2gray(handles.CIm))/256;
            im2 = double(rgb2gray(handles.NIm))/256;
            
            sigmaXY = handles.SIGMA;
            sampleStep = handles.S_STEP;

            % determine affine flow on region of interest (roi) and convert
            % this into a flow warp to be applied to fascicle end points
            af = affine_flow('image1', im1, 'image2', im2, 'sigmaXY',sigmaXY,'sampleStep',sampleStep,'regionOfInterest',handles.ROI);
            af = af.findFlow;
            flow = af.flowStruct;
            w = affine_flow.warp(flow);
            
            % apply flow warp (affine transformation) to the fascicle end
            % points of current frame
            pos = [handles.current_xy(:,1) handles.current_xy(:,2) ones(length(handles.current_xy),1)] * w;
            
            % save these as the fascicle end points of new frame and the
            % current points for the next iteration
            handles.current_xy = pos;

            handles.fas_x{frame_no}(1) = handles.current_xy(1,1);
            handles.fas_y{frame_no}(1) = handles.current_xy(1,2);
            handles.fas_x{frame_no}(2) = handles.current_xy(2,1);
            handles.fas_y{frame_no}(2) = handles.current_xy(2,2);
            
            % make the new frame the current frame for the next iteration
            handles.CIm = handles.NIm;
            
            % calculate the length and pennation for the current frame
            handles.fas_pen(frame_no) = atan2(abs(diff(handles.fas_y{frame_no})),abs(diff(handles.fas_x{frame_no})));
            scalar = image_depth/handles.movObj.Height;
            handles.fas_length(frame_no) = scalar*sqrt(diff(handles.fas_y{frame_no}).^2 + diff(handles.fas_x{frame_no}).^2);
                    
        else % if the tracking has been performed already then just display 
            % the fascicle for this frame and make the fascicle end points
            % for this frame current
            handles.current_xy(1,1) = handles.fas_x{frame_no}(1);
            handles.current_xy(1,2) = handles.fas_y{frame_no}(1);
            handles.current_xy(2,1) = handles.fas_x{frame_no}(2);
            handles.current_xy(2,2) = handles.fas_y{frame_no}(2);
            handles.fas_pen(frame_no) = atan2(abs(diff(handles.fas_y{frame_no})),abs(diff(handles.fas_x{frame_no})));
            scalar = image_depth/handles.movObj.Height;
            handles.fas_length(frame_no) = scalar*sqrt(diff(handles.fas_y{frame_no}).^2 + diff(handles.fas_x{frame_no}).^2);
        end
        
        % display the current fascicle on the current frame image
        hold on
        plot(handles.axes1,handles.fas_x{frame_no},handles.fas_y{frame_no},'ro-','LineWidth',2)
        hold off
                        
    end
        
    % Update handles structure
    guidata(hObject, handles);
              
end

Contact us