Skip to content

Commit 00f843f

Browse files
author
Nick Holschuh
committed
Functions generated during the 2023-24 Antarctic Field Season
1 parent fff971a commit 00f843f

19 files changed

+481
-63
lines changed

‎Untitled.ipynb

100644100755
File mode changed.

‎calculate_flowlines.py

100644100755
+203-33
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,22 @@
11
import numpy as np
22
import matplotlib.pyplot as plt
3+
import tqdm
4+
import xarray as xr
35

4-
def calculate_flowlines(x,y,u,v,seed_points,max_error=0.00001):
6+
################ This is the import statement required to reference scripts within the package
7+
import os,sys,glob
8+
ndh_tools_path_opts = [
9+
'/mnt/data01/Code/',
10+
'/mnt/l/mnt/data01/Code/',
11+
'/home/common/HolschuhLab/Code/'
12+
]
13+
for i in ndh_tools_path_opts:
14+
if os.path.isfile(i): sys.path.append(i)
15+
################################################################################################
16+
17+
import NDH_Tools as ndh
18+
19+
def calculate_flowlines(input_xr,seed_points,uv_varnames=['u','v'],xy_varnames=['x','y'],steps=20000,ds=2,forward0_both1_backward2=1):
520
"""
621
% (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
722
%
@@ -10,7 +25,7 @@ def calculate_flowlines(x,y,u,v,seed_points,max_error=0.00001):
1025
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1126
% The inputs are:
1227
%
13-
% input_array -- array of data to analyze
28+
% input_xr -- xarray dataarray that has the gradient objects in it
1429
%
1530
%%%%%%%%%%%%%%%
1631
% The outputs are:
@@ -20,34 +35,189 @@ def calculate_flowlines(x,y,u,v,seed_points,max_error=0.00001):
2035
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2136
"""
2237

