function MotionAnalyst_tc_multi_points(working_dir)
%clear;
%clc;
%format compact;
%memory;
%pack;
% Author:	Robb Colbrunn, Cleveland Clinic BioRobotics and Mechanical Testing Core
% Date:		2-10-10
% Last change:  5-13-12
% Version:	1.3
% Description: This software automatically tracks 2D motion of points in a
% series of images using a texture correlation algorithm.  

%Revision History
% Revision  Date        Author              Comment
% 1.0                   Mike Bey            Created tc_multi_points.m
% 1.1       2-10-10     Robb Colbrunn       Adapted for use with Motion
                                          % Analyst Suite.  This version 
                                          % switches the working directory of Matlab while
                                          % running.  In addition, it outputs a Result2D.mat
                                          % file for a binary reporting of the texture
                                          % correlation data.     
% 1.2       11-10-11    Nicholas Mealey     Added line 181, which takes all NaN that are
                                          % in corr_matrix and sets them equal to zero to
                                          % keep the program from stopping due to NaNs. 
% 1.3       5-13-12     Robb Colbrunn       Revised to accomodate other
                                          % image adjustments of brightness,
                                          % contrast, and edge detection.
                                          % Includes modifications to control file.
                                          % Added sub-pixel motion
                                          % algorithm to prevent creep.
                                          % Made cross platform compatible
                                          % for file paths. (Mac and Unix)
                                          % Fixed new bug with Matlab 2011 with
                                          % opening files with 'versn'.
%License
%     Copyright (C) 2010-2012  Bey, Colbrunn
% 
%     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 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

