function compare_result = compare_contours(dicom_path,manual_contour_path,auto_contour_path,para)

%COMPARE_CONTOURS Compare manual rawn contours with auto contours

%   COMPARE_CONTOURS(DICOM_PATH,MANUAL_CONTOUR_PATH,AUTO_CONTOUR_PATH)

%

%   Copyright: Imaging Research, Sunnybrook Health Sciences Centre, Toronto, ON, Canada

%   Author: Perry Radau and Yingli Lu

%   Email: perry.radau@gmail.com

%   Version: 1.0

%   Date: 2009/06/29

%   OUTPUT: compare_result is a struct with fields of

%-----------------------------------------------------

%   auto_number_i:  total number of auto inside contours

%   auto_number_o:  total number auto outside contours;

%   manual_number_i: total number manual inside contours

%   manual_number_o: total number manual outside contours

%   detect_percent_i: percent of detected number of auto inside contours compared to total number of  manual inside contours

%   detect_percent_o: percent of detected number of auto outside contours compared to total number of  manual outside contours

%   good_percent_i: percent of good auto inside contours (a good contour has an average distance smaller than para.dist_limit)

%   good_percent_o: percent of good auto outside contours

%   auto_ef_pic: auto contour's ejection fraction, '_pic' means Papillary Included in the LV Cavity

%   auto_ef_pim: auto contour's ejection fraction, '_pim' means Papillary Included in the Myocardium

%   auto_lvm_pic: auto contour's lv mass

%   auto_lvm_pim: 

%   manual_ef_pic: auto contour's ejection fraction

%   manual_ef_pim: 

%   manual_lvm_pic: manual contour's lv mass

%   manual_lvm_pim: 

%   avg_dist_i: average distance of inside contours

%   avg_dist_o: average distance of outside contours

%   avg_dm_i: average dice metric of inside contours

%   avg_dm_o: average dice metric of outside contours

%   auto_missing_index_i: missing auto inside contours' index

%   auto_missing_index_o: missing auto outside contours' index

%   auto_bad_index_i: bad auto inside contours' index (a bad contour has an  average distance larger than para.dist_limit)

%   auto_bad_index_o: bad auto outside contours' index 

%-----------------------------------------------------

%-initialize

compare_result = [];

%-check path

if ~exist(dicom_path,'dir')

    disp([dicom_path ' :NOT exist!'])

    return;

end

if ~exist(manual_contour_path,'dir')

    disp([manual_contour_path ' :NOT exist!'])

    return;

end

if ~exist(auto_contour_path,'dir')

    disp([auto_contour_path ' :NOT exist!'])

    return;

end

%-check files

dicom_files = dir([dicom_path filesep '*.dcm']);

if isempty(dicom_files)

    disp([dicom_path ' :NO dicom files!'])

    return;

end

manual_contour_files = dir([manual_contour_path filesep '*.txt']);

if isempty(manual_contour_files)

    disp([manual_contour_path ' :NO manual contours!'])

    return;

end

auto_contour_files = dir([auto_contour_path filesep '*.txt']);

if isempty(auto_contour_files)

    disp([auto_contour_path ' :NO auto contours!'])

    return;

end

%-dicominfo

try

  dicom_filename = dicom_files(1).name; %use the first dicom file.

  full_dicom_filename = [dicom_path filesep dicom_filename];

  dicom_meta= dicominfo(full_dicom_filename);

  para.width = dicom_meta.Width;%image width

  para.height = dicom_meta.Height;%image height

  para.pixel_spacing = dicom_meta.PixelSpacing; % mm

  para.thickness = dicom_meta.SliceThickness;% mm

  para.gap = dicom_meta.SpacingBetweenSlices - para.thickness; % mm

  para.phase_number = dicom_meta.CardiacNumberOfImages; %numer of phases

  para.image_number = length(dicom_files);%number of images

catch

  s = lasterror;

  disp(s.message);

  return

end

%-calc percent of auto contours compared to manual contours

manual_number_i = 0; %number of manual inside contours

manual_number_o = 0; %number of manual outside contours

manual_index_i = []; %index of manual inside contours

manual_index_o = []; %index of manual outside contours

%count manual contours

for ix = 1:para.image_number 

    sindex = add_zero_index(ix,para.digit_length);

    %inside contour

    full_icontour_filename = get_contour_filename(manual_contour_path, para.name_prefix, sindex,para.inside_contour_mode,para.manual_seg_mode);

    if exist(full_icontour_filename,'file') 

       fileinfo_i = dir(full_icontour_filename);

       if fileinfo_i.bytes ~= 0 % not null

          manual_number_i = manual_number_i + 1;

          manual_index_i = [manual_index_i; ix]; %#ok<AGROW>

       end

    end

    %outside contour

    full_ocontour_filename = get_contour_filename(manual_contour_path,para.name_prefix, sindex, para.outside_contour_mode,para.manual_seg_mode);

    if exist(full_ocontour_filename,'file') 

       fileinfo_o = dir(full_ocontour_filename);

       if fileinfo_o.bytes ~= 0 % not null

          manual_number_o = manual_number_o + 1;

          manual_index_o = [manual_index_o; ix]; %#ok<AGROW>

       end

    end

