vBar (OCM) Updates

README.md

@@ -1,2 +1,16 @@
 # Video-Currents-Toolbox
-This repository contains code and documentation to quantify surf-zone currents from video imagery.  It is under development and has not had an official release yet.
+This repository contains code and documentation to quantify longshore surf-zone currents from video imagery.
+
+The code is written in MATLAB, and requires version 2023b
+
+Based on the original work of: Chickadel, C.C., Holman, R.A. and Freilich, M.H., 2003. An optical technique for the measurement of longshore currents. Journal of Geophysical Research: Oceans, 108(C11).
+
+Toolbox Input: This code has been written for vbar pixel arrays with naming convention: 1506873540.Sun.Oct.01_15_59_00.GMT.2017.argus02b.cx.vbar125.mat, for example. Input file data should contain (1) pixel data (2) time (3) camera numbers (4) xyz position. 
+Toolbox Output: This code will generate a table for each vbar transect listed in params.transects, containing 'x', 'y', 'vC', and 'wV'. 'vC' is a  structure variable that contains the timeseries analysis variables generated for the window, and 'wV' is the representative weighted velocity over the timeseries.  
+More detail on data requirements and output structure is provided in the user manual document. 
+
+The "videoCurrentsDemoNew.m" code will run you through the program starting from parameter selection ("vidCurrentsParams.m" file). 
+
+If longshore current direction is unknown, params.vBounds = [], and a direction will be interpreted from the data using the radonVbarDir.m function.
+
+"tcolor.m", "wmean.m", "redblue.m" are support functions. 

Support Functions/tcolor.m

