cycling/squish2.m

   1 function [tign,tign_large] = squish2(ps)
   2 %function takes path structure and makes new estimated tign
   3 %inputs  ps = graph_dets(w) or ps = cluster_paths(w), struct with points, paths, etc.
   4
   5 %copy variables to shorter names
   6 pts = ps.grid_pts;
   7 pts(:,3) = ps.points(:,3);
   8 new_pts = pts;
   9 tign = ps.red.tign;
  10 t_max = max(ps.points(:,3));
  11 tign(tign>=t_max)=t_max;
  12 idx = ps.idx;
  13 new_idx = idx;
  14 fignum = 111;
  15
  16 %make new points by interpolation along the paths
  17 for i = 1:length(ps.paths)
  18     p = ps.paths(i).p;
  19     for j = 1:length(p)
  20 %         p_i = idx(p(j),1);
  21 %         p_j = idx(p(j),2);
  22         % new point
  23         if j > 1
  24             w1 = 0.5;
  25             w2 = 0.5;
  26 %             frp1 = ps.points(p(j-1),5);
  27 %             frp2 = ps.points(p(j),5);
  28 %             w1 = frp1/(frp1+frp2);
  29 %             w2 = frp2/(frp1+frp2);
  30             new_lat = w1*pts(p(j-1),1)+w2*pts(p(j),1);
  31             new_lon = w1*pts(p(j-1),2)+w2*pts(p(j),2);
  32             new_t = w1*pts(p(j-1),3)+w2*pts(p(j),3);
  33             [new_i,new_j,new_lat,new_lon]= fixpt(ps.red,[new_lat,new_lon]);
  34             new_pts=[new_pts;[new_lat new_lon new_t]];
  35             new_idx = [new_idx;[new_i,new_j]];
  36         end
  37     end % for j
  38 end %for i
  39 %new points ready to work with
  40 [new_pts,u1,u2]=unique(new_pts,'rows');
  41 % will this uniques command give good results used separately ???
  42 idx = new_idx(u1,1:2);
  43 times = new_pts(:,3);
  44 lats = new_pts(:,1);
  45 lons = new_pts(:,2);
  46
  47 %sort data by time
  48 [st1,st2] = sort(times);
  49 times = times(st2);
  50 lats = lats(st2);
  51 lons=lons(st2);
  52 idx(:,1:2) = idx(st2,1:2);
  53
  54 n = length(new_pts);
  55 end_time = max(max(times),ps.red.end_datenum);
  56 start_time = min(times);
  57 total_time = ceil(end_time-start_time);
  58
  59 tign = end_time*ones(size(ps.red.fxlong));
  60 temp_tign = tign;
  61 pts_tign = tign;
  62 lons_set = [];
  63 lats_set = [];
  64 times_set = [];
  65
  66 %make random patch around detections
  67 rm = 0;
  68 for i = 1:n;% total_time %half day interval
  69     if rm ~= 0
  70         pts_tign(idx(i,1)-round(rm*rand):idx(i,1)+round(rm*rand),idx(i,2)-round(rm*rand):idx(i,2)+round(rm*rand))=times(i);
  71     end
  72     %draw polygon around every 20 detections
  73     if mod(i,20) == 0 || i == n
  74     pt_set = times <= times(i);
  75     lons_set = [lons_set;lons(pt_set)];
  76     lats_set = [lats_set;lats(pt_set)];
  77     times_set = [times_set;times(pt_set)];
  78 %     figure,scatter3(lons_set,lats_set,times_set)
  79     %convex hull operations
  80 %     dt = delaunayTriangulation(lats_set,lons_set);
  81 %     c = convexHull(dt);
  82 %     [in,on] = inpolygon(ps.red.fxlat,ps.red.fxlong,lats_set(c),lons_set(c));
  83 %     figure,plot(lons_set(c),lats_set(c));
  84     [in,on] = inpolygon(ps.red.fxlat,ps.red.fxlong,lats_set,lons_set);
  85     %in = logical(in+on);
  86     temp_tign(in) = max(times_set);
  87     temp_tign(~in) = end_time;
  88     tign = min(tign,temp_tign);
  89     %tign = imgaussfilt(tign,1);
  90     figure(fignum)
  91     mesh(ps.red.fxlong,ps.red.fxlat,tign)
  92     end
  93 end
  94 %tign=min(tign,pts_tign);
  95 tign = imgaussfilt(tign,1/13,'FilterDomain','frequency');
  96 fprintf('Blurring ... \n')
  97 figure(fignum),mesh(ps.red.fxlong,ps.red.fxlat,tign)
  98 xlabel('Lon'),ylabel('Lat'),zlabel('Time'),title('Interpolated by Polygons')
  99 %resize, if small matrices were used
 100 if isfield(ps.red,'red')
 101     if size(ps.red.fxlong) ~= size(ps.red.red.fxlong)
 102         fprintf('resizing to original size of grid \n');
 103         [n,m] = size(ps.red.red.tign);
 104         la = linspace(ps.red.red.min_lat,ps.red.red.max_lat,m);
 105         lo = linspace(ps.red.red.min_lon,ps.red.red.max_lon,n);
 106         [lat,lon] = meshgrid(la,lo);
 107         F = scatteredInterpolant(ps.red.fxlat(:),ps.red.fxlong(:),tign(:));
 108         tign_large = F(lat,lon);
 109         %griddata seems to work the same as the scattererd interpolant
 110         %here...
 111         %tign = griddata(ps.red.fxlat(:),ps.red.fxlong(:),tign(:),lat,lon);
 112     end
 113 end
 114
 115
 116 end