end

if (manual_number_i == 0)

   disp([manual_contour_path ' :NO inside contours!'])

   return; 

end

if (manual_number_o == 0)

   disp([manual_contour_path ' :NO outside contours!'])

   return; 

end

para.manual_index_i = manual_index_i;

para.manual_index_o = manual_index_o;

%count auto inside contours according to manual_index_i 

auto_number_i = 0;

auto_missing_index_i = [];

for ix = 1:length(manual_index_i)

    sindex = add_zero_index(manual_index_i(ix),para.digit_length);

    full_icontour_filename = get_contour_filename(auto_contour_path,para.name_prefix, sindex, para.inside_contour_mode, para.auto_seg_mode);

    if exist(full_icontour_filename,'file') 

       fileinfo_i = dir(full_icontour_filename);

       if fileinfo_i.bytes ~= 0 % not null

          auto_number_i = auto_number_i + 1;  

       else

          auto_missing_index_i = [auto_missing_index_i  manual_index_i(ix)]; %#ok<AGROW> 

       end  

    else

       auto_missing_index_i = [auto_missing_index_i  manual_index_i(ix)]; %#ok<AGROW>

    end

end

if (auto_number_i == 0)

   disp([auto_contour_path  ' :NO inside contours!'])

   return; 

end

%count auto outside contours according to manual_index_o 

auto_number_o = 0;

auto_missing_index_o = [];

for ix = 1:length(manual_index_o)

    sindex = add_zero_index(manual_index_o(ix),para.digit_length);

    full_ocontour_filename = get_contour_filename(auto_contour_path,para.name_prefix, sindex, para.outside_contour_mode,para.auto_seg_mode);

    if exist(full_ocontour_filename,'file') 

       fileinfo_o = dir(full_ocontour_filename);

       if fileinfo_o.bytes ~= 0 % not null

          auto_number_o = auto_number_o + 1;

       else

           auto_missing_index_o = [auto_missing_index_o  manual_index_o(ix)]; %#ok<AGROW>

       end

    else

       auto_missing_index_o = [auto_missing_index_o  manual_index_o(ix)]; %#ok<AGROW>

    end

end

if (auto_number_o == 0)

   disp([auto_contour_path ' :NO outside contours!'])

   return; 

end

%-calc perpendicular distance & Dice Metric

%inside contours

avg_dist_i=ones(length(manual_index_i),1)*para.init_value;

dm_i = ones(length(manual_index_i),1)*para.init_value;

auto_bad_index_i = [];

for ix = 1:length(manual_index_i)

    sindex = add_zero_index(manual_index_i(ix),para.digit_length);

    % manual contour

    manual_icontour_filename = get_contour_filename(manual_contour_path,para.name_prefix,sindex, para.inside_contour_mode,para.manual_seg_mode);

    % auto contour

    auto_icontour_filename = get_contour_filename(auto_contour_path,para.name_prefix,sindex, para.inside_contour_mode,para.auto_seg_mode);

    % if both exist, compare

    if (exist(manual_icontour_filename,'file') &&  exist(auto_icontour_filename,'file'))

         manual_xy =  load(manual_icontour_filename);

         auto_xy =  load(auto_icontour_filename);

         max_auto_x = max(auto_xy(:,1));

         max_auto_y = max(auto_xy(:,2));

         min_auto_xy = min(auto_xy(:));

         if (size(auto_xy,1) > 15  && min_auto_xy > 1 && max_auto_x < para.width && max_auto_y < para.height) %% if "point number of auto contour" < 1, may lead to difficulty to find a subset of reference that correspond to the current target

            try

                %perpendicular distance

                if (para.auto_based_noraml == 1) %calc normal based on auto contours, manual constour is the reference contour

                   avg_dist_i(ix) =calc_dist(auto_xy,manual_xy,auto_icontour_filename,para);

                else

                   avg_dist_i(ix) =calc_dist(manual_xy,auto_xy,auto_icontour_filename,para); 

                end

                %dice metric

                dm_i(ix) = calc_dm(auto_xy,manual_xy,para);

            catch

                s = lasterror;

                disp(s.message);

                continue;

            end

         end

     end

end

if max(avg_dist_i) == para.init_value

    return;

end

auto_bad_ix_i =find(avg_dist_i >= para.dist_limit);

if ~isempty(auto_bad_ix_i)

    auto_bad_index_i = manual_index_i(auto_bad_ix_i);

end

auto_good_ix_i =find(avg_dist_i < para.dist_limit & avg_dist_i ~= para.init_value);

auto_good_index_i = manual_index_i(auto_good_ix_i);

para.auto_good_index_i = auto_good_index_i ;

%outside contours

auto_bad_index_o = [];

avg_dist_o=ones(length(manual_index_o),1)*para.init_value;

dm_o = ones(length(manual_index_o),1)*para.init_value;

