# Surface tracker

Back to CReSIS Toolbox Guide Tracking 2D and 3D Images or Maintenance and Utility Functions sections.

# Overview

The surface tracker is used to track the surface of the Earth in sounding data. It is generally called from the qlook.m process. It may be called directly from run_track_layer.m.

The surface tracker guide is one part of the CReSIS Toolbox Guide.

It supports basic-snake (tracker_snake_simple.m), max (tracker_max.m), and threshold (tracker_threshold.m) trackers. It may be run in a debug mode so that each tracked frame can be monitored as the tracking is done.

If the layer just needs to be shifted up/down or left/right and the layer data are stored in the OPS, use runOpsShiftTwtt.m or runOpsShiftGpsTime.m respectively.

If the echogram needs to be shifted up/down then correct param.radar.wfs(wf).Tadc_adjust and either recreate the data products (ideal) or use run_update_surface_twtt_delta.m to update existing data products (more efficient, but may affect SAR focusing if the time-shift is large). If the echogram needs to be shifted left/right, then correct param.vectors.gps.time_offset and run run_update_records.m.

To view tracked layers on echograms use imb.picker if the layer data are in OPS or run_picker is the layer data are in layerData files.

## layer_tracker.m

This function updates the surface layer. There are several choices of automated trackers, but "threshold" is usually the best. The function is called from run_layer_tracker.m. The parameter spreadsheet structure is used to determine which segments and frames are processed and how they are processed.

The param.layer_tracker field controls how the surface is tracked. The layer_tracker structure has the following fields:

