How to find distance between two connected points through white pixels?

I want to find the two endpoints of a curve.

I used convex hull to find 8 boundary-points and then for every point, distances with all other points are calculated. The points with maximum distances are assumed to be the endpoints. It works for most of the cases when the shape is simple, like below:

I'm using this code to find the two points:

% Calculate 8-boundary points using Convex Hull
% (points consists of all white pixels for the curve)
out=convhull(points);
extremes=points(out,:);

% Finding two farthest points among above boundary points
arr_size = size(extremes);
max_distance = 0;
endpoints = [extremes(1,:); extremes(2,:)];
for i=1:arr_size(1)
    p1 = extremes(i,:);
    for j=i+1:arr_size(1)
        p2 = extremes(j,:);
        dist = sqrt((p2(1)-p1(1))^2 + (p2(2)-p1(2))^2);
        if dist>max_distance
            max_distance = dist;
            endpoints = [p1; p2];
        end
    end
end

disp(endpoints);

In most of the cases it works, but the issue comes when it is applied on U-typed shape.

I want to find out the number of white pixels between the two points. I followed this Counting white pixels between 2 points in an image but it helps in case of only straight lines not for curve. Any help appreciated!

EDIT 1 :

I have updated my code based on the answer by @Rotem and it fixed my problem.

Updated Code:

for i=1:arr_size(1)
    p1 = extremes(i,:);
    distanceMap = calculate_distance_map(arr, p1);
    for j=i+1:arr_size(1)
        p2 = extremes(j,:);
        dist = distanceMap(p2(2), p2(1));
        if dist>max_distance
            max_distance = dist;
            endpoints = [p1; p2];
        end
    end
end

However it takes significantly much time. Any suggestions on how to reduce the calculation time?

I found a solution inspired by dynamic programming algorithms.

Explanations are in the comments:

function n_points = BlueDotsDist()

I = imread('BlueDots.png');

%Fix the image uploaded to the site - each dot will be a pixel.
I = imresize(imclose(I, ones(3)), 1/6, 'nearest');

%Convert each color channel to a bynary image:
R = imbinarize(I(:,:,1));G = imbinarize(I(:,:,2));B = imbinarize(I(:,:,3));
%figure;imshow(cat(3, double(R), double(G), double(B)));impixelinfo

%Find blue dost:
[blue_y, blue_x] = find((B == 1) & (R == 0));

%Compute distance map from first blue point
P = CalcDistMap(R, blue_y(1), blue_x(1));

%Compute distance map from second blue point
Q = CalcDistMap(R, blue_y(2), blue_x(2));

%Get 3x3 pixels around second blue point.
A = P(blue_y(2)-1:blue_y(2)+1, blue_x(2)-1:blue_x(2)+1);
dist_p = min(A(:)); %Minimum value is the shortest distance from first point (but not the number of white points).

T = max(P, Q); %Each element of T is maximum distance from both points.
T(T > dist_p) = 10000; %Remove points that are more far than dist_p.

%Return number of white points between blue points.
n_points = sum(T(:) < 10000);




function P = CalcDistMap(R, blue_y, blue_x)
%Each white pixel in R gets the distance from the blue dot in coordinate (blue_x, blue_y).

%Initialize array with values of 10000 (high value - marks infinite distance).
P = zeros(size(R)) + 10000; %P - distance from blue dot

P(blue_y, blue_x) = 0; %Distance from itself.

prvPsum = 0;
while (sum(P(:) < 10000) ~= prvPsum)
    prvPsum = sum(P(:) < 10000); %Sum of "marked" dots.
    for y = 2:size(R, 1)-1
        for x = 2:size(R, 2)-1
            p = P(y, x);
            if (p < 10000) %If P(y,x) is "marked"
                A = P(y-1:y+1, x-1:x+1); %3x3 neighbors.
                A = min(A, p+1); %Distance is the minimum of current distance, and distance from p pixel + 1.
                A(R(y-1:y+1, x-1:x+1) == 0) = 10000; %Ignore black pixels.
                A(P(y-1:y+1, x-1:x+1) == 0) = 0; %Restore 0 in blue point.
                P(y-1:y+1, x-1:x+1) = A; %Update distance map.
            end
        end
    end
end

---

Illustration of distance maps (10000 replaced with zero):

P (distances from first blue point):

Q (distances from second blue point):

max(P, Q):

T: