perimeter_new/propagate.m

   1 function [t,d]=propagate(t,d,dir,fire_area,fire_mask,distance,ros,time_end,print)
   2 % [t,d]=propagate(t,d,dir,fire_area,fire_mask,distance,ros,time_end,print)
   3 %   t(i,j,2,2)  fire arrival time location (i,j)
   4 %   t(i,j,a,b)  fire arrival time at a point on the line connecting (i,j) and (i+a-1,j+b-1)
   5 %   t(i,j,:,:) =fire arrival time given at (i,j) at start
   6 %   d(i,j,a,b)  distance remaining to (i+a-1,j+b-1), in proper units
   7 %   dir  - 1 time forward, -1 backward
   8 %   fire_area - 1 what can burn, 0 what not
   9 %   fire_mask - 1 can propagate from there, 0 -ignore
  10 %   distance(m,n,a,b) - distance on (i,j) and (i+a-1,j+b-1)
  11 %   ros(m,n,a,b) - rate of spread
  12 %   time_end - the latest time to propagate to (earliest if dir=-1)
  13 %   print  1 for tracing steps, 2 for detailed
  14
  15 if abs(dir) ~= 1, error('dir must be +-1'),end
  16 if print>1,,tign=t(:,:,2,2),end
  17 m=size(t,1);n=size(t,2);
  18 if ~exist('print','var'),
  19     print=0;
  20 end
  21 active=dir*(time_end-t(:,:,2,2))>0;
  22 for step=1:10*max(m,n),
  23     t_old=t;
  24     for i=[1:m,m-1:-1:1]
  25         for j=[1:n,n-1:-1:1]
  26             if active(i,j) & fire_mask(i,j),
  27                 for a=1:3,
  28                     for b=1:3,
  29                         if a ~=2 | b~=2,
  30                             if print>1,
  31                                 fprintf('step %i point %i %i direction %i %i time %g ',step,i,j,a-2,b-2,t(i,j,a,b));
  32                             end
  33                             dt = max(dir*(time_end-t(i,j,a,b)),0);    % time available to time)end
  34                             if dt>0,
  35                                 dd = dt.*ros(i,j,a,b);          % distance to travel until time_end
  36                                 if d(i,j,a,b)> dd,              % positive distance remains
  37                                     t(i,j,a,b)=time_end;        % the end of the segment traveled is at time_now
  38                                     d(i,j,a,b)= d(i,j,a,b)-dd;  % decrease the distances remaining
  39                                     if print>1,
  40                                         fprintf('distance remaining %g time %g',d(i,j,a,b),t(i,j,a,b));
  41                                     end
  42                                 elseif d(i,j,a,b)>0,
  43                                     t_end=t(i,j,a,b)+dir*d(i,j,a,b)./ros(i,j,a,b); % time at the segment end                                  if print>1,
  44                                     if print>1,
  45                                         fprintf('time at end %g ',t_end);
  46                                     end
  47                                     ii=i+a-2; % the grid point this end point coincides with
  48                                     jj=j+b-2;
  49                                     if ii>=1 & ii<=m & jj>=1 & jj<=n & fire_area(ii,jj)
  50                                         if dir>0,
  51                                             val=min(t(ii,jj,2,2),t_end);
  52                                         else
  53                                             val=max(t(ii,jj,2,2),t_end);
  54                                         end
  55                                         if print>1,
  56                                             fprintf('setting %i %i from %g to %g',ii,jj,t(ii,jj,2,2),val);
  57                                         end
  58                                         t(ii,jj,:,:)=val;               % reinitialize the point ii, jj
  59                                         d(ii,jj,:,:)=distance(ii,jj,:,:); % no distance traveled
  60                                         active(ii,jj)=true;             % can propagate from this
  61                                     end
  62                                     t(i,j,a,b)=t_end;
  63                                     d(i,j,a,b)=0;
  64                                 else
  65                                 end
  66                             end
  67                             if print>1,
  68                                 fprintf('\n');
  69                             end
  70                         end
  71                     end
  72                 end
  73             end
  74         end
  75     end
  76     change=norm(t(:)-t_old(:),1);
  77     if print>0,
  78         maxt=max(t(:));
  79         mint=min(t(:));
  80         fprintf('step %i tign min %g max %g change %g\n',step,mint,maxt,change)
  81     end
  82     if print>1,tign=t(:,:,2,2),end
  83     if change<eps,
  84         break
  85     end
  86 end
  87 end