bugfix: use plot_config in plot_cmap

track_analyzer/plotting.py

@@ -990,7 +990,7 @@ def plot_all_traj(data_dir, df, image={'image_fn': None, 't_dim': None, 'z_dim':
             label = tpr.make_param_label(color_code, l_unit=info["length_unit"], t_unit=info["time_unit"])
         if traj_parameters['cmap_lim'] is not None:
             plot_cmap(plot_dir, label, traj_parameters['cmap'], traj_parameters['cmap_lim'][0],
-                  traj_parameters['cmap_lim'][1],suffix='_traj')
+                  traj_parameters['cmap_lim'][1],suffix='_traj',plot_config=plot_config)
     
     # make a colmuns of indices to be used for color_cycle
     elif color_code == "group": 
@@ -1061,7 +1061,7 @@ def plot_all_scalar_fields(data_dir, df, data, field, image={'image_fn': None, '
         label = tpr.make_param_label(None, manual_symbol=custom_var[field]['name'], manual_unit=custom_var[field]['unit'])
     else: 
         label = tpr.make_param_label(field, l_unit=info['length_unit'], t_unit=info['time_unit'])
-    plot_cmap(plot_dir, label, map_param_['cmap'], map_param_['vlim'][0], map_param_['vlim'][1],suffix='_'+field)
+    plot_cmap(plot_dir, label, map_param_['cmap'], map_param_['vlim'][0], map_param_['vlim'][1],suffix='_'+field,plot_config=plot_config)
 
     stack = []  # stack to store images
     for frame in np.sort(df['frame'].unique()):
@@ -1118,7 +1118,7 @@ def plot_all_vector_fields(data_dir, df, data, field, plot_on_field=None, dim=3,
                 label = tpr.make_param_label(None, manual_symbol=custom_var[plot_on_field['plot_on']]['name'], manual_unit=custom_var[plot_on_field['plot_on']]['unit'])
             else: 
                 label = tpr.make_param_label(plot_on_field['plot_on'], l_unit=info['length_unit'], t_unit=info['time_unit'])
-            plot_cmap(plot_dir, label, plot_on_field['cmap'], plot_on_field['vlim'][0], plot_on_field['vlim'][1],suffix='_vector_'+plot_on_field['plot_on'])
+            plot_cmap(plot_dir, label, plot_on_field['cmap'], plot_on_field['vlim'][0], plot_on_field['vlim'][1],suffix='_vector_'+plot_on_field['plot_on'],plot_config=plot_config)
 
     stack = []  # stack to store images
     for frame in np.sort(df['frame'].unique()):
@@ -1190,7 +1190,7 @@ def plot_all_Voronoi(data_dir, df, data, show_local_area=True, df_mean = None,
 
     if show_local_area:
         label = tpr.make_param_label('area', l_unit=info['length_unit'], t_unit=info['time_unit'])
-        plot_cmap(plot_dir, label, map_param['cmap'], map_param['vlim'][0], map_param['vlim'][1],suffix='_area')
+        plot_cmap(plot_dir, label, map_param['cmap'], map_param['vlim'][0], map_param['vlim'][1],suffix='_area',plot_config=plot_config)
 
     stack = []  # stack to store images
     for frame in np.sort(df['frame'].unique()):
@@ -1425,7 +1425,7 @@ def plot_total_traj(data_dir, df, dim=3, plot_dir=None, plot_fn=None, plot_confi
         else: 
             label = tpr.make_param_label(color_code, l_unit=info["length_unit"], t_unit=info["time_unit"])
         if cmap_lim is not None:
-            plot_cmap(plot_dir, label, cmap, cmap_lim[0],cmap_lim[1],suffix='_total_traj')
+            plot_cmap(plot_dir, label, cmap, cmap_lim[0],cmap_lim[1],suffix='_total_traj',plot_config=plot_config)
     
     elif color_code == 'group':
         if 'subset' in df.columns:

track_analyzer/scripts/make_synthetic_data.py

@@ -49,22 +49,59 @@ frame_num = 5  # number of frame
 dim = 2  # dimension
 lims = [[0,500],[0,500]]  # space limits
 periodic = False  # space periodicity
-space = {'dim':dim,'lims':lims,'periodic':periodic}
+random = False  # initial position spacing
+space = {'dim':dim,'lims':lims,'periodic':periodic,'random':random}
 pos0_lims = [[100,400],[100,400]]  # initial position limits
 diff_field = {'cst':1}  # constant diffusion 
-vel_field = [{'a':20,'p':2,'b':0.5,'c':-20,'q':2,'d':0.5},  # vx = x^2 - y^2 (centered)
-             {'a':20,'p':2,'b':0.5,'c':-20,'q':2,'d':0.5},  # vy = x^2 - y^2 (centered)
-            ] 
+vel_field = [{'a':0.01,'p':2,'b':250,'c':-0.01,'q':2,'d':250},  # vx = x^2 - y^2 (centered)
+             {'a':0.01,'p':2,'b':250,'c':-0.01,'q':2,'d':250},  # vy = x^2 - y^2 (centered)
+            ]
 
+config_default = {'traj_num':traj_num, 
+                'frame_num':frame_num,
+                'space':space,
+                'pos0_lims':pos0_lims,
+                'diff_field':diff_field,
+                'vel_field':vel_field,    
+                }
 
-def make_pos0(pos0_lims=pos0_lims,traj_num=traj_num): 
+def get_config(data_dir): 
+    """
+    Load from config file 'config.csv' if it exists, else return None
+    """
+    config_fn = osp.join(data_dir,'config.csv')
+    if osp.exists(config_fn):
+        config = tpr.load_dict(config_fn)
+    else: 
+        config = None
+
+    return config
+
+
+def make_pos0(pos0_lims=pos0_lims,traj_num=traj_num,random=True): 
     """
     Generate list of initial positions at random position within boudaries given by pos0_lims 
     """
     pos0_list = []
-    for i in range(traj_num): 
-        pos0 = np.array([np.random.uniform(pos0_lims_[0], pos0_lims_[1]) for pos0_lims_ in pos0_lims])
-        pos0_list.append(pos0)
+    
+    if random:
+        for i in range(traj_num): 
+            pos0 = np.array([np.random.uniform(pos0_lims_[0], pos0_lims_[1]) for pos0_lims_ in pos0_lims])
+            pos0_list.append(pos0)
+    else:
+        # create a rectangle array of N_ positions evenly spaced 
+        # height/width = N_h/N_w and N_h*N_w = N_ 
+        width = np.abs(pos0_lims[0][1] - pos0_lims[0][0])
+        height = np.abs(pos0_lims[1][1] - pos0_lims[1][0])
+        N_ = int(np.sqrt(traj_num))**2  # number of positions needs to be a square number
+        N_h = int(np.sqrt(N_) * height/width)
+        N_w = int(N_ / N_h)
+
+        # meshgrid
+        X,Y = np.meshgrid(np.linspace(pos0_lims[0][0],pos0_lims[0][1],N_w),np.linspace(pos0_lims[1][0],pos0_lims[1][1],N_h))
+        for i in range(X.shape[0]):
+            for j in range(X.shape[1]):
+                pos0_list.append(np.array([X[i,j],Y[i,j]]))
     
     return pos0_list
 
@@ -72,8 +109,8 @@ def make_pos0(pos0_lims=pos0_lims,traj_num=traj_num):
 def polynom_mapping(pos,params={'cst':0},space=space):
     """
     Evaluate polynomial at position pos.
-    Polynomial function: a*(x-b)^p + c*(y-d)^q + e*(z-f)^r + cst
-    Each coordinate are normalized by the space size. 
+    Polynomial function: (a*(x-b))^p + (c*(y-d))^q + (e*(z-f))^r + cst
+    Useful to generate a velocity field with polynomial expression
     """
 
     # put parameter to either if missing
@@ -86,24 +123,19 @@ def polynom_mapping(pos,params={'cst':0},space=space):
     if pos_dim > space['dim']:
         raise Exception('position dimension larger than space dimension')
 
-    # normalize position
-    norm_pos = np.zeros(pos_dim)
-    for i in range(pos_dim):
-        norm_pos[i] =  pos[i] / space['lims'][i][1]
-
     value = params['cst']
 
     if pos_dim == 1:
-        value += params['a'] * (norm_pos[0] - params['b']) ** params['p']
+        value += (params['a'] * (pos[0] - params['b'])) ** params['p']
 
     elif pos_dim == 2:
-        value += params['a'] * (norm_pos[0] - params['b']) ** params['p']
-        value += params['c'] * (norm_pos[1] - params['d']) ** params['q']
+        value += (params['a'] * (pos[0] - params['b'])) ** params['p']
+        value += (params['c'] * (pos[1] - params['d'])) ** params['q']
 
     elif pos_dim == 3:
-        value += params['a'] * (norm_pos[0] - params['b']) ** params['p']
-        value += params['c'] * (norm_pos[1] - params['d']) ** params['q']
-        value += params['e'] * (norm_pos[2] - params['f']) ** params['r']
+        value += (params['a'] * (pos[0] - params['b'])) ** params['p']
+        value += (params['c'] * (pos[1] - params['d'])) ** params['q']
+        value += (params['e'] * (pos[2] - params['f'])) ** params['r']
 
     return value
 
@@ -124,6 +156,70 @@ def n_polynom_mapping(pos,params_list=[{'cst':0}],space=space):
     return vector
 
 
+def export_vel_fields(vel_field,pos_list,data_dir,space=space):
+    """
+    Export velocity curl and div fields evaluated at an list of positions.
+    vel_field (generated by n_polynom_mapping) is a list of n elements (n=2 or 3), each element being a spatial component
+    """
+
+    # check 2D velocity field is 2D
+    if len(vel_field) < 2: 
+        raise Exception("Need a 2D velocity field")
+
+    # check missing parameters in vel_field
+    for i in range(len(vel_field)):
+        for p in list('abcdefpqr'):
+            if p not in vel_field[i].keys():
+                vel_field[i][p] = 0
+
+    # get velocity list
+    vel_list = []
+    for pos in pos_list: 
+        vel_list.append(n_polynom_mapping(pos,params_list=vel_field,space=space))
+    
+    # calculate 2D curl and div
+    curl_list = []
+    div_list = []
+    for pos in pos_list: 
+        # curl
+        curl = dx_polynom(vel_field[1],pos) - dy_polynom(vel_field[0],pos)  # Dx_vy - Dy_vx 
+        curl_list.append(curl)
+        # div
+        div =  dx_polynom(vel_field[0],pos) + dy_polynom(vel_field[1],pos)  # Dx_vx + Dy_vy
+        div_list.append(div)
+
+    # save to df 
+    data = np.concatenate((np.array(pos_list),np.array(vel_list),np.array([curl_list]).T,np.array([div_list]).T),axis=1)
+
+    columns = ['x','y','vx','vy','vz','curl','div'] if len(vel_field) == 3 else ['x','y','vx','vy','curl','div']
+    df_out = pd.DataFrame(data, columns=columns)
+    df_out.to_csv(osp.join(data_dir,'fields.csv'))
+
+
+def dx_polynom(polynom,pos):
+    """
+    x-component of spatial derivative of polynom_mapping at position pos
+    """
+    if polynom['p'] > 0:
+        value = polynom['a']*polynom['p']*(pos[0]-polynom['b'])**(polynom['p']-1)
+    else:
+        value = 0
+
+    return value
+
+
+def dy_polynom(polynom,pos):
+    """
+    y-component of spatial derivative of polynom_mapping at position pos
+    """
+    if polynom['q'] > 0:
+        value = polynom['c']*polynom['q']*(pos[1]-polynom['d'])**(polynom['q']-1)
+    else:
+        value = 0
+
+    return value
+
+
 def make_traj(pos0,frame_num=frame_num,diff_field=diff_field,vel_field=vel_field,space=space,track_id=0):
     """
     Generate a trajectory of frame_num positions from a an initial position pos0. 
@@ -220,7 +316,7 @@ def parse_args(args=None):
     parser.add_argument('-p', '--plot',
                         action="store_true",
                         default=False,
-                        help='refresh database')
+                        help='plot positions and save to tif stack')
 
     parsed_args = parser.parse_args(args)
 
@@ -236,15 +332,23 @@ def main(args=None):
     data_dir = osp.realpath(parsed_args.data_dir)
     plot = parsed_args.plot
 
-    pos0_list = make_pos0(pos0_lims=pos0_lims,traj_num=traj_num)
+    # get config 
+    config_ = get_config(data_dir)
+    config = config_default if config_ is None else config_
+
+    traj_num = config["traj_num"]  # number of tracks
+    frame_num = config["frame_num"]  # number of frame
+    space = config["space"]
+    pos0_lims = config["pos0_lims"]  # initial position limits
+    diff_field = config["diff_field"]  # constant diffusion 
+    vel_field = config["vel_field"]
+
+    pos0_list = make_pos0(pos0_lims=pos0_lims,traj_num=traj_num,random=space['random'])
     df = make_dataset(pos0_list,frame_num=frame_num,diff_field=diff_field,vel_field=vel_field,space=space)
     df.to_csv(osp.join(data_dir,'positions.csv'), index=False)
 
-    #save parameters
-    parameters = {'dim':dim,
-                'lims':lims,
-                'periodic':periodic,
-                'space':space,
+    # save parameters
+    parameters = {'space':space,
                 'pos0_lims':pos0_lims,
                 'traj_num':traj_num,
                 'diff_field':diff_field,