for ix = 1:length(manual_index_o)

    sindex = add_zero_index(manual_index_o(ix),para.digit_length);

    % manual contour

    manual_ocontour_filename = get_contour_filename(manual_contour_path,para.name_prefix,sindex, para.outside_contour_mode,para.manual_seg_mode);

    % auto contour

    auto_ocontour_filename = get_contour_filename(auto_contour_path,para.name_prefix,sindex, para.outside_contour_mode,para.auto_seg_mode);

    % if both exist, compare

    if (exist(manual_ocontour_filename,'file') &&  exist(auto_ocontour_filename,'file'))

         manual_xy =  load(manual_ocontour_filename);

         auto_xy =  load(auto_ocontour_filename);

         max_auto_x = max(auto_xy(:,1));

         max_auto_y = max(auto_xy(:,2));

         min_auto_xy = min(auto_xy(:));

         if (size(auto_xy,1) > 15  && min_auto_xy > 1 && max_auto_x < para.width && max_auto_y < para.height)  % if "point number of auto contour" < 1, may lead to difficulty to find a subset of reference that correspond to the current target

            try

                %perpendicular distance

                if (para.auto_based_noraml == 1) %calc normal based on auto contours, manual constour is the reference contour

                   avg_dist_o(ix) =calc_dist(auto_xy,manual_xy,auto_ocontour_filename,para);

                else

                   avg_dist_o(ix) =calc_dist(manual_xy,auto_xy,auto_ocontour_filename,para); 

                end

                %dice metric

                dm_o(ix) = calc_dm(auto_xy,manual_xy,para);

            catch

               s = lasterror;

               disp(s.message);

               continue;

            end

         end

     end

end

if max(avg_dist_o) == para.init_value

    return;

end

auto_bad_ix_o =find(avg_dist_o >= para.dist_limit);

if ~isempty(auto_bad_ix_o)

    auto_bad_index_o = manual_index_o(auto_bad_ix_o);

end

% if auto_good_ix_o is empty, mean(dm_o(auto_good_ix_o)) can be a NAN (undefined)

auto_good_ix_o = find(avg_dist_o < para.dist_limit & avg_dist_o ~= para.init_value);

auto_good_index_o = manual_index_o(auto_good_ix_o);

para.auto_good_index_o = auto_good_index_o;

%-calc ejection fraction and lv mass

para.auto = false;

lv_manual = calc_clinical_para(dicom_path, manual_contour_path, para);

para.auto = true;

lv_auto = calc_clinical_para(dicom_path, auto_contour_path, para,1,para.image_number/para.phase_number,lv_manual.es,lv_manual.ed);

%record result

compare_result.auto_number_i = auto_number_i; %total auto inside contours

compare_result.auto_number_o = auto_number_o;

compare_result.manual_number_i = manual_number_i;

compare_result.manual_number_o = manual_number_o;

compare_result.detect_percent_i = auto_number_i / manual_number_i * 100; %percent of detected auto inside contours

compare_result.detect_percent_o = auto_number_o / manual_number_o * 100;

compare_result.auto_missing_index_i = auto_missing_index_i; % missing auto inside contours

compare_result.auto_missing_index_o = auto_missing_index_o;

compare_result.auto_bad_index_i = auto_bad_index_i; %bad auto inside contours' index (average distance larger than para.dist_limit)

compare_result.auto_bad_index_o = auto_bad_index_o;

compare_result.good_percent_i = (auto_number_i - length(auto_bad_index_i)) /manual_number_i * 100; %percent of good auto inside contours

compare_result.good_percent_o = (auto_number_o - length(auto_bad_index_o)) /manual_number_o * 100;

compare_result.auto_ef_pic = lv_auto.ef_pic;%auto's ef, '_pic' means Papillary Included in the LV Cavity 

compare_result.auto_ef_pim = lv_auto.ef_pim; %'_pim' means Papillary Included in the Myocardium

compare_result.auto_lvm_pic = lv_auto.lvm_pic; %auto's lv mass

compare_result.auto_lvm_pim = lv_auto.lvm_pim; 

compare_result.manual_ef_pic = lv_manual.ef_pic; %manual's ef

compare_result.manual_ef_pim = lv_manual.ef_pim; 

compare_result.manual_lvm_pic = lv_manual.lvm_pic; %manual's lv mass

compare_result.manual_lvm_pim = lv_manual.lvm_pim; 

compare_result.avg_dist_i = mean(avg_dist_i(auto_good_ix_i)); %average distance

compare_result.avg_dist_o = mean(avg_dist_o(auto_good_ix_o));

compare_result.avg_dm_i = mean(dm_i(auto_good_ix_i)); %average dice metric

compare_result.avg_dm_o = mean(dm_o(auto_good_ix_o));

end %compare contour

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%sub functions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function sindex = add_zero_index(index,digit_length)

    %add zero before index with total length of digit_length, for instance, change 20 to 0020.

    sindex=int2str(index);

    while length(sindex) < digit_length

           sindex = ['0', sindex]; %#ok<AGROW>

    end

end 

function full_contour_filename = get_contour_filename(contour_path, name_prefix, sindex, contour_mode, seg_mode)

    %gef contour filename with full path

    full_contour_filename = [contour_path filesep name_prefix, sindex, '-',contour_mode,'-',seg_mode, '.txt'];

end

function varargout = calc_dist(target_xy,reference_xy,contour_filename,para)

%calc perpendicular distance between target and reference contours

    %plot contours

    if para.save_figure

        %create distance figure folder

        if (para.auto_based_noraml == 1)

            str = strrep(contour_filename, para.auto_contour_foldername, [para.distance_figure_foldername '_auto_based_normal']);

        else

            str = strrep(contour_filename, para.auto_contour_foldername, [para.distance_figure_foldername '_manual_based_normal']);

        end

        [pathstr, name, ext, versn] = fileparts(str);

        if (~exist(pathstr,'dir'))

            mkdir(pathstr)

        end

        %read dicom image

        str_dicom = strrep(contour_filename, [para.auto_contour_foldername filesep  para.auto_contour_subfoldername],para.dicom_foldername);

        dicom_filename = [str_dicom(1:end-18),'.dcm']; %IM-0001-XXXX-icontour-auto.txt

        raw_image = dicomread(dicom_filename);

        %show box

        min_showbox_x = round(min([target_xy(:,1); reference_xy(:,1)]));

        max_showbox_x = round(max([target_xy(:,1); reference_xy(:,1)]));

        min_showbox_y = round(min([target_xy(:,2); reference_xy(:,2)]));

        max_showbox_y = round(max([target_xy(:,2); reference_xy(:,2)]));

        %validate

        if min_showbox_x <1

           min_showbox_x = 1;

        end

        if min_showbox_y <1

           min_showbox_y = 1;

        end

        showbox_offset= 3;

        if ((min_showbox_x(1)-showbox_offset)<1 || (min_showbox_y(1)-showbox_offset)<1 || (max_showbox_x(1)+showbox_offset)>size(raw_image,2) || (max_showbox_y(1)+showbox_offset)>size(raw_image,1))

            showbox_offset = 0;

        end

        %crop image

        croped_image = raw_image(min_showbox_y(1)-showbox_offset:max_showbox_y(1)+showbox_offset,min_showbox_x(1)-showbox_offset:max_showbox_x(1)+showbox_offset);

        %show image

        h = figure;

        set(h,'Visible','off');

        axis image; clf;

        imshow(croped_image,[],'InitialMagnification',1000,'Border','tight')

        %plot contours

         if (para.auto_based_noraml == 1) %calc normal based on auto contours, manual constour is the reference contour

             hold on; plot(target_xy(:,1)-min_showbox_x(1)+showbox_offset,target_xy(:,2)-min_showbox_y(1)+showbox_offset,'r.-');

             hold on; plot(reference_xy(:,1)-min_showbox_x(1)+showbox_offset,reference_xy(:,2)-min_showbox_y(1)+showbox_offset,'g.-');

         else

             hold on; plot(reference_xy(:,1)-min_showbox_x(1)+showbox_offset,reference_xy(:,2)-min_showbox_y(1)+showbox_offset,'r.-');

             hold on; plot(target_xy(:,1)-min_showbox_x(1)+showbox_offset,target_xy(:,2)-min_showbox_y(1)+showbox_offset,'g.-');

         end

     end

    %check closed contours

    if target_xy(1,1) == target_xy(end,1) && target_xy(1,2) == target_xy(end,2)

        target_xy = target_xy(1:end-1,:);

    end

    if reference_xy(1,1) == reference_xy(end,1) && reference_xy(1,2) == reference_xy(end,2)

        reference_xy = reference_xy(1:end-1,:);

    end

    %for finding a subset of reference that correspond to the current target

    target_mask = poly2mask (target_xy(:,1),target_xy(:,2),double(para.width),double(para.height));

    reference_mask = poly2mask (reference_xy(:,1),reference_xy(:,2),double(para.width),double(para.height));

    %centroid of refernece mask

    reference_cen = regionprops(bwlabel(reference_mask),'Centroid');

    reference_cen_x = round(reference_cen.Centroid(1));

    reference_cen_y = round(reference_cen.Centroid(2));

    %check if target mask's centroid is in refernece's mask

    if  target_mask(reference_cen_y,reference_cen_x) ==1

        target_cen = reference_cen;

    else

        target_cen = regionprops(bwlabel(target_mask),'Centroid');

    end

    %calc degrees between contours points and centroid

    target_degree = atan2( target_xy(:,2)-target_cen.Centroid(2),target_xy(:,1)-target_cen.Centroid(1))*180/pi;

    reference_degree = atan2( reference_xy(:,2)-reference_cen.Centroid(2),reference_xy(:,1)-reference_cen.Centroid(1))*180/pi;

    %initialize

    distance=zeros(size(target_xy,1),1);

    for idx = 1: size(target_xy,1)

        %Get 3 points,[left current right]

        if idx == 1

            target_xy_3 = [target_xy(end,:); target_xy(idx:idx+1,:)];

        elseif (idx == size(target_xy,1))

            target_xy_3 = [target_xy(end-1:end,:); target_xy(1,:)];

        else

            target_xy_3 = target_xy(idx-1:idx+1,:);

        end

        %current target to reference distance

        target_curr_x = target_xy_3(2,1);

        target_curr_y = target_xy_3(2,2);

        %find a subset of reference that correspond to the current target

        if ( abs(target_degree(idx)) >= 150 )

            reference_xy_subset = reference_xy(find(abs(reference_degree)>=130),:);

        else

            reference_xy_subset = reference_xy(find(abs(reference_degree-target_degree(idx))<50),:);

        end

        %plot reference_xy_subset

        %hold on; plot(reference_xy_subset(:,1)-min_showbox_x(1)+showbox_offset, reference_xy_subset(:,2)-min_showbox_y(1)+showbox_offset ,'r*')

        %fit line

        target_X = target_xy_3(:,1);

        target_Y = target_xy_3(:,2);

        %3 points are the same

        if ( sum(abs(target_X - target_X(2))) + sum(abs(target_Y - target_Y(2))) == 0 )

            continue;

        end

        target_LINE = [target_X ones(size(target_X))] \ target_Y; %left division-least squares,  %y=LINE(1)*x+LINE(2)

        target_LINE_inv = [target_Y ones(size(target_Y))] \ target_X; %left division-least squares,  %x=LINE(1)*y+LINE(2)  Reverse the coordinates.

        %if the line is very close to parallel with the Y axis, and the residuals for the fit

        %improve by reversing the coordinates, then use the line found by

        %reversing the coordinates.

        if ( abs(target_LINE_inv(1))< 0.1 && sum(abs([target_Y ones(size(target_Y))] * target_LINE_inv - target_X)) < sum(abs([target_X ones(size(target_X))] * target_LINE - target_Y)))

            target_LINE = 1./(target_LINE_inv + 0.00001) ;

            %hold on; plot(target_X- min_showbox_x(1) + showbox_offset,target_Y - min_showbox_y(1) + showbox_offset,'Y*')

        end

        %Determine the line normal to the above tangent line (target_LINE).

        %Normal line: Ax+By+C=0;

        if abs(target_LINE(1)) < 0.01 %fit line parallel to X axis, normal parallel to Y axis

            A_N = 1;

            B_N = 0;

            C_N = -target_curr_x;

            reference_xy_subset_to_normal_distance = abs(reference_xy_subset(:,1)- target_curr_x); 

        elseif abs(target_LINE(1)) > 100 %fit line parallel to Y axis, normal parallel to X axis

            A_N = 0;

            B_N = 1;

            C_N = -target_curr_y;

            reference_xy_subset_to_normal_distance = abs(reference_xy_subset(:,2)- target_curr_y); 

        else

            A_N = -1/target_LINE(1);

            B_N = -1;

            C_N = target_curr_y - A_N * target_curr_x;

            reference_xy_subset_to_normal_distance = abs( (A_N*reference_xy_subset(:,1) + B_N*reference_xy_subset(:,2) + C_N))/sqrt(A_N^2 + B_N^2);%normal: y=N(1)*x+N(2);

        end

        %find reference point most close to normal

        [min_dist_subset, min_idx_subset] = min(reference_xy_subset_to_normal_distance);

        %if there are more than one min distance points, choose the one most close to current target

        min_idx_subset_1= find(reference_xy_subset_to_normal_distance <= 0.5); 

        if (length(min_idx_subset_1)>=2)

             min_dist_subset_to_target_curr_distance = sqrt((target_curr_x -  reference_xy_subset(min_idx_subset_1,1)).^2 + (target_curr_y - reference_xy_subset(min_idx_subset_1,2)).^2);

             [min_dist_subset, min_idx_subset_2] = min(min_dist_subset_to_target_curr_distance);

             reference_to_min_dist_subset_distance = sqrt((reference_xy(:,1) -  reference_xy_subset(min_idx_subset_1(min_idx_subset_2),1)).^2 + (reference_xy(:,2) - reference_xy_subset(min_idx_subset_1(min_idx_subset_2),2)).^2);

        else

             reference_to_min_dist_subset_distance = sqrt((reference_xy(:,1) -  reference_xy_subset(min_idx_subset,1)).^2 + (reference_xy(:,2) - reference_xy_subset(min_idx_subset,2)).^2);

        end

        min_idx_temp = find(reference_to_min_dist_subset_distance == 0);

        min_idx = min_idx_temp(1);

        %Get 3 points of reference,[left most-close-to-normal right]

        if min_idx == 1

            reference_xy_3 = [reference_xy(end,:); reference_xy(min_idx:min_idx+1,:)];

        elseif (min_idx == size(reference_xy,1))

            reference_xy_3 = [reference_xy(end-1:end,:); reference_xy(1,:)];

        else

            reference_xy_3 = reference_xy(min_idx-1:min_idx+1,:);

        end

        %fit reference line: Ax+By+C=0; 

        reference_X = reference_xy_3(:,1);

        reference_Y = reference_xy_3(:,2);

        reference_LINE = [reference_X ones(size(reference_X))] \ reference_Y; %left division-least squares,         %y=LINE(1)*x+LINE(2)

        reference_LINE_inv = [reference_Y ones(size(reference_Y))] \ reference_X; %left division-least squares,  %y=LINE(1)*x+LINE(2)

        %fit line parallel to Y axis,

        if ( abs(reference_LINE_inv(1))< 0.1 && sum(abs([reference_Y ones(size(reference_Y))] * reference_LINE_inv - reference_X)) < sum(abs([reference_X ones(size(reference_X))] * reference_LINE - reference_Y)))

            reference_LINE = 1./(reference_LINE_inv + 0.00001) ;

            %hold on; plot(target_X- min_showbox_x(1) + showbox_offset,target_Y - min_showbox_y(1) + showbox_offset,'Y*')

        end

        if abs(reference_LINE(1)) < 0.01 % reference line parallel to X axis

            A_R = 0;

            B_R = 1;

            C_R = -reference_Y(2);

        elseif abs(reference_LINE(1)) >100  % reference line parallel to Y axis

            A_R = 1;

            B_R = 0;

            C_R = -reference_X(2);

        else

            A_R = reference_LINE(1);

            B_R = -1;

            C_R = reference_Y(2) - A_R * reference_X(2);

        end

        %intersection of normal and reference line

        AA=[A_N B_N; A_R B_R];

        BB =[-C_N; -C_R];

        XX =AA\BB; %  AA is a square matrix, AA\BB is roughly the same as inv(AA)*BB

        %validate intersection point

        XX_reference_distance = sqrt((XX(1) -  reference_xy_subset(:,1)).^2 + (XX(2) - reference_xy_subset(:,2)).^2  );

        if min(XX_reference_distance)<1.5

            %distance

            distance(idx) = sqrt((target_curr_x - XX(1))^2 + (target_curr_y - XX(2))^2);

            if para.save_figure

                %plot intersection point

                hold on; plot(XX(1)-min_showbox_x(1)+showbox_offset,XX(2)-min_showbox_y(1)+showbox_offset,'m*') 

                %plot distance line

                hold on; plot([target_curr_x - min_showbox_x(1) + showbox_offset XX(1) - min_showbox_x(1) + showbox_offset], [target_curr_y - min_showbox_y(1) + showbox_offset XX(2) - min_showbox_y(1) + showbox_offset],'b'); 

            end

        else

            distance(idx) = -999;

        end

    end

    %in mm

    distance = distance * para.pixel_spacing(1);

    distance_effective = distance(distance>0);

    varargout{1} =mean(distance_effective); 

    if para.save_figure

        %plot minimum and maximum distance point

        min_dist_idx = find(distance == min(distance_effective));

        max_dist_idx = find(distance == max(distance_effective));

        hold on; plot(target_xy(min_dist_idx(1),1)-min_showbox_x(1)+showbox_offset,target_xy(min_dist_idx(1),2)-min_showbox_y(1)+showbox_offset,'yo')

        hold on; plot(target_xy(max_dist_idx(1),1)-min_showbox_x(1)+showbox_offset,target_xy(max_dist_idx(1),2)-min_showbox_y(1)+showbox_offset,'yd')

        %text minimum and maximum distance

        text(target_xy(min_dist_idx(1),1)-min_showbox_x(1)+showbox_offset+0.5,target_xy(min_dist_idx(1),2)-min_showbox_y(1)+showbox_offset+0.5,[num2str(distance(min_dist_idx(1))),'(mm)'],'color','g','FontSize',11) %+0.5, not polt overlap

        text(target_xy(max_dist_idx(1),1)-min_showbox_x(1)+showbox_offset+0.5,target_xy(max_dist_idx(1),2)-min_showbox_y(1)+showbox_offset+0.5,[num2str(distance(max_dist_idx(1))),'(mm)'],'color','g','FontSize',11)

        %text mean and std of distances

        text(2,2,['Mean:' num2str(mean(distance_effective)) '(mm)' ', Std:' num2str(std(distance_effective)) '(mm)'],'color','g','FontSize',11);

        %legend

        legend('target','reference','intersection','distance line')

        %save figure as png

        saveas(h,  [str(1:end-9), '.png']); %IM-0001-XXXX-icontour-auto.txt

        close(h)

    end

