other/Matlab/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>0,step=0,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=squeeze(dir*(time_end-t(:,:,2,2))>0);
  22 for step=1:100,
  23     t_old=t;
  24     for i=1:m, for j=1:n,
  25         if active(i,j) & fire_mask(i,j),
  26             for a=1:3, for b=1:3, if a ~=2 | b~=2,
  27                 if print>1,fprintf('step %i point %i %i direction %i %i time %g ',step,i,j,a-2,b-2,t(i,j,a,b));end
  28                 dt = max(dir*(time_end-t(i,j,a,b)),0);    % time available to propagate
  29                 if dt>0,
  30                     dd = dt.*ros(i,j,a,b);        % distance traveled to tnow
  31                     if d(i,j,a,b)> dd,              % positive distance remains
  32                         t(i,j,a,b)=time_end;        % the end of the segment traveled is at time_now
  33                         d(i,j,a,b)= d(i,j,a,b)-dd;  % decrease the distances remaining
  34                         if print>1,fprintf('distance remaining %g time %g',d(i,j,a,b),t(i,j,a,b));end
  35                     elseif d(i,j,a,b)>0,
  36                         t_end=t(i,j,a,b)+dir*d(i,j,a,b)./ros(i,j,a,b); % time at the end point
  37                         if print>1,fprintf('time at end %g ',t_end);end
  38                         ii=i+a-2; % the grid point this end point coincides with
  39                         jj=j+b-2;
  40                         if ii>=1 & ii<=m,if jj>=1 & jj<=n & fire_area(ii,jj)
  41                             val=dir*min(dir*t(ii,jj,2,2),dir*t_end);
  42                             if print>1,fprintf('setting %i %i from %g to %g',ii,jj,t(ii,jj,2,2),val);end
  43                             t(ii,jj,:,:)=val;               % reinitialize the point ii, jj
  44                             d(ii,jj,:,:)=distance(ii,jj,:,:); % no distance traveled
  45                             active(ii,jj)=true;             % can propagate from this
  46                         end, end
  47                         t(i,j,a,b)=t_end;
  48                         d(i,j,a,b)=0;
  49                     end
  50                 end
  51                 if print>1,fprintf('\n');end
  52             end, end, end
  53         end
  54     end, end
  55     if print>0,step,tign=t(:,:,2,2),end
  56     done=~any(t(:)-t_old(:));
  57     if done,
  58         break
  59     end
  60 end
  61 tign=t(:,:,2,2);
  62