3
\$\begingroup\$

Requirement

After a user draws two circles with their mouse successively, the program should draw all common tangents, if any, to them.

Here is my implementation

function InteractiveCommonTangent
tolerance = 0.001;
color = 'b';
style = '--';
width = 0.5;

circleCount = 1;
buttonDown = 0;
x1 = [0 0]; y1 = [0 0]; x2 = [0 0]; y2 = [0 0];

figure('WindowButtonDownFcn', @getBeginPoint, ...
    'WindowButtonMotionFcn', @updateCircle, ...
    'WindowButtonUpFcn', @getEndPoint);
ah = axes('SortMethod', 'childorder');
circles = [rectangle('Position', [0 0 0 0], 'Curvature', 1), ...
    rectangle('Position', [0 0 0 0], 'Curvature', 1)];
hold on;
axis equal;
grid on;
axis ([0 1 0 1]);

    function getBeginPoint(src, ~)
        if strcmp(get(src, 'SelectionType'), 'normal')
            buttonDown = 1;
            [x1(circleCount), y1(circleCount)] = get_point(ah);
            if circleCount == 1
                assets = findobj('Type', 'Line', ...
                    '-or', 'Type', 'Transform', ...
                    '-or', 'Type', 'Text');
                delete(assets);
                set(circles, 'Position', [0 0 0 0]);
            end
        end
    end

    function updateCircle(~, ~)
        if buttonDown
            [x, y] = get_point(ah);
            x0 = x1(circleCount);
            y0 = y1(circleCount);
            xx = (x + x0) / 2;
            yy = (y + y0) / 2;
            r = norm([x-x0, y-y0]) / 2;
            set(circles(circleCount), 'Position', [xx-r yy-r 2*r 2*r]);
            axis ([0 1 0 1]);
        end
    end

    function getEndPoint(~, ~)
        buttonDown = 0;
        [x2(circleCount), y2(circleCount)] = get_point(ah);
        circleCount = circleCount + 1;
        if circleCount > 2
            rawData = [x1' y1' x2' y2'];
            drawCommonTangent(rawData);
            axis ([0 1 0 1]);
            circleCount = 1;
        end
    end

    function [x, y] = get_point(ah)
        cp = get(ah, 'CurrentPoint');
        x = cp(1,1);
        y = cp(1,2);
    end

    function drawCommonTangent(rawCircles)
        r1 = norm(rawCircles(1, 3:4) - rawCircles(1, 1:2)) / 2;
        r2 = norm(rawCircles(2, 3:4) - rawCircles(2, 1:2)) / 2;
        if (r1 <= tolerance || r2 <= tolerance)
            set(circles, 'Position', [0 0 0 0]);
            return;
        end
        % make r1 >= r2
        if (r1 < r2)
            rawCircles = flip(rawCircles);
            [r2, r1] = deal(r1, r2);
        end

        xx1 = (rawCircles(1,1) + rawCircles(1,3)) / 2;
        yy1 = (rawCircles(1,2) + rawCircles(1,4)) / 2;
        xx2 = (rawCircles(2,1) + rawCircles(2,3)) / 2;
        yy2 = (rawCircles(2,2) + rawCircles(2,4)) / 2;
        vX = [xx1 xx2];
        vY = [yy1 yy2];
        d = norm([xx2-xx1, yy2-yy1]);
        unitTangent = [xx2-xx1, yy2-yy1] / d;
        unitNormal = [yy1-yy2, xx2-xx1] / d;

        set(circles(1), 'Position', [xx1-r1 yy1-r1 2*r1 2*r1]);
        set(circles(2), 'Position', [xx2-r2 yy2-r2 2*r2 2*r2]);

        if (abs(r1-r2) <= tolerance ...
                && abs(xx1-xx2) <= tolerance ...
                && abs(yy1-yy2) <= tolerance)
            % I didn't have much education. Don't try to fool me.
            text(xx1, yy1, 'THIS MAKES NO SENSE!', ...
                'HorizontalAlignment', 'center');

        else
            % Internal Common Tangents
            if (d + tolerance >= r1 + r2)
                if (d - tolerance <= r1 + r2)
                    center = deal([xx1 yy1] + unitTangent * r1);
                    theta = pi/2;
                    makeTransformedLine(vX, vY, center, theta, 3 * r1);
                else
                    D = min(realmax, d * r1/(r1+r2));
                    center = deal([xx1 yy1] + unitTangent * D);
                    theta = asin((r1+r2)/d);
                    makeTransformedLine(vX, vY, center, theta, 3 * D);
                    makeTransformedLine(vX, vY, center, -theta, 3 * D);
                end
            end

            % External Common Tangents
            if (d + tolerance >= r1 - r2)
                if (r1 - r2 <= tolerance)
                    [X1, Y1] = deal([xx1 yy1] - unitTangent * d);
                    [X2, Y2] = deal([xx2 yy2] + unitTangent * d);
                    delta = unitNormal * r1;
                    line([X1, X2] + delta(1), [Y1, Y2] + delta(2), ...
                        'Color', color, ...
                        'LineStyle', style, ...
                        'LineWidth', width);
                    line([X1, X2] - delta(1), [Y1, Y2] - delta(2), ...
                        'Color', color, ...
                        'LineStyle', style, ...
                        'LineWidth', width);
                elseif (d - tolerance <= r1 - r2)
                    center = deal([xx1 yy1] + unitTangent * r1);
                    theta = pi/2;
                    makeTransformedLine(vX, vY, center, theta, 3 * r1);
                else
                    D = min(realmax, d * r1/(r1-r2));
                    center = deal([xx1 yy1] + unitTangent * D);
                    theta = asin((r1-r2)/d);
                    makeTransformedLine(vX, vY, center, theta, 3 * D);
                    makeTransformedLine(vX, vY, center, -theta, 3 * D);
                end
            end
        end
    end

    function makeTransformedLine(vX, vY, center, theta, length)
        oldLength = norm([vX(2)-vX(1), vY(2)-vY(1)]);
        scale = min(realmax, length / oldLength);
        lineCenter = [mean(vX), mean(vY)];

        ht = hgtransform;
        line(vX, vY, 'Parent', ht, ...
            'Color', color, ...
            'LineStyle', style, ...
            'LineWidth', width);

        N = makehgtform('translate', -[lineCenter 0]);
        R = makehgtform('zrotate', theta);
        S = makehgtform('scale', scale);
        T = makehgtform('translate', [center 0]);

        set(ht, 'Matrix', T*S*R*N);
    end
