octave_packages/geometry-1.5.0/geom2d/crackPattern.m

   1 %% Copyright (c) 2011, INRA
   2 %% 2007-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} = } crackPattern (@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 = crackPattern(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 two '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   % There is one 'end position' for each side of the crack.
  74   lines = [points speed.*cos(alpha) speed.*sin(alpha) Inf*ones(size(points, 1), 2)];
  75
  76   % initialize lines for borders, but assign a very high speed, to be sure
  77   % borders will stop all cracks.
  78   dx = (box([2 3 4 1],1)-box([1 2 3 4],1))*max(speed)*5;
  79   dy = (box([2 3 4 1],2)-box([1 2 3 4],2))*max(speed)*5;
  80
  81   % add borders to the lines set
  82   lines = [lines ; createLine(box, dx, dy) Inf*ones(4,2)];
  83
  84   edges = zeros(0, 4);
  85
  86
  87   while true
  88       modif = 0;
  89
  90       % try to update each line
  91           for i=1:size(points, 1)
  92
  93           % compute intersections with all other lines
  94           pi = intersectLines(lines(i,:), lines);
  95
  96           % compute position of all intersection points on the current line
  97           pos = linePosition(pi, lines(i,:));
  98
  99           % consider points to the right (positive position), and sort them
 100           indr = find(pos>=0 & pos~=Inf);
 101           [posr, indr2] = sort(pos(indr));
 102
 103
 104           % look for the closest intersection to the right
 105           for i2=1:length(indr2)
 106
 107               % index of intersected line
 108               il = indr(indr2(i2));
 109
 110               % position of point relative to intersected line
 111               pos2 = linePosition(pi(il, :), lines(il, :));
 112
 113               % depending on the sign of position, tests if the line2 can
 114               % stop the current line, or if it was stopped before
 115               if pos2>0
 116                   if pos2<abs(posr(i2)) && pos2<lines(il, 5)
 117                       if lines(i, 5) ~= posr(i2)
 118                           edges(i, 3:4) = pi(il,:);
 119                           lines(i, 5) = posr(i2);
 120                           modif = 1;
 121                       end
 122                       break;
 123                   end
 124               else
 125                    if abs(pos2)<abs(posr(i2)) && abs(pos2)<lines(il, 6)
 126                       if lines(i, 5) ~= posr(i2);
 127                           edges(i, 3:4) = pi(il,:);
 128                           lines(i, 5) = posr(i2);
 129                           modif = 1;
 130                       end
 131                       break;
 132                   end
 133               end
 134
 135           end   % end processing of right points of the line
 136
 137
 138
 139            % consider points to the left (negative position), and sort them
 140           indl = find(pos<=0 && pos~=Inf);
 141           [posl, indl2] = sort(abs(pos(indl)));
 142
 143           % look for the closest intersection to the right
 144           for i2=1:length(indl2)
 145               % index of intersected line
 146               il = indl(indl2(i2));
 147
 148               % position of point relative to intersected line
 149               pos2 = linePosition(pi(il, :), lines(il, :));
 150
 151               % depending on the sign of position, tests if the line2 can
 152               % stop the current line, or if it was stopped before
 153               if pos2>0
 154                   if pos2<abs(posl(i2)) && pos2<lines(il, 5)
 155                       if lines(i, 6) ~= abs(posl(i2));
 156                           edges(i, 1:2) = pi(il, :);
 157                           lines(i, 6) = abs(posl(i2));
 158                           modif = 1;
 159                       end
 160                       break;
 161                   end
 162               else
 163                    if abs(pos2)<abs(posl(i2)) && abs(pos2)<lines(il, 6)
 164                       if lines(i, 6) ~= abs(posl(i2));
 165                           edges(i, 1:2) = pi(il, :);
 166                           lines(i, 6) = abs(posl(i2));
 167                           modif = 1;
 168                       end
 169                       break;
 170                   end
 171               end
 172
 173           end % end processing of left points of the line
 174
 175
 176           end % end processing of all lines
 177
 178       % break the infinite loop if no more modification was made
 179       if ~modif
 180           break;
 181       end
 182   end
 183
 184   % add edges of the surronding box.
 185   edges = [edges ; box(1,:) box(2,:) ; box(2,:) box(3,:); ...
 186                    box(3,:) box(4,:) ; box(4,:) box(1,:)  ];
 187
 188 endfunction
 189