function varargout = calc_dist(target_xy,reference_xy,para)
%calc perpendicular distance between target and reference contours
%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(logical(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_mask=clean_segs(target_mask);
target_cen = regionprops(logical(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);
else
distance(idx) = -999;
end
end
%in mm
distance = distance * para.pixel_spacing(1);
distance_effective = distance(distance>0);
varargout{1} =mean(distance_effective);
end