end

function dm = calc_dm(autoPoints,manualPoints,para)

    %calc dice metric

    auto_mask = poly2mask (autoPoints(:,1),autoPoints(:,2),double(para.width),double(para.height));

    manual_mask = poly2mask (manualPoints(:,1),manualPoints(:,2),double(para.width),double(para.height));

    auto_size = sum(auto_mask(:)>0);

    manual_size = sum(manual_mask(:)>0);

    intersect_size = sum((auto_mask(:) + manual_mask(:))==2);

    dm = 2 * intersect_size / (auto_size + manual_size);

end

function lv = calc_clinical_para(dicomPath, contour_path, para,varargin)

    %calc ejection fraction and lv mass

    [dicomCount, sliceCount] = dicom_counter(dicomPath,para.phase_number);

    contTable_i = zeros(sliceCount, para.phase_number);

    contTable_o = zeros(sliceCount, para.phase_number);

    contTable_p1 = zeros(sliceCount, para.phase_number);

    contTable_p2 = zeros(sliceCount, para.phase_number);

    if length(varargin) > 3

        startingSlice = varargin{1};

        endingSlice = varargin{2};

        systolePhase = varargin{3};

        diastolePhase = varargin{4};

    elseif length(varargin) > 1

        startingSlice = varargin{1};

        endingSlice = varargin{2};

        systolePhase = 0;

        diastolePhase = 0;

    else 

        startingSlice = 1;

        endingSlice = sliceCount;

        systolePhase = 0;

        diastolePhase = 0;    

    end    

    if startingSlice < 1 || startingSlice > sliceCount

        startingSlice = 1;

    end

    if endingSlice < 1 || endingSlice > sliceCount

        endingSlice = sliceCount;

    end

    if startingSlice > endingSlice

        endTemp = startingSlice;

        startingSlice = endingSlice;

        endingSlice = endTemp;

    end   

    %list all contours

    contour_list = dir([contour_path filesep '*contour*.txt']);

    for contourIdx = 1:length(contour_list)

        contourFile = contour_list(contourIdx).name;

        %inside contour

        if (strcmp(contourFile(14),'i'))

            %manual contours

            if para.auto == false

               good_icontour = true;

            end

            %auto contours

            if para.auto == true

               good_icontour = any(ismember(para.auto_good_index_i, str2num(contourFile(9:12)))); 

            end

            if ~good_icontour

                continue;

            end

            match = regexp((contourFile),'-','start');

            imNum = contourFile(match(2) + 1:match(3) - 1);

            currSlice = getSlice(str2double(imNum),para.phase_number);

            currPhase = getPhase(str2double(imNum),para.phase_number);

            full_contour_filename = [contour_path filesep contourFile];

           try

               xy = load(full_contour_filename);

               area = polyarea(xy(:,1),xy(:,2));

               area_mm = area* para.pixel_spacing(1)^2;

               vol_cm3 = area_mm * (para.thickness + para.gap) *(1/10)^3; %change to cm3

               contTable_i(currSlice, currPhase) = vol_cm3; 

           catch

               s = lasterror;

               disp(s.message);

               return;

            end

        end

        %outside contour

        if (strcmp(contourFile(14),'o') )

            if para.auto == false

                good_ocontour = true;

            end

            if para.auto == true

               good_ocontour = any(ismember(para.auto_good_index_o, str2num(contourFile(9:12)))); 

            end

            if ~good_ocontour

                continue;

            end

            match = regexp((contourFile),'-','start');

            imNum = contourFile(match(2) + 1:match(3) - 1);

            currSlice = getSlice(str2double(imNum),para.phase_number);

            currPhase = getPhase(str2double(imNum),para.phase_number);

            full_contour_filename = [contour_path filesep contourFile];

            try

               xy = load(full_contour_filename);

               area = polyarea(xy(:,1),xy(:,2));

               area_mm = area* para.pixel_spacing(1)^2;

               vol_cm3 = area_mm * (para.thickness + para.gap) *(1/10)^3;%change to cm3

               contTable_o(currSlice, currPhase) = vol_cm3; 

           catch

               s = lasterror;

               disp(s.message);

               return;

            end

        end

        %pap contour

        if (strcmpi(contourFile(14), 'p') )

            if para.auto == false

                good_pcontour = true;

            end

            if para.auto == true

                good_pcontour = any(ismember([ para.auto_good_index_i; para.auto_good_index_o] , str2num(contourFile(9:12))));

            end

            if ~good_pcontour

                continue;

             end

            match = regexp((contourFile),'-','start');

            imNum = contourFile(match(2) + 1:match(3) - 1);

            currSlice = getSlice(str2double(imNum),para.phase_number);

            currPhase = getPhase(str2double(imNum),para.phase_number);

            full_contour_filename = [contour_path filesep contourFile];

            try

               xy = load(full_contour_filename);

               area = polyarea(xy(:,1),xy(:,2));

               area_mm = area* para.pixel_spacing(1)^2;

               vol_cm3 = area_mm * (para.thickness + para.gap) *(1/10)^3;%change to cm3

               if ( strcmpi(contourFile(15), '1') )

                    contTable_p1(currSlice, currPhase) = vol_cm3; 

               end

               if ( strcmpi(contourFile(15), '2') )

                    contTable_p2(currSlice, currPhase) = vol_cm3; 

               end

            catch

               s = lasterror;

               disp(s.message);

               return;

            end

        end %end if

    end     

    maxPhase = length(contTable_i(1,:));

    if systolePhase < 1 || systolePhase > maxPhase || diastolePhase < 1 || diastolePhase > maxPhase

        [systolePhase, diastolePhase] = esedDetermine(contTable_i);

    end

    [contTableChecked_i, zeroed, slicesExcluded, slicesIncluded] = efCheck(contTable_i, systolePhase, diastolePhase, startingSlice, endingSlice);

    %calc ef - use only inside contour info.

    esvol_pic = sum(contTableChecked_i(:, systolePhase)); %'_pic' means pap included in the LV cavity

    esvol_pim = sum(contTableChecked_i(:, systolePhase))- sum(contTable_p1(slicesIncluded,systolePhase))- sum(contTable_p2(slicesIncluded,systolePhase));%'_pim' means pap included in the myocardium

    edvol_pic = sum(contTableChecked_i(:, diastolePhase));

    edvol_pim = sum(contTableChecked_i(:, diastolePhase)) - sum(contTable_p1(slicesIncluded,diastolePhase))- sum(contTable_p2(slicesIncluded,diastolePhase));

    strokeVol_pic = edvol_pic - esvol_pic;

    strokeVol_pim = edvol_pim - esvol_pim;

    ef_pic = strokeVol_pic / edvol_pic * 100;

    ef_pim = strokeVol_pim / edvol_pim * 100;

    lv.esv_pic = esvol_pic;

    lv.esv_pim = esvol_pim;

    lv.edv_pic = edvol_pic;

    lv.edv_pim = edvol_pim;

    lv.sv_pic = strokeVol_pic;

    lv.sv_pim = strokeVol_pim;

    lv.ef_pic = ef_pic;

    lv.ef_pim = ef_pim;

    lv.es = systolePhase;

    lv.ed = diastolePhase;

    lv.zeroed = zeroed;

    lv.slice_excluded = slicesExcluded;

    %calc lv mass 

    combined_table(:,1) = contTable_i(:,diastolePhase);

    combined_table(:,2) = contTable_o(:,diastolePhase);

    common_slice = ((combined_table(:,1) ~= 0) + (combined_table(:,2) ~= 0)) == 2 ;

    edv_i_pic = sum(combined_table(common_slice, 1));

    edv_o_pic = sum(combined_table(common_slice, 2));

    lvm_pic = (edv_o_pic - edv_i_pic) * 1.05; %1.05 (g/cm3)

    lvm_pim = (edv_o_pic - edv_i_pic + sum(contTable_p1(common_slice,diastolePhase)) + sum(contTable_p2(common_slice,diastolePhase)) ) * 1.05; 

    lv.lvm_pic = lvm_pic;

    lv.lvm_pim = lvm_pim;

