octave_packages/geometry-1.5.0/geom2d/crackPattern2.m

   1 %% Copyright (c) 2011, INRA
   2 %% 2004-2011, David Legland <david.legland@grignon.inra.fr>
   3 %% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   4 %%
   5 %% All rights reserved.
   6 %% (simplified BSD License)
   7 %%
   8 %% Redistribution and use in source and binary forms, with or without
   9 %% modification, are permitted provided that the following conditions are met:
  10 %%
  11 %% 1. Redistributions of source code must retain the above copyright notice, this
  12 %%    list of conditions and the following disclaimer.
  13 %%
  14 %% 2. Redistributions in binary form must reproduce the above copyright notice,
  15 %%    this list of conditions and the following disclaimer in the documentation
  16 %%    and/or other materials provided with the distribution.
  17 %%
  18 %% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19 %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20 %% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21 %% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22 %% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23 %% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24 %% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25 %% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26 %% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27 %% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28 %% POSSIBILITY OF SUCH DAMAGE.
  29 %%
  30 %% The views and conclusions contained in the software and documentation are
  31 %% those of the authors and should not be interpreted as representing official
  32 %% policies, either expressed or implied, of copyright holder.
  33
  34 %% -*- texinfo -*-
  35 %% @deftypefn {Function File} {@var{e} = } crackPattern2 (@var{box}, @var{points}, @var{alpha})
  36 %% Create a (bounded) crack pattern tessellation
  37 %%
  38 %%   E = crackPattern2(BOX, POINTS, ALPHA)
  39 %%   create a crack propagation pattern wit following parameters :
  40 %%   - pattern is bounded by area BOX which is a polygon.
  41 %%   - each crack originates from points given in POINTS
  42 %%   - directions of each crack is given by a [NxM] array ALPHA, where M is
  43 %%   the number of rays emanating from each seed/
  44 %%   - a crack stop when it reaches another already created crack.
  45 %%   - all cracks stop when they reach the border of the frame, given by box
  46 %%   (a serie of 4 points).
  47 %%   The result is a collection of edges, in the form [x1 y1 x2 y2].
  48 %%
  49 %%   E = crackPattern2(BOX, POINTS, ALPHA, SPEED)
  50 %%   Also specify speed of propagation of each crack.
  51 %%
  52 %%
  53 %%   See the result with :
  54 %%     figure;
  55 %%     drawEdge(E);
  56 %%
  57 %%   @seealso{drawEdge}
  58 %% @end deftypefn
  59
  60 function edges = crackPattern2(box, points, alpha, varargin)
  61
  62   if ~isempty(varargin)
  63       speed = varargin{1};
  64   else
  65       speed = ones(size(points, 1), 1);
  66   end
  67
  68   % Compute line equations for each initial crack.
  69   % The 'Inf' at the end correspond to the position of the limit.
  70   % If an intersection point is found with another line, but whose position
  71   % is after this value, this means that another crack stopped it before it
  72   % reach the intersection point.
  73   NP = size(points, 1);
  74   lines = zeros(0, 5);
  75   for i=1:size(alpha, 2)
  76       lines = [lines; points speed.*cos(alpha(:,i)) speed.*sin(alpha(:,i)) Inf*ones(NP, 1)];
  77   end
  78   NL = size(lines, 1);
  79
  80   % initialize lines for borders, but assign a very high speed, to be sure
  81   % borders will stop all cracks.
  82   dx = (box([2 3 4 1],1)-box([1 2 3 4],1))*max(speed)*5;
  83   dy = (box([2 3 4 1],2)-box([1 2 3 4],2))*max(speed)*5;
  84
  85   % add borders to the lines set
  86   lines = [lines ; createLine(box, dx, dy) Inf*ones(4,1)];
  87
  88   edges = zeros(0, 4);
  89
  90
  91   while true
  92       modif = 0;
  93
  94       % try to update each line
  95           for i=1:NL
  96
  97           % initialize first point of edge
  98           edges(i, 1:2) = lines(i, 1:2);
  99
 100           % compute intersections with all other lines
 101           pi = intersectLines(lines(i,:), lines);
 102
 103           % compute position of all intersection points on the current line
 104           pos = linePosition(pi, lines(i,:));
 105
 106
 107           % consider points to the right (positive position), and sort them
 108           indr = find(pos>1e-12 & pos~=Inf);
 109           [posr, indr2] = sort(pos(indr));
 110
 111
 112           % look for the closest intersection to the right
 113           for i2=1:length(indr2)
 114
 115               % index of intersected line
 116               il = indr(indr2(i2));
 117
 118               % position of point relative to intersected line
 119               pos2 = linePosition(pi(il, :), lines(il, :));
 120
 121               % depending on the sign of position, tests if the line2 can
 122               % stop the current line, or if it was stopped before
 123               if pos2>0
 124                   if pos2<abs(posr(i2)) && pos2<lines(il, 5)
 125                       if lines(i, 5) ~= posr(i2)
 126                           edges(i, 3:4) = pi(il,:);
 127                           lines(i, 5) = posr(i2);
 128                           modif = 1;
 129                       end
 130                       break;
 131                   end
 132               end
 133           end   % end processing of right points of the line
 134
 135
 136           end % end processing of all lines
 137
 138       % break the infinite loop if no more modification was made
 139       if ~modif
 140           break;
 141       end
 142   end
 143
 144   % add edges of the surronding box.
 145   edges = [edges ; box(1,:) box(2,:) ; box(2,:) box(3,:); ...
 146                    box(3,:) box(4,:) ; box(4,:) box(1,:)  ];
 147
 148 endfunction