octave_packages/geometry-1.5.0/geom2d/beltproblem.m

   1 %% Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   2 %%
   3 %%    This program is free software: you can redistribute it and/or modify
   4 %%    it under the terms of the GNU General Public License as published by
   5 %%    the Free Software Foundation, either version 3 of the License, or
   6 %%    any later version.
   7 %%
   8 %%    This program is distributed in the hope that it will be useful,
   9 %%    but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 %%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 %%    GNU General Public License for more details.
  12 %%
  13 %%    You should have received a copy of the GNU General Public License
  14 %%    along with this program. If not, see <http://www.gnu.org/licenses/>.
  15
  16 %% -*- texinfo -*-
  17 %% @deftypefn {Function File} {[@var{tangent},@var{inner}] = } beltproblem (@var{c}, @var{r})
  18 %% Finds the four lines tangent to two circles with given centers and radii.
  19 %%
  20 %% The function solves the belt problem in 2D for circles with center @var{c} and
  21 %% radii @var{r}.
  22 %%
  23 %% @strong{INPUT}
  24 %% @table @var
  25 %% @item c
  26 %% 2-by-2 matrix containig coordinates of the centers of the circles; one row per circle.
  27 %% @item r
  28 %% 2-by-1 vector with the radii of the circles.
  29 %%@end table
  30 %%
  31 %% @strong{OUPUT}
  32 %% @table @var
  33 %% @item tangent
  34 %% 4-by-4 matrix with the points of tangency. Each row describes a segment(edge).
  35 %% @item inner
  36 %% 4-by-2 vector with the point of intersection of the inner tangents (crossed belts)
  37 %% with the segment that joins the centers of the two circles. If
  38 %% the i-th edge is not an inner tangent then @code{inner(i,:)=[NaN,NaN]}.
  39 %% @end table
  40 %%
  41 %% Example:
  42 %%
  43 %% @example
  44 %% c         = [0 0;1 3];
  45 %% r         = [1 0.5];
  46 %% [T inner] = beltproblem(c,r)
  47 %% @result{} T =
  48 %%  -0.68516   0.72839   1.34258   2.63581
  49 %%   0.98516   0.17161   0.50742   2.91419
  50 %%   0.98675  -0.16225   1.49338   2.91888
  51 %%  -0.88675   0.46225   0.55663   3.23112
  52 %%
  53 %% @result{} inner =
  54 %%   0.66667   2.00000
  55 %%   0.66667   2.00000
  56 %%       NaN       NaN
  57 %%       NaN       NaN
  58 %%
  59 %% @end example
  60 %%
  61 %% @seealso{edges2d}
  62 %% @end deftypefn
  63
  64 function [edgeTan inner] = beltproblem(c,r)
  65
  66   x0 = [c(1,1) c(1,2) c(2,1) c(2,2)];
  67   r0 = r([1 1 2 2]);
  68
  69   middleEdge = createEdge(c(1,:),c(2,:));
  70
  71   ind0 = [1 0 1 0; 0 1 1 0; 1 1 1 0; -1 0 1 0; 0 -1 1 0; -1 -1 1 0;...
  72           1 0 0 1; 0 1 0 1; 1 1 0 1; -1 0 0 1; 0 -1 0 1; -1 -1 0 1;...
  73           1 0 1 1; 0 1 1 1; 1 1 1 1; -1 0 1 1; 0 -1 1 1; -1 -1 1 1;...
  74           1 0 -1 0; 0 1 -1 0; 1 1 -1 0; -1 0 -1 0; 0 -1 -1 0; -1 -1 -1 0;...
  75           1 0 0 -1; 0 1 0 -1; 1 1 0 -1; -1 0 0 -1; 0 -1 0 -1; -1 -1 0 -1;...
  76           1 0 -1 -1; 0 1 -1 -1; 1 1 -1 -1; -1 0 -1 -1; 0 -1 -1 -1; -1 -1 -1 -1];
  77   nInit = size(ind0,1);
  78   ir = randperm(nInit);
  79   edgeTan = zeros(4,4);
  80   inner = zeros(4,2);
  81   nSol = 0;
  82   i=1;
  83
  84   %% Solve for 2 different lines
  85   while nSol<4 && i<nInit
  86       ind = find(ind0(ir(i),:)~=0);
  87       x = x0;
  88       x(ind)=x(ind)+r0(ind);
  89       [sol f0 nev]= fsolve(@(x)belt(x,c,r),x);
  90       if nev~=1
  91            perror('fsolve',nev)
  92       end
  93
  94       for j=1:4
  95          notequal(j) = all(abs(edgeTan(j,:)-sol) > 1e-6);
  96       end
  97       if all(notequal)
  98        nSol = nSol+1;
  99        edgeTan(nSol,:) = createEdge(sol(1:2),sol(3:4));
 100        % Find innerTangent
 101        inner(nSol,:) = intersectEdges(middleEdge,edgeTan(nSol,:));
 102       end
 103
 104       i = i+1;
 105   end
 106
 107 endfunction
 108
 109 function res = belt(x,c,r)
 110   res = zeros(4,1);
 111
 112   res(1,1) = (x(3:4) - c(2,1:2))*(x(3:4) - x(1:2))';
 113   res(2,1) = (x(1:2) - c(1,1:2))*(x(3:4) - x(1:2))';
 114   res(3,1) = sumsq(x(1:2) - c(1,1:2)) - r(1)^2;
 115   res(4,1) = sumsq(x(3:4) - c(2,1:2)) - r(2)^2;
 116
 117 end
 118
 119 %!demo
 120 %!  c         = [0 0;1 3];
 121 %!  r         = [1 0.5];
 122 %!  [T inner] = beltproblem(c,r)
 123 %!
 124 %!  figure(1)
 125 %!  clf
 126 %!  h = drawEdge(T);
 127 %!  set(h(find(~isnan(inner(:,1)))),'color','r');
 128 %!  set(h,'linewidth',2);
 129 %!  hold on
 130 %!  drawCircle([c(1,:) r(1); c(2,:) r(2)],'linewidth',2);
 131 %!  axis tight
 132 %!  axis equal
 133 %!
 134 %! % -------------------------------------------------------------------
 135 %! % The circles with the tangents edges are plotted. The solution with
 136 %! % crossed belts (inner tangets) is shown in red.