Merge pull request #192 from CReSIS/master

Master

cresis-toolbox/display/echo_detrend.m

@@ -225,7 +225,7 @@
       if top_bin > Nt || bottom_bin < 1
         continue;
       end
-      bins = top_bin:bottom_bin;
+      bins = round(top_bin):round(bottom_bin);
 
       if length(bins) >= 2
         mask(bins,rline) = isfinite(data(bins,rline));
@@ -243,8 +243,8 @@
     trend = zeros(Nt,1);
     for rline = 1:Nx
       % Section 2: Evaluate the polynomial
-      top_bin = param.layer_top(rline);
-      bottom_bin = param.layer_bottom(rline);
+      top_bin = round(param.layer_top(rline));
+      bottom_bin = round(param.layer_bottom(rline));
       if top_bin < 1 && bottom_bin >= 1
         top_bin = 1;
       end
@@ -305,9 +305,9 @@
         figure(1); clf;
         plot(data(:,rline));
         if param.roll_comp_en
-          title(sprintf('Rangeline %d block %d, Roll= %.3f',rline,param.block,roll(rline) ),'Interpreter','none');
+          title(sprintf('Rangeline %d Roll= %.3f',rline,roll(rline) ),'Interpreter','none');
         else
-          title(sprintf('Rangeline %d block %d',rline,param.block),'Interpreter','none');
+          title(sprintf('Rangeline %d',rline),'Interpreter','none');
         end
         hold on;
         plot(trend);

cresis-toolbox/display/echo_norm.m

@@ -40,10 +40,10 @@
 end
 
 % valid_max_range_dB: 2 element numeric vector; specifies the valid range
-% for the max value; default is [-inf inf] which effectively disables
-% this valid max value range constraint;
+% for the max value; default is [-inf inf]. Use [-inf inf] to effectively
+% disable this valid max value range constraint;
 if ~isfield(param,'valid_max_range_dB') || isempty(param.valid_max_range_dB)
-  param.valid_max_range_dB = [30 inf];
+  param.valid_max_range_dB = [-inf inf];
 end
 
 if ~isfield(param,'scale') || isempty(param.scale)
@@ -67,4 +67,4 @@
   max(data(:))));
 
 % Scale and offset data
-data = param.scale(1) + (data-noise) / (max_scalar - noise) * param.scale(2);
+data = param.scale(1) + (data-noise) / (max_scalar - noise) * (param.scale(2)-param.scale(1));

cresis-toolbox/matlab_speed_test.m

@@ -45,11 +45,17 @@
   % gpuDeviceTable
   for gpu_idx = 1:gpuDeviceCount('available')
     gpu_dev = gpuDevice(gpu_idx);
-    fprintf('GPU\t%s\n', gpu_dev.Name);
-
+    % 8 bytes per sample, 3 copies in memory, 10e3 rows, 80% utilization
+    cols = min(10e3,floor(gpu_dev.TotalMemory / 10e3 / 8 / 3 * 0.8));
+    if cols == 10e3
+      fprintf('GPU\t%s\n', gpu_dev.Name);
+    else
+      fprintf('GPU\t%s: only using %d instead of 10000 columns due to VRAM limitation\n', gpu_dev.Name, cols);
+    end
+
     start_time = tic;
     fprintf('fft_start\t%g\n', toc(start_time));
-    A = randn(10e3,10e3) + 1i*randn(10e3,10e3);
+    A = randn(10e3,10e3,'single') + 1i*randn(10e3,10e3,'single');
     % gpurng(0, 'Philox');
     % B = randn(10e3,5e3,'single','gpuArray') + 1i*randn(10e3,5e3,'single','gpuArray');
     B = gpuArray(A);

cresis-toolbox/processing/collate_equal.m

@@ -29,7 +29,6 @@
 if ~isfield(param.collate_equal,'cmd_idx') || isempty(param.collate_equal.cmd_idx)
   param.collate_equal.cmd_idx = 1;
 end
-cmd = param.analysis.cmd{param.collate_equal.cmd_idx};
 
 if ~isfield(param.collate_equal,'debug_out_dir') || isempty(param.collate_equal.debug_out_dir)
   param.collate_equal.debug_out_dir = 'collate_equal';
@@ -156,6 +155,7 @@
   dt = waveform.dt;
   fc = waveform.fc;
   ref_wf_adc_idx = param.collate_equal.ref;
+  cmd = waveform.param_analysis.analysis.cmd{param.collate_equal.cmd_idx};
 
   % Taper off end of record to reduce circular convolution effects that may
   % show up during time delay compensation.

cresis-toolbox/tracker/block_data.m

@@ -4,25 +4,71 @@
 % param: Parameter structure from read_param_xls parameter spreadsheet
 %
 % param.block_data: Structure which controls the size of each block