Field Description
debug_plots Cell array of strings. Default is an empty cell array. If cell array contains "debug" string, then the tracker will stop in debug mode after each frame is tracked and display the echogram with the tracked surface. Type "dbcont" to continue.
echogram_img Specifies which echogram image to use for tracking. Set to 0 to use the combined image (e.g. "Data_20170330_01_001.mat') or to a positive integer to use one of the images (e.g. setting to 2 would load "Data_img_02_20170330_01_001.mat').
echogram_source String containing the source to use. This source will be passed into ct_filename_out. Setting "qlook" will look in the CSARP_qlook directory for the echogram. Setting "CSARP_post/standard" will look in the CSARP_post/CSARP_standard directory for the echogram.
layer_params Structure that specifies where the tracked surface will be stored. The structure is similar to what opsLoadLayers and opsCopyLayers uses.
track The track structure controls the tracking.

### track structure

The track structure specifies how the surface tracker operates. The fields are:

Field Description
data_noise_en Logical scalar. Default is false. If true, then a matrix data_noise will be created which will go through all the same operations as data except no "sidelobe" and no "feedthru" masking will be done. This is sometimes necessary to enable when sidelobe or feedthru remove too much data so that the noise estimation process in the tracker does not function properly. Currently only tracker_threshold makes use of data_noise.
debug_time_guard Vertical band around layer in debug plot
detrend Either a string or a positive integer scalar. If an integer, then the average range line power is fit to a polynomial of that order and this is subtracted from each range line. If a string, then the detrend file generated by run_get_echogram_stats.m is used for detrending (e.g. ct_tmp/echogram_stats/snow/2011_Greenland_P3/stats_20110329_02.mat).
en Set to true to enable the tracker. Most relevant when running from qlook. Default is true (enabled).
feedthru Structure array if to apply a mask to the data (usually a mask to ignore the direct path and/or feed through signal). Default is empty. The feedthru structure specifies the feed through power level. This is useful when the surface return comes at the same time as the feed through. This assumes that the feed through signal is the same for each range line. Default is disabled, [].
filter Two element vector that specifies the number of incoherent averages (multilooks) in the fast-time (first element) and along-track (second element) dimensions. This is especially useful when the SNR is low. The fast-time incoherent averaging is often set to 1 (disabled) and the along-track is set to values between 3 and 7. Default is [1 1] or no filtering.
filter_trim Two element vector that specifies the number of range bins to trim off the top and bottom of the image after filtering.
fixed_value The value to be used for the surface two way travel time when the track.method == fixed.
init Structure specifying the initial surface estimate.
max_bin Sets the maximum two way travel time that the surface is allowed to be at. This is useful when there are large signals beneath the surface that confuse the tracker.m. Default is inf.
max_rng Specifies a search range around the tracked layer to search for the peak/max value. The units are specified by max_rng_units. Default is [0 0] which causes the automated tracker result to be returned. For threshold tracking, usually this is set to [0 T] where T is chosen based on the time step, dt, and the the track.filter duration. For example, if td is the automated result for a particular range line, then td will be adjusted to the peak value between [td+0 and td+T].
max_rng_units String specifying the units of max_rng. Default is 'time' and max_rng should be specified in seconds. If 'bins', then the max_rng should be in range bins.
medfilt Scalar positive integer. Layer is median filtered with a robust version of medfilt1.
medfilt_threshold Scalar double. Only updates the output if the median filtered surface if different by >medfilt_threshold bins. The default is zero so that all range lines are updated by the median filter. Setting to inf would cause no range lines to be updated by the median filter.
method String describing the surface tracking method to use. Default is ''.
min_bin Sets the minimum two way travel time that the surface is allowed to be at. This is useful when there is a strong feed through signal at the beginning of the record/range-line, but the levels are not known precisely enough to use the feedthru mask feature. If the surface is closer than min_bin, layer_tracker.m will not be able to track the surface properly. This is usually set to values from 1e-6 to 2e-6 for low altitude terrain following at 1500 ft AGL. Should be zero for ground operation. Default is zero.
profile String. Default is "default". This loads default settings for the track structure. See layer_tracker_profile.m to see which profiles exist. Typical values are "RDS_OIB" for rds radars, "ACCUM" for accum radars, and "SNOW_AWI" for ka-band, ku-band, and snow radars.
search_rng After the surface is tracked, the surface is adjusted to find the first peak at or after the tracked value. This specifies the number of bins after the tracked bin to look for a peak. If search_rng is set to 0, then no adjustment is made. Usually set to values between "0" and "0:9". Default is zero.
trim 2 element vector [Nstart Nend]. Default is [0 0]. The first nonzero Nstart bins will be set to zero and the last nonzero Nend bins will be set to zero.

### Surface tracking methods

The surface tracking methods are described below:

Method String Function Description
'' NA Uses the results of the initial surface estimate.
'fixed' NA Sets the layer to a fixed two way travel time.
'max' tracker_max.m Chooses the bin with the maximum power.
'snake' tracker_snake_simple.m A very simple max power tracker.
'threshold' tracker_threshold.m Uses a noise estimator to determine the noise floor and then finds the first range bin that exceeds the noise floor by a threshold subject to a number of constraints.
'viterbi' tracker_viterbi.m Uses the Viterbi algorithm to track the surface.

### snake method

The snake surface tracking parameters are described below:

Field Description
search_rng Range of bins relative to the last bin to search for the new surface. Should be symmetric (e.g. -90:90 or -5:5).

### threshold method

The threshold surface tracking parameters are described below:

Field Description
init The init structure describes the initialization step of the tracker. This step tries to find an initial surface. Default is 'max'.
noise_rng Three element vector which specifies the region relative to the peak power to estimate the noise power from. It is a 1 by 3 integer array.

first element: All data points that are zero at the beginning of the record will be ignored in the noise calculation. This specifies a buffer beyond this in sample bins that will also be ignored.
second and third elements: Specifies a range of bins relative to the max bin that will be used to estimate the noise (usually these are negative such as [-50 -10] so that the noise estimate uses data before the peak).

threshold Scalar value in dB. This is added to the noise floor estimate to determine the threshold. The tracker looks for the first bin in each range line that exceeds the threshold. The search bins are limited to those bins within "max_diff" seconds of the initial surface though. If the initial surface is off by a lot and max_diff is too small, then the tracker may not be able to track the surface.

### Viterbi method

The Viterbi surface tracking parameters are described below. Note that all arguments involving indices should be one-indexed. Viterbi.cpp will handle the shift to zero-indexed arrays for internal use and then handle the return back to one-indexing with the output layer passed back to Matlab.

Field Description
input_img (image) MxN matrix of doubles representing rows and columns of the data on which to run Viterbi.
surface_gt (sgt) 1xN vector of doubles representing the range bins of surface ground truth (one-indexed).
extra_gt (egt) 2xNgt (Ngt is arbitrary) matrix of doubles containing ground truth points. First row corresponds to range lines while the second row corresponds to range bins (both one-indexed).
ice_mask (mask) 1xN (from image) vector of doubles representing the proximity of ice for each column. 0 indicates no ice while inf indicates ice presence in entire vicinity. Other positive numbers indicate a transition between icy and non-icy regions.
img_mag_weight A double representing the weight by which to multiply the image magnitude cost. Greater weight = prefer greater image magnitude.
smooth_slope 1xN-1 vector of doubles representing the slope of the target layer. Slope generally occurs due to the plane changing altitude. 0 for no expected slope.
max_slope A double representing the maximum allowed slope of the target layer based on distance between range bin indices. -1 for no max.
bounds 1x2 vector of ints representing the beginning and ending indices (one-indexed). Default for each index, [1 Nx] (actually [0 Nx-1] in Viterbi.cpp but the input bounds are expected to be one-indexed), is set respectively for each negative index passed or if passed an empty vector.
gt_weights 1xN vector of doubles representing weighting factor for all points. Applied to extra ground truth points. Greater weight = prefer closer to ground truth.
gt_cutoffs 1xN vector of doubles representing max range bin distance from gt allowed. Applied to extra ground truth points. -1 for any distance from gt.
mask_dist 1xN vector of doubles representing the distance from the nearest ice-mask. Areas with no ice (mask of zero) should contain distance to nearest ice.
costmatrix MxN matrix of doubles representing the costs for the Distance-to-Ice-Margin model.
transition_weights 1xN-1 vector of doubles representing the weight by which to multiply the binary cost. Greater weight = prefer less sloping.
surf_weight A double representing surface repulsion. A greater surf_weight means a greater cost for points on the surface layer.
mult_weight A double representing surface multiple repulsion. A greater mult_weight means a greater cost for points on expected surface multiple layers.
mult_weight_decay A double which multiplies the mult_weight for each subsequent multiple. i.e. the rate of cost decay for following multiples to compensate for a fainter return.
mult_weight_local_decay A double which multiplies the mult_weight for each subsequent range bin away from the nearest multiple. e.g. when five range bins past the second expected multiple, the multiple cost will have been multiplied by mult_weight_decay^1 and then by mult_weight_local_decay^5.
[zero_bin] An int64 which signifies the range bin occupied by the plane (one-indexed). Used in anticipating the location of surface multiples. Defaults to 0 if argument left out.

### Initial surface tracker

The initial surface tracker fields are:

Field Description
method String describing the initial tracker method. See table below for valid methods. If empty, then just the max value in each range line will be used.
medfilt Scalar used by medfilt method. Specifies the length of the medfilt1 filter.
lidar_source String specifying the lidar source to use (similar to opsLoadLayers). Valid LIDAR sources are "atm", "awi".

Initial tracker methods are:

Method Description
dem Uses land and sea surface digital elevation models and LIDAR data (if lidar_source specified) to determine the initial surface. If elevation or LIDAR data are not available, this initial tracker will perform very badly. The LIDAR data has highest priority.
lidar Uses LIDAR data (if lidar_source specified) to determine the initial surface. If LIDAR data are not available, this initial tracker will perform very badly.
medfilt Same as max except that the result is passed through a median filter (medfilt1).
snake A very simple max power tracker.

### Feed through

All values in the radar image which exceed the power mask set by the time and power_dB vectors will be set to zero power. The power mask is linearly interpolated. To get these values a typical procedure is:

plot(mdata.Time, lp(mean(mdata.Data,2)))
[surf.feedthru.time,surf.feedthru.power_dB] = ginput; % Press enter to end the input


The fields of the feed through structure are:

Method Description
power_dB N length vector of power in dB corresponding to each time.
time N length vector of two way travel times