23-
################# This uses a modified plt.streamline to pass through a user-editable keyword
24-
################# argument "max_error", which goes into the interpolater to guarantee
25-
################# accurate streammline calculation. Copy this version of streamplot
26-
################# into your matplotlib directory to enable the use of streamline
27-
################# calculation from NDH_Tools
28-
29-
30-
sls = []
31-
32-
if isinstance(seed_points,list):
33-
seed_points = np.array(seed_points)
34-
35-
if len(seed_points.shape) == 1:
36-
seed_points = np.expand_dims(seed_points,axis=0)
37-
38-
fig = plt.figure()
39-
40-
for ind0, sp in enumerate(seed_points[:,0]):
41-
streamlines = plt.streamplot(x,y,u,v,start_points=[seed_points[ind0,:]], max_error=max_error)
42-
43-
44-
########### Here we extract the coordinate information along the streamline
45-
sl = [streamlines.lines.get_paths()[0].vertices[0]]
46-
for i in streamlines.lines.get_paths():
47-
sl.append(i.vertices[1])
48-
49-
sls.append(np.array(sl))
50-
51-
plt.close(fig)
52-
53-
return sls
38+
##################### Here, we standardize the naming convention within the xarray object
39+
input_xr = input_xr.rename({xy_varnames[0]:'x',xy_varnames[1]:'y'})
40+
uv_scalar = np.sqrt(input_xr[uv_varnames[0]].values**2 + input_xr[uv_varnames[1]].values**2)
41+
input_xr[uv_varnames[0]] = (('y','x'),input_xr[uv_varnames[0]].values/uv_scalar)
42+
input_xr[uv_varnames[1]] = (('y','x'),input_xr[uv_varnames[1]].values/uv_scalar)
43+
44+
45+
#################### We initialize the objects for the flowline calculation
46+
flowlines = []
47+
48+
#################### Here is the forward calculation
49+
if forward0_both1_backward2 <= 1:
50+
temp_xs = np.expand_dims(seed_points[:,0],0)
51+
temp_ys = np.expand_dims(seed_points[:,1],0)
52+
53+
for ind0 in tqdm.tqdm(np.arange(steps)):
54+
x_search = xr.DataArray(temp_xs[-1,:],dims=['vector_index'])
55+
y_search = xr.DataArray(temp_ys[-1,:],dims=['vector_index'])
56+
new_u = input_xr[uv_varnames[0]].sel(x=x_search,y=y_search,method='nearest')
57+
new_v = input_xr[uv_varnames[1]].sel(x=x_search,y=y_search,method='nearest')
58+
59+
######### This is an order of magnitude slower
60+
#new_u = input_xr[uv_varnames[0]].interp(x=x_search,y=y_search)
61+
#new_v = input_xr[uv_varnames[1]].interp(x=x_search,y=y_search)
62+
63+
temp_xs = np.concatenate([temp_xs,temp_xs[-1:,:]+new_u.values.T*ds])
64+
temp_ys = np.concatenate([temp_ys,temp_ys[-1:,:]+new_v.values.T*ds])
65+
66+
xs = temp_xs
67+
ys = temp_ys
68+
else:
69+
xs = np.empty([2,1])
70+
ys = np.empty([2,1])
71+
72+
73+
#################### Here is the backward calculation
74+
if forward0_both1_backward2 >= 1:
75+
temp_xs = np.expand_dims(seed_points[:,0],0)
76+
temp_ys = np.expand_dims(seed_points[:,1],0)
77+
78+
for ind0 in tqdm.tqdm(np.arange(steps)):
79+
x_search = xr.DataArray(temp_xs[-1,:],dims=['vector_index'])
80+
y_search = xr.DataArray(temp_ys[-1,:],dims=['vector_index'])
81+
new_u = input_xr[uv_varnames[0]].sel(x=x_search,y=y_search,method='nearest')
82+
new_v = input_xr[uv_varnames[1]].sel(x=x_search,y=y_search,method='nearest')
83+
84+
######### This is an order of magnitude slower
85+
#new_u = input_xr[uv_varnames[0]].interp(x=x_search,y=y_search)
86+
#new_v = input_xr[uv_varnames[1]].interp(x=x_search,y=y_search)
87+
88+
temp_xs = np.concatenate([temp_xs,temp_xs[-1:,:]-new_u.values.T*ds])
89+
temp_ys = np.concatenate([temp_ys,temp_ys[-1:,:]-new_v.values.T*ds])
90+
91+
xs = np.concatenate([np.flipud(temp_xs),xs])
92+
ys = np.concatenate([np.flipud(temp_ys),ys])
93+
94+
95+
96+
flowlines = []
97+
for ind0 in np.arange(len(xs[0,:])):
98+
xy = np.stack([xs[:,ind0],ys[:,ind0]]).T
99+
flowlines.append(xy)
100+
101+
return flowlines
102+
103+
##########################################################################################
104+
#### This version of the code doesn't work quite right...
105+
##########################################################################################
106+
##def calculate_flowlines(x,y,u,v,seed_points,max_error=0.00001,retry_count_threshold=10):
107+
## """
108+
## % (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
109+
## %
110+
## % This function prints out the minimum and maximum values of an array
111+
## %
112+
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113+
## % The inputs are:
114+
## %
115+
## % input_array -- array of data to analyze
116+
## %
117+
## %%%%%%%%%%%%%%%
118+
## % The outputs are:
119+
## %
120+
## % output -- the min and max in a 1x2 array
121+
## %
122+
## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123+
## """
124+
##
125+
## ################# This uses a modified plt.streamline to pass through a user-editable keyword
126+
## ################# argument "max_error", which goes into the interpolater to guarantee
127+
## ################# accurate streamline calculation. Copy the updated version of streamplot
128+
## ################# into your matplotlib directory to enable the use of streamline
129+
## ################# calculation from NDH_Tools (in streamline.py, which calls _integrate_rk12)
130+
##
131+
## if isinstance(seed_points,list):
132+
## seed_points = np.array(seed_points)
133+
##
134+
## if len(seed_points.shape) == 1:
135+
## seed_points = np.expand_dims(seed_points,axis=0)
136+
##
137+
## ###################### Initialize the returned object
138+
## final_sls = []
139+
## for ind0 in np.arange(len(seed_points[:,0])):
140+
## final_sls.append([])
141+
##
142+
## retry_count = 0
143+
## retry_inds = np.arange(0,len(seed_points[:,0]))
144+
## seed_subset = seed_points
145+
##
146+
## while len(retry_inds) > 0:
147+
##
148+
## sls = []
149+
##
150+
## ################# Calculate the streamlines for all unfound seed points
151+
## fig = plt.figure()
152+
## if retry_count == 0:
153+
## print('The initial streamline calculation -- this can be slow. Finding '+str(len(seed_subset[:,0]))+' streamlines')
154+
## try:
155+
## streamlines = plt.streamplot(x,y,u,v,start_points=seed_points, max_error=max_error, density=100)
156+
## except:
157+
## streamlines = plt.streamplot(x,y,u,v,start_points=seed_points, density=100)
158+
## if retry_count == 0:
159+
## print('Note: You need to update your matplotlib streamline.py and reduce the max error for this to work properly')
160+
## plt.close(fig)
161+
##
162+
## ################# Here we extract the coordinate info from the streamlines
163+
## sl_deconstruct = []
164+
## for i in streamlines.lines.get_paths():
165+
## sl_deconstruct.append(i.vertices[1])
166+
## sl_deconstruct = np.array(sl_deconstruct)
167+
##
168+
## ################ Here we separate the streamlines based on large breaks in distance
169+
## sl_dist = ndh.distance_vector(sl_deconstruct[:,0],sl_deconstruct[:,1],1)
170+
## dist_mean = np.mean(sl_dist)
171+
## breaks = np.where(sl_dist > (dist_mean+1)*50)[0]
172+
## if len(breaks) > 0:
173+
## breaks = np.concatenate([np.array([-1]),breaks,np.array([len(sl_deconstruct[:,0])])])+1
174+
## else:
175+
## breaks = np.array([0,len(sl_deconstruct[:,0])+1])
176+
##
177+
## for ind0 in np.arange(len(breaks)-1):
178+
## sls.append(sl_deconstruct[breaks[ind0]:breaks[ind0+1],:])
179+
##
180+
## ################ Here we identify which streamline goes with which seed_point
181+
## matching = []
182+
## for ind0 in np.arange(len(seed_subset[:,0])):
183+
## dists = []
184+
## for ind1,sl in enumerate(sls):
185+
## comp_vals = ndh.find_nearest_xy(sl,seed_subset[ind0,:])
186+
## dists.append(comp_vals['distance'][0])
187+
## best = np.where(np.array(dists) < 1e-8)[0]
188+
## try:
189+
## matching.append(best[0])
190+
## except:
191+
## matching.append(-1)
192+
##
193+
## ################# populate the final object
194+
## for ind0,i in enumerate(matching):
195+
## if i != -1:
196+
## final_sls[retry_inds[ind0]] = sls[i]
197+
##
198+
## ################# Finally, we identify the new set of streamlines that need to be computed, based on which have no match
199+
## new_retry_inds = np.where(np.array(matching) == -1)[0]
200+
## seed_subset = seed_points[retry_inds[new_retry_inds],:]
201+
## retry_inds = retry_inds[new_retry_inds]
202+
##
203+
## if len(retry_inds) > 0:
204+
## retry_count = retry_count+1
205+
## print('Recalculating for nearly overlapping points -- try '+str(retry_count)+'. Finding '+str(len(seed_subset[:,0]))+' streamlines')
206+
##
207+
## if retry_count > retry_count_threshold:
208+
## break
209+
##
210+
## if 0:
211+
## plt.figure()
212+
## plt.plot(test_dist)
213+
## plt.axhline(dist_median,c='orange')
214+
##
215+
## if 0:
216+
## plt.figure()
217+
## plt.plot(test[:,0],test[:,1],c='blue')
218+
## for i in final_sls:
219+
## plt.plot(i[:,0],i[:,1],c='red')
220+
## plt.plot(seed_points[:,0],seed_points[:,1],'o')
221+
##
222+
##
223+
## return final_sls