+%
 %  .block_size: number of columns in each block
+%
 %  .block_overlap: the percentage of overlap between each block
-%  .top_gap: number of rows before the first layer
-%  .bottom_pad : number of rows after the deepest layer
-%  .surface_flat_en:	Enable/Disable surface flattening
-%  .surface_rel_layers_flat_en:	Optional feature when surface filtering is enabled. Enable this feature to flatten the layers relative to the filtered surface.
-%  .surface_filter_len:	Specifies the length of the filter for filtering the surface
-%  .pre_detrend_filter_en:	Enable/Disable filtering before detrending
-%  .post_detrend_filter_en:	Enable/Disable filtering after detrending (before normalization)
-%  .uncompress_en:	Depending on the echogram data product used (e.g qlook, post), the echogram may be compressed. This flag when true uncompresses the compressed data using uncompress_echogram function prior to any processing.
-%  .early_trunc:	Truncate data immediately after surface flattening (before detrending and normalizing)
-%  .late_trunc:	Truncate data after all data manipulation( i.e detrending and normalizing ) is done.
+%
 %  .debug_plot:	Set to true for debug plots.
+%
 %  .detrend_debug:	Set to true for detrend debug plots.
-%  .echo_path:	Path to echogram data, typically an argument of ct_filename_out function e.g 'CSARP\standard' => ct_filename_out(param,'CSARP\standard').
-%  .out_fn:	Path where output blocks and files are saved. Currently, this is passed as an argument to ct_filename_tmp to save the outputs in KU user's scratch
-%  .layers_source:	This specifies where the layer data is loaded from(e.g layerdata, records, lidar, etc). This forms a field of the layer_params struct passed into opsLoadLayers. See runOpsLoadLayers.m
-%  .layerdata_source:	When layers_source is layerdata, this string specifies the layerdata (e.g layer_koenig, layer, post) to be loaded. This field is also one of the fields of the layer_params struct passed into opsLoadLayers. See runOpsLoadLayers.m
-%  .regexp:	When layers_source is layerdata, all the layers with layer names that match this regular expression pattern are loaded. This field is also one of the fields of the layer_params struct passed into opsLoadLayers. See runOpsLoadLayers.m
+%
+%  .bottom_pad : number of rows after the deepest layer
+%
+%  .early_trunc:	Truncate data immediately after surface flattening
+%  (before detrending and normalizing)
+%
+%  .echo_path:	Path to echogram data, typically an argument of
+%  ct_filename_out function e.g 'CSARP\standard' =>
+%  ct_filename_out(param,'CSARP\standard').
+%
+%  .late_trunc:	Truncate data after all data manipulation( i.e detrending
+%  and normalizing ) is done.
+%
+%  .layer_params: opsLoadLayers.m structure describing which layers to load
+%  and store in the block files
+%
+%  .layers_source:	This specifies where the layer data is loaded from(e.g
+%  layerdata, records, lidar, etc). This forms a field of the layer_params
+%  struct passed into opsLoadLayers. See runOpsLoadLayers.m
+%
+%  .layerdata_source:	When layers_source is layerdata, this string
+%  specifies the layerdata (e.g layer_koenig, layer, post) to be loaded.
+%  This field is also one of the fields of the layer_params struct passed
+%  into opsLoadLayers. See runOpsLoadLayers.m
+%
+%  .out_fn:	Path where output blocks and files are saved. Currently, this
+%  is passed as an argument to ct_filename_tmp to save the outputs in KU
+%  user's scratch
+%
+%  .post_detrend_filter_en:	Enable/Disable filtering after detrending
+%  (before normalization)
+%
+%  .pre_detrend_filter_en:	Enable/Disable filtering before detrending
+%
+%  .regexp:	When layers_source is layerdata, all the layers with layer
+%  names that match this regular expression pattern are loaded. This field
+%  is also one of the fields of the layer_params struct passed into
+%  opsLoadLayers. See runOpsLoadLayers.m
+%
+%  .surf_param: opsLoadLayers.m structure describing which surface to load
+%  and store in the preprocessing and in the block files
+%
+%  .surface_flat_en:	Enable/Disable surface flattening
+%
+%  .surface_rel_layers_flat_en:	Optional feature when surface filtering is
+%  enabled. Enable this feature to flatten the layers relative to the
+%  filtered surface.
+%
+%  .surface_filter_len:	Specifies the length of the filter for filtering
+%  the surface
+%
+%  .top_gap: number of rows before the first layer
+%
+%  .uncompress_en:	Depending on the echogram data product used (e.g qlook,
+%  post), the echogram may be compressed. This flag when true uncompresses
+%  the compressed data using uncompress_echogram function prior to any
+%  processing.
 %
 % Authors: Ibikunle ( Adapted from John Paden's koenig_mat_loader )
 %
@@ -55,6 +101,10 @@
 %% Input Checks: block_data
 % =========================================================================
 
+if ~isfield(param,'block_data') || isempty(param.block_data)
+  param.block_data = [];
+end
+
 % block_data.block_along_track: scalar double, along-track length of each
 % block in meters, default is 5000 m
 if ~isfield(param.block_data,'block_along_track') || isempty(param.block_data.block_along_track)
@@ -102,6 +152,17 @@
   param.block_data.file.mat_en = true;
 end
 
+% block_data.flatten: structure controlling echo_flatten.m
+if ~isfield(param.block_data,'flatten') || isempty(param.block_data.flatten)
+  param.block_data.flatten = [];
+end
+if ~isfield(param.block_data.flatten,'resample_field') || isempty(param.block_data.flatten.resample_field)
+  param.block_data.flatten.resample_field = [];
+end
+if ~isfield(param.block_data.flatten,'interp_method') || isempty(param.block_data.flatten.interp_method)
+  param.block_data.flatten.interp_method = [];
+end
+
 % Incoherent decimation (inc_dec, inc_B_filter) input check
 % Setting inc_dec = 0: returns coherent data
 % Setting inc_dec = 1: returns power detected data with no decimation
@@ -197,6 +258,12 @@
 ops_param.cmd.frms = frm_list;
 [layers,layer_params] = opsLoadLayers(ops_param, param.block_data.layer_params);
 
+%% Load surface layer
+% =========================================================================
+ops_param = param;
+ops_param.cmd.frms = frm_list;
+[surf,surf_param] = opsLoadLayers(ops_param, param.block_data.surf_param);
+
 %% Block Loop
 % =========================================================================
 frm_mask = false(size(frames.frame_idxs));
@@ -239,13 +306,12 @@
     end
   end
   Nx_original = length(cat_data.GPS_time);
+  if ~isempty(param.block_data.surf_param)
+    cat_data.Surface = interp_finite(interp1(surf.gps_time,surf.twtt,cat_data.GPS_time),0);
+  end
 
   %% Block: Echogram layer flattening
   % =======================================================================
-  param.block_data.flatten.resample_field = [];
-  param.block_data.flatten.resample_field.name = 'surface';
-  param.block_data.flatten.resample_field.source = 'layerdata';
-  param.block_data.flatten.interp_method = [];
   if isfield(param.block_data,'flatten') && ~isempty(param.block_data.flatten)
     [cat_data.Data,resample_field] = echo_flatten(cat_data, ...
       param.block_data.flatten.resample_field, false, ...
@@ -316,9 +382,28 @@
   cat_data.Surface = interp1(dec_idxs,cat_data.Surface,new_axis);
   resample_field = interp1(dec_idxs.',resample_field.',new_axis).';
 
-  %% Block: Extract layers
+  %% Block: Extract surface
   % =======================================================================
   Nx = length(cat_data.GPS_time);
+  if isempty(param.block_data.surf_param)
+    surf_bin = interp1(cat_data.Time,1:length(cat_data.Time),cat_data.Surface);
+  else
+    % Loaded an update to the surface
+    cat_data.Surface = interp_finite(interp1(surf.gps_time,surf.twtt,cat_data.GPS_time),0);
+    twtt = interp1(surf.gps_time,surf.twtt,cat_data.GPS_time);
+    surf_bin = zeros(size(twtt));
+    if isfield(param.block_data,'flatten') && ~isempty(param.block_data.flatten)
+      for rline = 1:Nx
+        % time_flat: vector of twtt associated with each pixel for the
+        % particular column
+        time_flat = interp1(1:length(cat_data.Time),cat_data.Time,resample_field(:,rline),'linear','extrap');
+        surf_bin(rline) = interp1(time_flat,1:length(time_flat),twtt(rline),'linear','extrap');
+      end
+    end
+  end
+
+  %% Block: Extract layers
+  % =======================================================================
   layers_bin_bitmap = zeros(size(cat_data.Data),'uint8');
   layers_bitmap = zeros(size(cat_data.Data),'uint16');
   layers_segment_bitmap = zeros(size(cat_data.Data),'uint16');
@@ -370,13 +455,10 @@
 
   %% Block: Detrend
   % =======================================================================
-  param.block_data.detrend
-  param.block_data.norm
   if isfield(param.block_data,'detrend') && ~isempty(param.block_data.detrend)
-    echo_detrend_data.Data = cat_data.Data;
-    echo_detrend_data.Time = cat_data.Time;
-    echo_detrend_data.Surface = cat_data.Surface;
-    data = echo_detrend(echo_detrend_data,param.block_data.detrend);
+    param.block_data.detrend.layer_top = surf_bin;
+    param.block_data.detrend.layer_bottom = nan(size(surf_bin));
+    cat_data.Data = echo_detrend(cat_data,param.block_data.detrend);
   end
 
   %% Block: TBD (roll compensation, etc.)

cresis-toolbox/tracker/layer_tracker_input_check.m

@@ -320,6 +320,19 @@
     error('Unsupported max diff method %s. Options are merge_vectors, interp_finite. The default is interp_finite unless a dem or reference layer is provided.', track.init.max_diff_method);
   end
 
+  if ~isfield(track,'lsm') || isempty(track.lsm)
+    track.lsm = [];
+  end
+  % lsm.use_mean_surf_en: logical scalar, default false. If true, tracker
+  % uses the mean surface to initialize lsm.y
+  if ~isfield(track.lsm,'use_mean_surf_en') || isempty(track.lsm.use_mean_surf_en)
+    track.lsm.use_mean_surf_en = false;
+  end
+  % lsm.y: initial surface layer bin
+  if ~isfield(track.lsm,'y') || isempty(track.lsm.y)
+    track.lsm.y = 1;
+  end
+
   if ~isfield(track,'max_bin') || isempty(track.max_bin)
     track.max_bin = inf;
   elseif isstruct(track.max_bin)

cresis-toolbox/tracker/layer_tracker_task.m

@@ -675,7 +675,7 @@
   track.max_bin = track.max_bin - min_min_bin;
   track.crossovers(2,:) = track.crossovers(2,:) - min_min_bin;
 
-  %% Track: Create Initial Surface
+  %% Track: Create Initial Layer
   if strcmpi(track.init.method,'dem')
     % Correct for min_bin removal
     track.dem = track.dem - min_min_bin + 1;
@@ -725,15 +725,10 @@
   end
 
   %% Track: Tracking
-  if strcmpi(track.method, 'lsm')
-
-    if track.flag == 1
-      surf = interp_finite(interp1(param.layer_tracker.gt_params.gps_time,param.layer_tracker.gt_params.twtt,mdata.GPS_time));
-      surf_bins = round(interp1(mdata.Time,1:length(mdata.Time),surf));
-      if track.lsm.y == 1
-        track.lsm.y = mean(surf_bins);
-      end
-    end
+  if strcmpi(track.method, 'lsm') && track.lsm.use_mean_surf_en
+    surf = interp_finite(interp1(param.layer_tracker.gt_params.gps_time,param.layer_tracker.gt_params.twtt,mdata.GPS_time));
+    surf_bins = round(interp1(mdata.Time,1:length(mdata.Time),surf));
+    track.lsm.y = mean(surf_bins);
   end
 
   if strcmpi(track.method,'threshold')

cresis-toolbox/tracker/run_block_data.m

@@ -17,18 +17,35 @@
 %params = read_param_xls(ct_filename_param('rds_param_2009_Antarctica_TO.xls'),'20091228_01');
 params = ct_set_params(params,'cmd.generic',0);
 params = ct_set_params(params,'cmd.generic',1,'day_seg','20160519_0[14]'); % 20120330_03
-params = ct_set_params(params,'cmd.frms',[34:500]);
+params = ct_set_params(params,'cmd.frms',[34:44]);
 
 param_override.block_data.block_along_track = 50e3; % Along-track length of each block
 param_override.block_data.block_Nx = 256; % Number of samples in each block
 param_override.block_data.block_overlap = 0.5; % Set the % of overlap between each block
-param_override.block_data.rows.t0_pad = 150;
+param_override.block_data.rows.t0_pad = 50;
 param_override.block_data.rows.t1_pad = 35;
 
+param_override.block_data.inc_B_filter = ones(1,41)/41;
+param_override.block_data.inc_dec = 20;
+param_override.block_data.nan_dec_normalize_threshold = [];
+param_override.block_data.nan_dec = false;
+param_override.block_data.norm.scale = [0.4 1];
+param_override.block_data.detrend.method = 'polynomial';
+
+param_override.block_data.flatten.resample_field = [];
+param_override.block_data.flatten.resample_field.name = 'surface';
+param_override.block_data.flatten.resample_field.source = 'layerdata';
+param_override.block_data.flatten.interp_method = [];
+
 param_override.block_data.echo_path = 'CSARP_post/qlook';
 
 param_override.block_data.out_path = 'block_data';
 
+param_override.block_data.surf_param = [];
+param_override.block_data.surf_param.name = 'surface';
+param_override.block_data.surf_param.source = 'layerdata';
+param_override.block_data.surf_param.layerdata_source = 'layer_overly2021';
+
 param_override.block_data.layer_params = [];
 param_override.block_data.layer_params.name = 'surface';
 param_override.block_data.layer_params(1).source = 'layerdata';