end

The result may look like this

enter image description here

I'm not very familiar with matrix manipulation in MATLAB, which is said to be very handy. As a result, this piece code contain many inelegant statements like rawData = [x1' y1' x2' y2'];. Please help me improve them.

Also, I’ve put some effort in making this program more robust, but if you find any case which isn’t handled well, please point it out.

\$\endgroup\$
1
\$\begingroup\$

This is a nice little program. But it does have a bug, stemming from your misunderstanding of the deal function. This one is difficult to trigger, though. I'm not surprised you didn't know about it. I'll deal with this at the end.

You make nice use of the MATLAB graphics capabilities, especially hgtransform, which I have not seen used before.

You display pretty good general programming practices, in terms of program structure, function naming, etc.

The one aspect where this program could be improved is the separation of x and y components for the coordinates. This separation leads to several problems in the code:

  • Variable names

    They are in general very good, except when you're naming coordinates. You start with x1, y1, x2 and y2 for the two coordinates that define each circle. Then you use x, xx1, X1, vX, ... and things get a little confusing to me there. I think this stems from the division of x and y components into separate variables. It is the x and y component separation that drives the choice of names, rather than the function of the particular coordinate. If you had kept x and y together in a single variable, you would have been forced to come up with a better naming scheme.

    For example, use circle_edge1 and circle_edge2 for the two points defining each circle, with circle_edge1(1,:) the coordinates for the first point of circle 1, etc. Then xx1, yy1, xx2 and yy2 could be centers, with centers(1,:) the coordinates for the center of circle 1, etc.

  • Code repetition

    In many places you repeat calculations for the x and y component. Using a single variable to hold both components would mean you only write the calculation once (MATLAB automatically "vectorizes" many calculations, make use of that!).

  • Length of drawCommonTangent

    This function is a little too long to easily follow where xx1 and friends go. But with less code repetition this would be solved as well, I think.