@@ -0,0 +1,74 @@
+function H = tcolor(x,y,c,varargin)
+% for use like pcolor but with c as a true color matrix (uses texturemap for speed)...
+%
+% H=tcolor(x,y,c[,method])
+%
+% valid methods are: 'corners','normal','triangles'
+%
+% The normal method is texture mapping unto the plane given by x and y 
+% (which may be distorted arbitrarily)
+%
+% The 'corners' method is the fastest way to draw. However it requires that
+% the area is non-distorted... E.g. that the box defined by the corners
+% defines the area.
+%
+% The slowest method is 'triangles'... (sort of like pcolor). But shading
+% interp works with it.
+%
+% c=imread('C:\Projects\My Pictures\peppermint_girl.jpg');
+% [x,y] = meshgrid(1:size(im,2),1:size(im,1));
+% x=x+y/10; %skew the image
+% H=tcolor(x,y,c,'corners')
+%
+% Aslak Grinsted - July 2003
+
+cax = newplot;
+hold_state = ishold;
+
+lims = [min(min(x)) max(max(x)) min(min(y)) max(max(y))];
+
+
+method=1; %1=normal,2=corners,3=triangles
+if length(varargin)>0
+    method=strmatch(lower(varargin{1}),strvcat({'normal' 'corners' 'triangles'}));
+end
+
+triangles=0;
+switch method
+    case 1
+    case 2
+        if length(size(x))==2
+            x=x([1 end],[1 end]);
+            y=y([1 end],[1 end]);
+        else
+            x=x([1 1;end end]);
+            y=y([1 end;1 end]);    
+        end
+    case 3
+        triangles=1;
+    otherwise
+        error('Unknown method?')
+        return
+end
+
+if length(size(x))==1
+    [x,y]=meshgrid(x,y);
+end
+
+if triangles
+    H=patch(surf2patch(x,y,zeros(size(x)),im2double(c),'triangles'),'edgecolor','none','facecolor','flat');
+else
+    H=surface(x,y,zeros(size(x)),c,'EdgeColor','none','FaceColor','texturemap');
+end
+
+
+if ~hold_state
+    set(cax,'View',[0 90]);
+    set(cax,'Box','on');
+    axis(lims);
+end
+
+
+if (nargout==0) clear H; end;
+
+% Aslak Grinsted (2024). tcolor (a fast pcolor that likes RGB images) (https://www.mathworks.com/matlabcentral/fileexchange/3777-tcolor-a-fast-pcolor-that-likes-rgb-images), MATLAB Central File Exchange.

Support Functions/wmean.m

@@ -0,0 +1,49 @@
+function y = wmean(x,w,dim)
+%WMEAN   Weighted Average or mean value.
+%   For vectors, WMEAN(X,W) is the weighted mean value of the elements in X
+%   using non-negative weights W. For matrices, WMEAN(X,W) is a row vector 
+%   containing the weighted mean value of each column.  For N-D arrays, 
+%   WMEAN(X,W) is the weighted mean value of the elements along the first 
+%   non-singleton dimension of X.
+%
+%   Each element of X requires a corresponding weight, and hence the size 
+%   of W must match that of X.
+%
+%   WMEAN(X,W,DIM) takes the weighted mean along the dimension DIM of X. 
+%
+%   Class support for inputs X and W:
+%      float: double, single
+%
+%   Example:
+%       x = rand(5,2);
+%       w = rand(5,2);
+%       wmean(x,w)
+
+if nargin<2
+    error('Not enough input arguments.');
+end
+
+% Check that dimensions of X match those of W.
+if(~isequal(size(x), size(w)))
+    error('Inputs x and w must be the same size.');
+end
+
+% Check that all of W are non-negative.
+if (any(w(:)<0))
+    error('All weights, W, must be non-negative.');
+end
+
+% Check that there is at least one non-zero weight.
+if (all(w(:)==0))
+    error('At least one weight must be non-zero.');
+end
+
+if nargin==2, 
+  % Determine which dimension SUM will use
+  dim = min(find(size(x)~=1));
+  if isempty(dim), dim = 1; end
+end
+
+y = sum(w.*x,dim)./sum(w,dim);
+
+% Adam Auton (2024). wmean (https://www.mathworks.com/matlabcentral/fileexchange/14416-wmean), MATLAB Central File Exchange.

loadVbarRawFile.m

@@ -0,0 +1,84 @@
+function [sampleStack, timex] = loadVbarRawFile(fileSearchPath, searchDate, transects)
+% Function to load vBar and timex files for specified transects on a specific date
+%
+% fileSearchPath: master directory to search
+% searchDate: the date to match (in datetime format)
+% transects: array of x-locations to match (e.g., [125, 150, 175, 200, 225])
+
+% Convert the target date to the required format (e.g., 'Tue.Oct.17')
+targetDateStr = datestr(searchDate, 'ddd.mmm.dd');
+targetTimeStr1 = datestr(searchDate, 'HH_MM');  % e.g., '11_59'
+targetTimeStr2 = datestr(searchDate-minutes(1), 'HH_MM');
+
+% Generate the path for all subfolders
+searchPath = genpath(fileSearchPath);
+
+% Split the search path into individual folders
+folders = strsplit(searchPath, pathsep);
+
+% Initialize an empty array to hold the .mat file information
+matFileNames = [];
+
+% Loop through each folder and search for .mat files
+for i = 1:length(folders)
+    folder = folders{i};
+    if ~isempty(folder)
+        % Get the .mat files in the current folder
+        files = dir(fullfile(folder, '*.mat'));
+        % Append the files to the matFileNames array
+        matFileNames = [matFileNames; files];
+    end
+end
+
+% Initialize variables to store the files for vbar and timex
+vbarFiles = cell(length(transects), 1);  % Initialize cell array for each transect
+timexFiles = [];
+
+for i = 1:length(matFileNames)
+    fileName = matFileNames(i).name;
+    if contains(fileName, targetDateStr) && contains(fileName, targetTimeStr1) && ...
+            contains(fileName, '.timex.merge')
+        timexFiles = [timexFiles; matFileNames(i)];
+    end
+
+    % Loop through the transects
+    for j = 1:length(transects)
+        transect = transects(j);
+        transectStr = sprintf('.vbar%d', transect);  % Create the string for the current transect (e.g., .vbar125)
+        % Check if the file name matches the desired pattern for vbar files
+        if contains(fileName, targetDateStr) && contains(fileName, targetTimeStr2) && ...
+                contains(fileName, 'vbar') && ...
+                contains(fileName, transectStr)  % Check the current transect
+            vbarFiles{j} = [vbarFiles{j}; matFileNames(i)];  % Store vbar files for each transect
+        end
+    end
+end
+
+% Load the specified variables from the filtered files
+% There should only be one timex
+if ~isempty(timexFiles)
+    timex = load(fullfile(timexFiles.folder, timexFiles.name), 'x', 'y', 'Ip');
+else
+    sprintf('No timex file found.');
+end
+
+% Initialize sampleStack for holding vBar data
+sampleStack = cell(length(transects), 1);
+
+% Loop through vBar files for each transect and load data into struct array
+for j = 1:length(transects)
+    for i = 1:length(vbarFiles{j})
+        filePath = fullfile(vbarFiles{j}(i).folder, vbarFiles{j}(i).name);
+        disp(['Loading file: ', filePath]);
+        % Load the variables and store them in the sampleStack struct
+        loadedData = load(filePath, 'XYZ', 'T', 'RAW', 'CAM');
+
+        % Assign the loaded data to the sampleStack struct at the j-th position
+        sampleStack{j}.XYZ = loadedData.XYZ;
+        sampleStack{j}.T = loadedData.T;
+        sampleStack{j}.RAW = loadedData.RAW;
+        sampleStack{j}.CAM = loadedData.CAM;
+    end
+end
+
+end

prepDataForInput.m

@@ -0,0 +1,40 @@
+function [inpDat] = prepDataForInput(sampleStack,params)
+
+% Grab stack: grayscale image of timestack, size [MxN]
+data = sampleStack.RAW;
+
+% Grab time: time line (starting from zero) of stack, size [1xM]
+time = sampleStack.T;
+
+% Grab xy: x,y position [Nx2] of each pixel in a dimension (equally spaced)
+xyz = sampleStack.XYZ;
+
+cam = sampleStack.CAM; 
+
+% create "list" of cameras.
+for j = 1:max(double(sampleStack.CAM))
+    %For each camera, find extent of XY and RAW. Store variable
+    fieldNameI = sprintf('cam%1.0d', j);
+    camList.(fieldNameI).XY = xyz(cam == j, 1:2);
+    camList.(fieldNameI).RAW = data(:,cam == j);
+    % temp vars for iteration of loop
+    tempXY = camList.(fieldNameI).XY;
+    tempRAW = camList.(fieldNameI).RAW;
+
+    minY = max(ceil(min(tempXY(:,2))), min(params.yLims));
+    maxY = min(floor(max(tempXY(:,2))), max(params.yLims));
+
+    % Use Y limits of camera (and T on x-axis) to create grid
+    yBounds = minY : params.delY : maxY;
+    [tGrid, yGrid] = meshgrid(time, yBounds);
+
+    % Define a new variable called "inputData" which stores the gridded
+    % data that will be used as input for "videoCurrentGen"
+    inpDat.(fieldNameI).tGrid = tGrid;
+    inpDat.(fieldNameI).yGrid = yGrid;
+
+    % Project 'RAW' data for this camera onto the grid
+    inpDat.(fieldNameI).rawGrid = griddata(tempXY(:,2), time, double(tempRAW), yGrid, tGrid, 'natural');
+end
+
+end

radonVbarDir.m

@@ -0,0 +1,64 @@
+function [params] = radonVbarDir(inpDat, params, plotFlag)
+%% Apply Radon to define the velocity search bounds (vBounds) to + or -
+theta = -85:85;  % the range of angles we're considering
+
+% Check if plotFlag is provided, if not, set it to 0 (no plotting)
+if nargin < 2
+    plotFlag = 0;
+end
+
+% Convert time to numeric 'datenum' format for plotting
+numeric_time = inpDat.tGrid(1,:) / (3600.0 * 24) + datenum(1970,1,1);
+plot_time = datetime(numeric_time, 'ConvertFrom', 'datenum');  % Use datetime for better interpretation
+
+% Perform Radon transform on all angles between -85 and +85 degrees
+[R, xp] = radon(inpDat.rawGrid, theta);
+
+% Exclude energy from stationary objects (~0 degrees = shore normal)
+excluded_angles = (theta >= -5 & theta <= 5);  % Exclude shore-normal angles
+R(:, excluded_angles) = 0;
+
+% Find peaks in the Radon transform to identify the dominant angle
+
+maxR = max(R, [], 1);
+[peaks, locations] = islocalmax(maxR);
+
+% Plot the results if plotFlag is set
+if plotFlag
+    figure();
+    % Plot the intensity data (raw pixel data)
+    pcolor(numeric_time, inpDat.yGrid(:,1), inpDat.rawGrid);
+    hold on;
+    shading flat;
+    colormap gray;
+    title('Detected Oblique Foam Traces');
+    xlabel('Time');
+    ylabel('Y-Position (m)');
+    datetick('x', 'HH:MM:SS', 'keeplimits');  % Convert numeric time back to readable format
+    hold off;
+
+    figure();
+    plot(theta, maxR, 'LineWidth', 2);
+    hold on;
+    scatter(theta(peaks), maxR(peaks), 'r*');  % Plot detected peaks
+    xlabel('Angle (degrees)');
+    ylabel('Max Radon Energy');
+    title('Radon Transform - Detected Peaks');
+    grid on;
+end
+
+% Detect the angle of the highest peak
+%detected_angle = theta(locations(peaks == max(peaks)));  % Angle of the highest peak
+neg_angle_energy = sum(maxR(theta >= -60 & theta >= -5));  % Sum energy over all negative angles
+pos_angle_energy = sum(maxR(theta >= 5 & theta >= 60));  % Sum energy over all positive angles
+
+% Interpretation of the detected angle for setting velocity bounds (vBounds)
+if neg_angle_energy > pos_angle_energy
+    disp('The foam is drifting south');
+    params.vBounds = [-2.5 -0.01];  % Adjust bounds for southward motion
+elseif neg_angle_energy < pos_angle_energy
+    disp('The foam is drifting north');
+    params.vBounds = [0.01 2.5];  % Adjust bounds for northward motion
+else
+    disp('The dominant angle is within the excluded range.');
+end

vcTableGen.m

@@ -0,0 +1,45 @@
+function [vcTable] = vcTableGen(inputDataX, params, xLocation)
+
+% Step 1: Initialize the table with y-values and cam-values
+%vcTable = table([], struct([]), [], 'VariableNames', {'y', 'vC', 'wV'});
+vcTable = table(); 
+count = 1; 
+
+for j = 1:params.numCams
+    fieldNameJ = sprintf('cam%1.0d', j);
+
+    minY = min(inputDataX.(fieldNameJ).yGrid, [],  'all'); 
+    maxY = max(inputDataX.(fieldNameJ).yGrid, [], 'all'); 
+
+    % Define the center points of the analysis windows
+    inputDataX.(fieldNameJ).yCenters = (minY+params.tileSize/2):params.tileSize:(maxY-params.tileSize/2);
+
+    for k = 1:length(inputDataX.(fieldNameJ).yCenters)
+        y = inputDataX.(fieldNameJ).yGrid(:,1); 
+
+        y1 = inputDataX.(fieldNameJ).yCenters(k) - params.tileSize/2;    % start point 
+        y2 = inputDataX.(fieldNameJ).yCenters(k) + params.tileSize/2;    % end point
+
+        i1 = find(y == y1, 1, 'first'); 
+        i2 = find(y == y2, 1, 'first'); 
+
+        % Run video-current-toolbox
+        % stack, time, xy, vBounds, fkBounds, Twin, Tstep {plotFlag})
+        stack = inputDataX.(fieldNameJ).rawGrid(i1:i2,:)'; 
+        xy = inputDataX.(fieldNameJ).yGrid(i1:i2,1); 
+        vC = videoCurrentGen(stack, params.mtime, xy, ...
+                params.vBounds, params.fkBounds, params.tWindow, params.tStep, params.plotFlag);
+
+        % Save vC to the table
+        vcTable.x(count) = xLocation; 
+        vcTable.y(count) = inputDataX.(fieldNameJ).yCenters(k);       % midpoint
+        vcTable.vC{count} = vC; 
+        indexKeep = find((vcTable.vC{count}.SNR > 5)&(vcTable.vC{count}.prob>0.8)&(vcTable.vC{count}.stdI>5));
+        vcTable.wV(count) = wmean(vC.meanV(indexKeep), 1./vC.stdV(indexKeep), 'omitnan');
+        vcTable.camNum(count) = j; 
+        count = count+1; 
+    end
+end
+
+vcTable = sortrows(vcTable, 'y');   
+end