‎compare_list.py

100644100755
File mode changed.

‎edgetrim_mask.py

100644100755
File mode changed.

‎find_pixelcoords.py

100644100755
+11-3
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0
4040
import NDH_Tools as ndh
4141
import numpy as np
4242
import cv2
43+
import tqdm
4344

4445

4546
im_handle = Image.open(im_filename)
@@ -50,16 +51,16 @@ def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0
5051
######################### This method only works if the image is perfectly white outside of the axes
5152
if 0:
5253
im_frame = np.where(np_frame[:,:,3] == 255)
53-
54+
5455
cinds = ndh.minmax(im_frame[1])
5556
rinds = ndh.minmax(im_frame[0])
5657
else:
5758
########################## This is meant to handle some spillover of lines and points outside the image
5859
row_sum = np.sum(np_frame[:,:,2] == 255,axis=1)
5960
col_sum = np.sum(np_frame[:,:,2] == 255,axis=0)
6061

61-
rinds = ndh.minmax(np.where(row_sum > np.mean(row_sum)))
62-
cinds = ndh.minmax(np.where(col_sum > np.mean(col_sum)))
62+
rinds = ndh.minmax(np.where(row_sum <= np.percentile(row_sum,80)-1))
63+
cinds = ndh.minmax(np.where(col_sum <= np.percentile(col_sum,80)-1))
6364