As an example for simplification of code as a result of using the more natural coordinate representation, this bit of code from updateCircle:

[x, y] = get_point(ah);
x0 = x1(circleCount);
y0 = y1(circleCount);
xx = (x + x0) / 2;
yy = (y + y0) / 2;
r = norm([x-x0, y-y0]) / 2;
set(circles(circleCount), 'Position', [xx-r yy-r 2*r 2*r]);

could instead be written as

pt2 = get_point(ah); % assuming it returns a single value
pt1 = circle_edge1(circleCount);
center = (pt1 + pt2) / 2;
d = norm(pt1 - pt2) / 2;
set(circles(circleCount), 'Position', [center-d/2, d, d]);

You could also write

r1 = norm(rawCircles(1, 3:4) - rawCircles(1, 1:2)) / 2;
r2 = norm(rawCircles(2, 3:4) - rawCircles(2, 1:2)) / 2;

as

d = rawCircles(:, 3:4) - rawCircles(:, 1:2);
% r = norm(d) / 2;               % doesn't work, norm does something different on matrices
r = sqrt(sum(d.^2, 2)) / 2;      % alternative 1
r = hypot(d(:, 1), 2(:, 2)) / 2; % alternative 2

and then use r(1) and r(2).

Deal

The deal function in MATLAB is not commonly used. You have one place where it's used correctly:

[r2, r1] = deal(r1, r2);

Another place where you could have used it:

x1 = [0 0]; y1 = [0 0]; x2 = [0 0]; y2 = [0 0];

can be written as:

[x1, y1, x2, y2] = deal([0, 0]);

Several places where it's superfluous:

center = deal([xx1 yy1] + unitTangent * r1);

Here you're "dealing" one value to one variable, instead just directly assign it.

And one place where it causes a bug:

[X1, Y1] = deal([xx1 yy1] - unitTangent * d);
[X2, Y2] = deal([xx2 yy2] + unitTangent * d);

Here you are calculating coordinates for two points, and assigning these coordinates to two variables each. That is, after these two assignments, X1 == Y1 and X2 == Y2. This bug is triggered only when the radii of the two circles are within tolerance of each other, so I had to be patient to trigger it. When I did manage to trigger it, the internal tangent lines were not drawn in the correct places.

In short, [x,y,z]=deal(a,b,c) is the same as x=a; y=b; z=c, and [x,y,z]=deal(a) is the same as x=a; y=a; z=a. This is true no matter what the shape or type of a is. You seem to assume that [x,y,z]=deal([a,b,c]) is the same as the first case, but it is actually the second case, because [a,b,c] is a single value (deal gets a single input argument). This is probably very different from how it works in Python.

Besides being a nice shortcut for swapping the values of two variables, use cases for deal are few and far between.

\$\endgroup\$
  • \$\begingroup\$ Thanks for your detailed answer! The r = norm(rawCircles(:, 3:4) - ... suggestion doesn’t seem to work though. In my version (2017a), I found that norm always returns a scalar value, even when a 2D vector is passed to it. \$\endgroup\$ – nalzok May 1 '18 at 16:03
  • 1
    \$\begingroup\$ @nalzok: You're right. Sorry, should have paid more attention there. I have edited the post with two alternatives. It's not such a good example any more, for reducing amount of code. :) \$\endgroup\$ – Cris Luengo May 1 '18 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.