viewfigure = 1; % initialize variable to display figure. A latch for closing
% to try to prevent trying to close it after it is already closed
% Read control file (*.csv)
csvpath = [working_dir filesep '*.csv'];
[control_filename,control_pathname] = uigetfile(csvpath,'Select Control File');
fprintf(1,'Reading control file...');
fid = fopen(strcat(control_pathname,control_filename),'rt');
line = fgets(fid);
row = 0;
trial_number = 0;
while (line ~= -1)
    keyword = line(1:findstr(line,',')-1);
    line = line(findstr(line,',')+1:length(line));
    switch lower(keyword)
        case 'directory'
            trial_number = trial_number + 1;
            directory(trial_number,:) = line(1:findstr(line,',')-1);
        case 'imagetype'
            image_type(trial_number,:) = lower(line(1:findstr(line,',')-1));           
        case 'initialxypoints'
            if ~isempty(strfind(line,',,'))
                line = line(1:findstr(line,',,')-1);
            end;
            xy_points_temp = str2num(line);
            number_of_points(trial_number) = length(xy_points_temp) / 2;
            xy_points(trial_number,1:length(xy_points_temp)) = xy_points_temp;
            clear xy_points_temp;
        case 'targetwidth'
            target_width(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'searchwidth'
            search_width(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'interpolate'
            interp(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'numframestodisplay'
            display(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'histogramadjust'
            histo_adjust(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'brightness'
            brightness(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'contrast'
            contrast(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'edgedetect'
            edge_detect(trial_number) = str2num(line(1:findstr(line,',')-1));
        case 'edgethreshold'
            edge_threshold(trial_number) = str2num(line(1:findstr(line,',')-1));
        otherwise
            fprintf(1,'Unknown keyword encountered in parameters file. \n');
    end;      
    line = fgets(fid);
end;
fclose(fid);
fprintf(1,'done\n');
number_of_trials = trial_number;

% figure out how big the xy_points 3D array needs to be by finding the
% directory with the largest # of files to analyze

max_number_of_files = 0;
for trial = 1:number_of_trials
    working_dir = directory(trial,:);
    % include only those files that match the image_type
    switch image_type(trial,:)
        case 'tiff'
            imgpath = [working_dir filesep '*.tiff'];
        case 'tif'
            imgpath = [working_dir filesep '*.tif'];
        case 'jpg'
            imgpath = [working_dir filesep '*.jpg'];
    end;
    dir_info = dir(imgpath);
    number_of_files(trial) = length(dir_info);
    clear dir_info;
end;
max_number_of_files = max(number_of_files);
xy_points(number_of_trials,2*max(number_of_points),max_number_of_files) = 0;


% begin processing all trials
tic;
for trial = 1:number_of_trials
    fprintf(1,'processing trial #%1i...\n',trial);
    % switch to the appropriate data directory
    working_dir = directory(trial,:);
    % get all the filenames from the directory
    switch image_type(trial,:)
        case 'tiff'
            imgpath = [working_dir filesep '*.tiff'];
        case 'tif'
            imgpath = [working_dir filesep '*.tif'];
        case 'jpg'
            imgpath = [working_dir filesep '*.jpg'];
    end;
    dir_info = dir(imgpath);
    number_of_files(trial) = length(dir_info);
        
    % read the "first" image
    fprintf(1,'    processing file %s...\n',dir_info(1).name);
    currentfile1 = [working_dir filesep dir_info(1).name];
    fileinfo=imfinfo(currentfile1);
    switch image_type(trial,:)
        case 'tiff'
            image1 = double(imread(currentfile1,'tiff'));
        case 'tif'
            image1 = double(imread(currentfile1,'tiff'));
        case 'jpg'
            image1 = double(imread(currentfile1,'jpeg'));
        otherwise
            fprintf(1,'Unknown image type specified in parameters file \n');
    end;
    image1 = image1(:,:,1);
   % interpolate the first image if necessary
    if (interp(trial) > 0)
        for i = 1:interp(trial)
            fprintf(1,'       interpolating image...\n');
            image1 = interp2(image1,'cubic');
            search_width(trial) = 2 * search_width(trial);
            target_width(trial) = 2 * target_width(trial);
            xy_points(trial,:) = 2 * xy_points(trial,:);            
        end;
    end;   
    % Apply additional image processing
    image1= MotionAnalyst_AdjustImage(image1, fileinfo.BitDepth, edge_detect(trial), edge_threshold(trial), histo_adjust(trial), brightness(trial), contrast(trial));
        
    % display the first image
    if (display(trial) > 0)
        im1_handle=figure('pos',[500 550 600 450]);
%     for i=1:number_of_points  % enable this to display the targets
%         target_handle(i)=figure('pos',[100*2*(i-1) 100 200 200]);
%     end
        figure(im1_handle);
        imagesc(image1);
        hold on;
        colormap(gray);
        for point = 1:number_of_points(trial)
            x = ((point - 1) * 2) + 1;
            y = x + 1;
            temp_handle = plot(xy_points(trial,x,1),xy_points(trial,y,1),'g.');
            set(temp_handle,'markersize',10);
        end;    
        pause(0.1);
    end;
    % initialize the sub-pixel creep prevention algorithm.  Assure that
    % first step is a move.
    for i=1:number_of_points
        moved(i)=1;
    end
    
    % loop through the remaining files in the directory
    for file_num = 2:length(dir_info)        
        % read the "second" image
        fprintf(1,'    processing file %s...\n',dir_info(file_num).name);
        currentfile2 = [working_dir filesep dir_info(file_num).name];
        switch image_type(trial,:)
            case 'tiff'
                image2 = double(imread(currentfile2,'tiff'));
            case 'tif'
                image2 = double(imread(currentfile2,'tiff'));
            case 'jpg'
                image2 = double(imread(currentfile2,'jpeg'));
        end;
        image2 = image2(:,:,1);
     
        % interpolate the "second" image if necessary
        if (interp(trial) > 0)
            for i = 1:interp(trial)
                fprintf(1,'       interpolating image...\n');
                image2 = interp2(image2,'cubic');
            end;
        end;
        % Apply additional image processing
        image2= MotionAnalyst_AdjustImage(image2,fileinfo.BitDepth, edge_detect(trial), edge_threshold(trial), histo_adjust(trial), brightness(trial), contrast(trial));
        
        % *** HERE IS WHERE WE LOOP THROUGH A VARIABLE NUMBER OF 'POINTS'
        
        %target(target_width(trial),target_width(trial),number_of_points)=0;
        for point = 1:number_of_points(trial)
            % extract the 'target' region from image1
            x = ((point - 1) * 2) + 1;
            y = x + 1;           
            x1 = xy_points(trial,x,file_num-1);
            y1 = xy_points(trial,y,file_num-1);
            
            % sub-pixel creep prevention algorithm. Use new target if
            % moved, use old target if not moved.  This will prevent creep
            % problems where the movement is sub-pixel and when the algorithm 
            % re-learns the new image sub-pixel motion gets tracked as no motion.
            % Interpolate function helps for this, but is not enough in some cases.
            if moved(point)==1
                target(:,:,point) = image1((y1-target_width(trial)):(y1+target_width(trial)),(x1-target_width(trial)):(x1+target_width(trial)));
%                 figure(target_handle(point)); % enable this to display the targets
%                 imagesc(target(:,:,point));
%                 colormap(gray);
            % else keep target as previously stored target
            end
           

            % search for the best match to the 'target'
            corr_matrix = zeros(size(image1));
            for row = y1-search_width(trial):y1+search_width(trial)
                for col = x1-search_width(trial):x1+search_width(trial)
                    image2_region(:,:,point) = image2(row-target_width(trial):row+target_width(trial),col-target_width(trial):col+target_width(trial));
                    corr_matrix(row,col) = corr2(target(:,:,point),image2_region(:,:,point));         
                end;
            end;
            %Sets all NaN in corr_matrix equal to zero
            corr_matrix(isnan(corr_matrix))= 0;
            % find highest correlation for point 1
            temp = ~((corr_matrix - (max(max(corr_matrix)))));
            [I1,J1] = find(temp);
            J1=round(sum(J1)/max(size(J1)));
            I1=round(sum(I1)/max(size(I1)));
            xy_points(trial,x,file_num) = J1;
            xy_points(trial,y,file_num) = I1;

            % figure out if the point moved. 
            if x1==J1 && y1==I1
                moved(point) = 0;
            else
                moved(point) = 1;
            end
        end;      
         
        % display the results on the "second" image
        if (file_num <= display(trial))        
            figure(im1_handle)
            imagesc(image2);
            colormap(gray);
            for point = 1:number_of_points(trial)
                x = ((point - 1) * 2) + 1;
                y = x + 1;
                if sum(moved)==0 %something did not move
                temp_handle = plot(xy_points(trial,x,file_num),xy_points(trial,y,file_num),'r.');
                else %everything moved
                temp_handle = plot(xy_points(trial,x,file_num),xy_points(trial,y,file_num),'g.');
                end
                set(temp_handle,'markersize',10);
            end;    
            pause(0.05);
        else
            pause(0.05);       
            if viewfigure == 1
                close(im1_handle);
                viewfigure =0;
            end
            pause(0.05);
        end;    
    % make the "second" image the "first" image for the next loop
            image1 = image2; 

    end;
        
    % save the results for the current trial
    fprintf(1,'    saving results to output file (results.csv)...\n');
    resultspath = [working_dir filesep 'results.csv'];
    %fid = fopen('results.csv','at');
    fid = fopen(resultspath,'at');
    % write date and time to output file
    fprintf(fid,'%s \n',datestr(now));
    % write column headings to output file
    for point=1:number_of_points(trial)
        temp_string = strcat('x',num2str(point),',');
        fprintf(fid,temp_string);
        temp_string = strcat('y',num2str(point));
        fprintf(fid,temp_string);
        if (point < number_of_points(trial))
            fprintf(fid,',');
        end;
    end;
   
    fprintf(fid,'\n');
    % write data to output file
    for file_num = 1:number_of_files(trial)
        for point = 1:number_of_points(trial)
            x = ((point - 1) * 2) + 1;
            y = x + 1;        
            fprintf(fid,'%4i,%4i',xy_points(trial,x,file_num),xy_points(trial,y,file_num));
            if (point < number_of_points(trial))
                fprintf(fid,',');
            end;
        end;
        fprintf(fid,'\n');
    end;
    fprintf(fid,'\n');
    fclose(fid);
    
    %write to .mat file
    sizes = size(xy_points);
    xypoints = (reshape(xy_points, sizes(2),sizes(3)))';
    Res2DFileHeader = {'Date/Time', datestr(now);
                    'Number of Markers', (sizes(2)/2);
                    'Number of Frames', sizes(3);
                    'Data format', 'x1 y1 x2 y2...';
                    'Data units', 'pixels'};
    resultsmatpath = [working_dir filesep 'results2D.mat'];
    save(resultsmatpath, 'Res2DFileHeader', 'xypoints');
    
    % clear some variables
    clear interm_image_handle;
    clear im1_handle;
    clear image1;
    clear image2;
    clear temp;
end;

% end processing all trials

elapsed_time = toc;
elapsed_time = elapsed_time / 60;
fprintf(1,'Elapsed Time: %1.2f minutes \n',elapsed_time);