6465
xrange = np.linspace(0,original_width,cinds[1]-cinds[0])
6566
yrange = np.linspace(0,original_height,rinds[1]-rinds[0])
@@ -89,6 +90,12 @@ def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0
8990
############ Here, we post-process the picks to deal with gaps / multiple layers, etc.
9091
#######################################################################################
9192
########### Finally, we loop through the edge_trims and the picks, and pull out just the relevant info
93+
94+
#print('--- Total number of countoured objects to process: %d' %len(process_contours))
95+
mod_step = 10**np.floor(np.log10(len(process_contours)/10));
96+
97+
98+
#for ind2 in tqdm.tqdm(range(len(process_contours))):
9299
for ind2 in range(len(process_contours)):
93100

94101
pick_img = np.zeros_like(process_im)
@@ -213,6 +220,7 @@ def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0
213220
pick_output.append(pick_temp)
214221
else:
215222
pick_output.append(pick_temp)
223+
216224

217225

218226
return pick_output

‎generate_pickingpdf.py

100644100755
+14-11
Original file line numberDiff line numberDiff line change
@@ -71,20 +71,21 @@ def generate_pickingpdf(fn,picking_root_dir,frame_spacing=25,surf_dir='CSARP_sur
7171
fn_list[-3] = surf_dir
7272
fn2 = '/'.join(fn_list)
7373

74-
try:
75-
surfdata = ndh.loadmat(fn2)
76-
surf_dims = surfdata['surf']['y'][1].shape
77-
bot_ind = ndh.find_nearest(data['Time'],np.max(surfdata['surf']['y'][1]))
78-
except:
79-
print(fn2+' could not be found')
80-
fn2 = 0
81-
bot_ind = {'index':[len(data['Time'])-25]}
74+
#try:
75+
surfdata = ndh.loadmat(fn2)
76+
surf_dims = surfdata['surf']['y'][1].shape
77+
bot_ind = ndh.find_nearest(data['Time'],np.array([np.max(surfdata['surf']['y'][1])]))
78+
bot_ind['index'] = [np.min([bot_ind['index'][0]+150,len(data['Time'])])]
79+
#except:
80+
# print(fn2+' could not be found')
81+
# fn2 = 0
82+
# bot_ind = {'index':[len(data['Time'])]}
8283

8384
#### Now we loop through the frames we want to plot and generate an image for
8485
frame_print = np.arange(0,len(xy['x']),frame_spacing)
8586
for ind1,i in enumerate(frame_print):
8687
ndh.remove_image(ax,1,verbose=0)
87-
ax.imshow(np.squeeze(np.log10(data['Tomo']['img'][:bot_ind['index'][0]+25,:,i])),cmap='bone_r')
88+
ax.imshow(np.squeeze(np.log10(data['Tomo']['img'][:bot_ind['index'][0],:,i])),cmap='bone_r')
8889
ax.set_aspect('auto')
8990
if fn2 != 0:
9091
ndh.remove_line(ax,1,verbose=0)
@@ -95,7 +96,7 @@ def generate_pickingpdf(fn,picking_root_dir,frame_spacing=25,surf_dir='CSARP_sur
9596
pass
9697

9798
plt.axis('off')
98-
plt.savefig('%s%s/%s/Frame_%0.4d_fs_%0.2d_crop_%0.4d.png' %(picking_root_dir,seg,frame,i,frame_spacing,bot_ind['index'][0]+25))
99+
plt.savefig('%s%s/%s/Frame_%0.4d_fs_%0.2d_crop_%0.4d.png' %(picking_root_dir,seg,frame,i,frame_spacing,bot_ind['index'][0]))
99100

100101

101102
print('Completed the image generation')
@@ -105,7 +106,7 @@ def generate_pickingpdf(fn,picking_root_dir,frame_spacing=25,surf_dir='CSARP_sur
105106

106107
########## This converts all the images to a single pdf
107108
pdfroot = picking_root_dir+'To_Pick/'+seg+'/'
108-
pdfend = '%s_fs_%0.2d_crop_%0.4d.pdf' %(frame,frame_spacing,bot_ind['index'][0]+25)
109+
pdfend = '%s_fs_%0.2d_crop_%0.4d.pdf' %(frame,frame_spacing,bot_ind['index'][0])
109110
pdfname=pdfroot+pdfend
110111

111112
frames = '%s%s/%s/*.png' % (picking_root_dir,seg,frame)
@@ -145,7 +146,9 @@ def generate_pickingpdf(fn,picking_root_dir,frame_spacing=25,surf_dir='CSARP_sur
145146
max_num = len(file_list)
146147
seg = file_list[-1].split('/')[-2]
147148
froot = seg
149+
148150
num_range = np.arange(1,max_num+1,1)
151+
print('Images from '+str(max_num)+' files')
149152

150153

151154

‎get_contour_data.py

100644100755
File mode changed.

0 commit comments

Comments
 (0)