end

function [systolePhase, diastolePhase] = esedDetermine(contTable)

    %determin ES and ED phase 

    contVolSum = sum(contTable);

    minVol = 9999; 

    maxVol = 0;

    for x = 1:length(contVolSum)

        valueUnderExamination = contVolSum(x);

        if valueUnderExamination > 0

            if minVol > valueUnderExamination

                minVol = contVolSum(x);

                systolePhase = x;

            end

            if maxVol < valueUnderExamination

                maxVol = contVolSum(x);

                diastolePhase = x;

            end

        end

    end

end

function sliceNum = getSlice(imageNum,phase_number)

    % calc slice number from the image number

     sliceNum = ceil(imageNum/phase_number);

end

function phaseNum = getPhase(imageNum,phase_number)

    %calc phase number from image number

    remainder = mod(imageNum,phase_number);

    if (remainder ~= 0)

        phaseNum = remainder;

    else

        phaseNum = phase_number;

    end 

end

function [volTableCropped, zeroed, slicesExcluded,slicesIncluded] = efCheck(volTable, systolePhase, diastolePhase, startSlice, endSlice)

    % modifies the contour table that is passed in by...

    % 1. removing the slices that are to be excluded using startSlice and

    % endSlice by zeroing all those values

    % 2. removing ED/ES pairs that are missing the other value.  

    % that is, if there is a zero in the systolic phase but not in the

    % diastolic phase in that slice, zero the other

    % 1. removing the non-included slices

    volTableCropped = zeros(length(volTable(:,1)), length(volTable(1,:)));

    volTableCropped(startSlice:endSlice, :) = volTable(startSlice:endSlice, :);

    % 2. finding the ES/ED slices with missing values

    sysZeros = find(volTable(:, systolePhase) == 0);

    diaZeros = find(volTable(:, diastolePhase) == 0);

    if ~isempty(sysZeros)

        volTableCropped(sysZeros, :) = 0;

    end

    if ~isempty(diaZeros)

        volTableCropped(diaZeros, :) = 0;

    end

    slicesExcluded = find(volTableCropped(:, systolePhase) == 0);

    slicesIncluded = find(volTableCropped(:, systolePhase) ~= 0);

    zeroed = length(slicesExcluded);

end

function [dicom_count, slice_count] = dicom_counter(dicom_path, phase_number)

    %calc dicom image number and slice number

    dicom_count = 0;

    slice_count = 0;

    if ~exist(dicom_path,'dir')

        return

    end

    dicom_files = dir([dicom_path filesep '*.dcm']);

    if isempty(dicom_files)

        disp([dicom_path ' :NO DICOM files!'])

        return;

    else 

        dicom_count = length(dicom_files);

        slice_count = ceil(dicom_count / phase_number);

    end

end