%% run06c_Himawari_SST_PIV_student.m
% MCC/PIV practice using two Himawari-9 SST images.
%
% This improved version has three modes:
%   PIV_MODE = 'standard' : stable broad-domain vectors
%   PIV_MODE = 'fine'     : smaller windows for coastal/small-eddy vectors
%   PIV_MODE = 'roi_eddy' : high-density PIV around 148E, 38.5N
%
% Sign convention:
%   dx > 0 : pattern moved eastward
%   dy > 0 : pattern moved northward
%
% Important:
%   These are apparent pattern velocities, not perfect ocean currents.
%   SST can change through heating/cooling, mixing, clouds, retrieval errors,
%   and advection.
%
% STUDENT FILL-IN VERSION. Fill the ________ parts.
% Matching input:
%   The MCC/PIV matching image is the original SST field after robust
%   normalization.  A broad high-pass field is NOT subtracted by default.
%   The normalized cross-correlation function removes the mean within each
%   small interrogation window, which is the appropriate local normalization.

% -------------------------------------------------------------------------
% STUDENT VERSION:
%   Complete the blanks marked by ________.
%   This script applies the same MCC/PIV idea to real Himawari SST.
% -------------------------------------------------------------------------

clear; close all; clc;

%% ------------------------------------------------------------------------
% 0. User settings
% -------------------------------------------------------------------------
RUN_ID = 'run06c';
PRODUCT_NAME = 'SST';

% Choose one: 'standard', 'fine', 'roi_eddy'
% For the coastal eddy near 148E, 38.5N, start with 'fine'.
% Use 'roi_eddy' if you want a dense vector field only around that eddy.
PIV_MODE = 'fine';

% Zoom/ROI used for the extra diagnostic figure and for roi_eddy mode.
ROI_LON = [147.0 149.8];
ROI_LAT = [37.2 39.9];
ROI_PAD_DEG = 0.30;

% Save a mode-tagged product and also overwrite the legacy filename used by run06e.
SAVE_LEGACY_FILENAMES = true;

OUTPUT_ROOT = fullfile(pwd, 'out06_Himawari_PIV');
FIG_DIR = fullfile(OUTPUT_ROOT, 'figures');
DATA_DIR = fullfile(OUTPUT_ROOT, 'data');
ensure_dir(OUTPUT_ROOT);
ensure_dir(FIG_DIR);
ensure_dir(DATA_DIR);

%% ------------------------------------------------------------------------
% 1. Load data
% -------------------------------------------------------------------------
file1 = 'HIMAWARI_SST_Hokkaido_20260601_0100_subset.mat';
file2 = 'HIMAWARI_SST_Hokkaido_20260601_0200_subset.mat';

if ~exist(file1,'file')
    error('Cannot find %s. Put this script in the same folder as the MAT files.', file1);
end
if ~exist(file2,'file')
    error('Cannot find %s. Put this script in the same folder as the MAT files.', file2);
end

D1 = load(file1);
D2 = load(file2);

lon = double(D1.lon(:)).';
lat = double(D1.lat(:));
SST1 = double(D1.SST);
SST2 = double(D2.SST);

[lat, SST1] = ensure_lat_ascending(lat, SST1);
[~,   SST2] = ensure_lat_ascending(double(D2.lat(:)), SST2);

SST1(~isfinite(SST1)) = NaN;
SST2(~isfinite(SST2)) = NaN;

% Time difference. The pair is 0100 UTC to 0200 UTC.
dt_hours = 1.0;
dt_seconds = dt_hours * 3600;

%% ------------------------------------------------------------------------
% 2. Analysis settings
% -------------------------------------------------------------------------
switch lower(PIV_MODE)
    case 'standard'
        WINDOW_SIZE = 41;
        STEP = 16;
        MAX_SHIFT = 5;
        MIN_VALID_FRACTION = 0.60;
        MIN_CORRELATION = 0.35;
        MIN_PEAK_GAP = 0.030;
        PEAK_EXCLUDE_RADIUS = 1;
        MIN_TEXTURE_STD = 0.15;
        MAX_SPEED_MS = 1.20;
        NEIGH_DIFF_GRID = 2.5;
        PLOT_ARROW_SCALE = 4.0;
        DO_ROI_CROP = false;

    case 'fine'
        % Smaller interrogation window and vector spacing.
        % This is better for coasts and small eddies than standard mode.
        WINDOW_SIZE = 25;
        STEP = 8;
        MAX_SHIFT = 3;
        MIN_VALID_FRACTION = 0.35;% 0.40;
        MIN_CORRELATION = 0.25; % 0.30;
        MIN_PEAK_GAP = 0.010; %0.015;
        PEAK_EXCLUDE_RADIUS = 1;
        MIN_TEXTURE_STD = 0.035; % 0.06;
        MAX_SPEED_MS = 1.50;
        NEIGH_DIFF_GRID = 2.5; % 2.0;
        PLOT_ARROW_SCALE = 3.0;
        DO_ROI_CROP = false;

    case 'roi_eddy'
        % Dense local PIV around the eddy near 148E, 38.5N.
        WINDOW_SIZE = 21;
        STEP = 5;
        MAX_SHIFT = 3;
        MIN_VALID_FRACTION = 0.35;
        MIN_CORRELATION = 0.28;
        MIN_PEAK_GAP = 0.012;
        PEAK_EXCLUDE_RADIUS = 1;
        MIN_TEXTURE_STD = 0.04;
        MAX_SPEED_MS = 1.50;
        NEIGH_DIFF_GRID = 2.5;
        PLOT_ARROW_SCALE = 2.5;
        DO_ROI_CROP = true;

    otherwise
        error('Unknown PIV_MODE: %s', PIV_MODE);
end

MODE_TAG = lower(PIV_MODE);
PREPROCESS_MODE = 'raw_robust_normalize';
LOCAL_DETREND_MODE = 'window_mean';
% Broad high-pass subtraction is intentionally disabled.
% Each interrogation window is mean-removed inside window_corr().

if DO_ROI_CROP
    lon_full = lon;
    lat_full = lat;
    [lon, lat, SST1] = crop_domain(lon_full, lat_full, SST1, ROI_LON, ROI_LAT, ROI_PAD_DEG);
    [~,   ~,   SST2] = crop_domain(lon_full, lat_full, SST2, ROI_LON, ROI_LAT, ROI_PAD_DEG);
end

dlon = mean(diff(lon));
dlat = mean(diff(lat));

fprintf('=== run06c Himawari SST MCC/PIV (%s mode) ===\n', PIV_MODE);
fprintf('Image size: %d x %d\n', size(SST1,1), size(SST1,2));
fprintf('lon range : %.3f to %.3f E\n', min(lon), max(lon));
fprintf('lat range : %.3f to %.3f N\n', min(lat), max(lat));
fprintf('SST1 range: %.3f to %.3f degC\n', nanmin2(SST1), nanmax2(SST1));
fprintf('SST2 range: %.3f to %.3f degC\n', nanmin2(SST2), nanmax2(SST2));
fprintf('WINDOW_SIZE=%d, STEP=%d, MAX_SHIFT=%d\n', WINDOW_SIZE, STEP, MAX_SHIFT);
fprintf('MCC input : raw SST, robust-normalized; no broad high-pass subtraction\n');
fprintf('Window corr: each interrogation window is mean-removed locally\n');
fprintf('MIN_VALID=%.2f, MIN_R=%.2f, MIN_PEAK_GAP=%.3f\n', MIN_VALID_FRACTION, MIN_CORRELATION, MIN_PEAK_GAP);

%% ------------------------------------------------------------------------
% 3. Make matching images
% -------------------------------------------------------------------------
% Do NOT subtract a broad background field here.
% The goal is to track the actual SST spatial pattern, including fronts and
% eddy-edge structures.  The normalized cross-correlation calculation below
% removes the mean within each interrogation window, so absolute offsets in
% each local window do not dominate the correlation.
I1 = robust_normalize(________, 3.0);
I2 = robust_normalize(________, 3.0);

%% ------------------------------------------------------------------------
% 4. MCC/PIV calculation
% -------------------------------------------------------------------------
P = mcc_piv_core(________, ________, lon, lat, WINDOW_SIZE, STEP, MAX_SHIFT, ...
                 MIN_VALID_FRACTION, PEAK_EXCLUDE_RADIUS);

[P.U, P.V, P.speed] = grid_displacement_to_velocity(________, ________, P.lon, dlon, dlat, dt_seconds);

basic_qc = P.valid_fraction >= ________ & ...
           P.rmax >= ________ & ...
           P.peak_gap >= ________ & ...
           P.texture >= ________ & ...
           P.speed <= ________ & ...
           ~P.edge_flag & ...
           isfinite(P.dx) & isfinite(P.dy);

neighbor_qc = neighbor_consistency_qc(P.dx, P.dy, ________, NEIGH_DIFF_GRID);
final_qc = ________ & ________;

zoom_mask = roi_mask(P, ROI_LON, ROI_LAT);

fprintf('Candidate windows  : %d\n', numel(P.dx));
fprintf('Raw finite vectors : %d\n', nnz(isfinite(P.dx) & isfinite(P.dy))); 
fprintf('Basic QC vectors   : %d\n', nnz(basic_qc));
fprintf('Final QC vectors   : %d\n', nnz(final_qc));
fprintf('Final vectors in ROI %.1f-%.1fE, %.1f-%.1fN: %d\n', ROI_LON(1), ROI_LON(2), ROI_LAT(1), ROI_LAT(2), nnz(final_qc & zoom_mask));
fprintf('Median speed final : %.3f m/s\n', nanmedian2(P.speed(final_qc)));
fprintf('Median peak gap    : %.3f\n', nanmedian2(P.peak_gap(final_qc)));

%% ------------------------------------------------------------------------
% 5. Figures
% -------------------------------------------------------------------------
fig1 = figure('Color','w','Position',[80 80 1500 950]);
tiledlayout(2,2,'TileSpacing','compact','Padding','compact');

nexttile;
plot_map(lon, lat, SST1, [0 30]);
title('Display image 1: SST, 2026-06-01 0100 UTC');
colorbar;

nexttile;
plot_map(lon, lat, SST2, [0 30]);
title('Display image 2: SST, 2026-06-01 0200 UTC');
colorbar;

nexttile;
plot_map(lon, lat, I1, [-3 3]);
title('MCC image 1: raw SST, robust-normalized');
colorbar;

nexttile;
plot_map(lon, lat, I2, [-3 3]);
title('MCC image 2: raw SST, robust-normalized');
colorbar;

sgtitle(sprintf('run06c: Himawari SST images and raw-normalized MCC input (%s mode)', PIV_MODE));

fig2 = figure('Color','w','Position',[80 80 1600 700]);
tiledlayout(1,2,'TileSpacing','compact','Padding','compact');

nexttile;
plot_map(lon, lat, SST1, [0 30]); hold on;
plot_vectors(P, isfinite(P.dx) & isfinite(P.dy), dlon, dlat, PLOT_ARROW_SCALE, 'k');
title('Before QC');
colorbar;

nexttile;
plot_map(lon, lat, SST1, [0 30]); hold on;
plot_vectors(P, final_qc, dlon, dlat, PLOT_ARROW_SCALE, 'k');
title('After final QC');
colorbar;

sgtitle(sprintf('run06c: displacement vectors in lon/lat units (%s mode)', PIV_MODE));

fig3 = figure('Color','w','Position',[80 80 1300 850]);
plot_map(lon, lat, SST1, [0 30]); hold on;
plot_vectors(P, final_qc, dlon, dlat, PLOT_ARROW_SCALE, 'k');
text(min(lon)+0.25, max(lat)-0.35, ...
    sprintf('%s mode: win=%d, step=%d, shift=+/%d px, vectors=%d', PIV_MODE, WINDOW_SIZE, STEP, MAX_SHIFT, nnz(final_qc)), ...
    'FontSize',15,'FontWeight','bold','BackgroundColor','w','Margin',6);
text(min(lon)+0.25, max(lat)-0.85, ...
    sprintf('dt=%.1f h, median speed=%.3f m s^{-1}; arrows x %.1f', dt_hours, nanmedian2(P.speed(final_qc)), PLOT_ARROW_SCALE), ...
    'FontSize',13,'BackgroundColor','w','Margin',6);
title('run06c: apparent SST pattern displacement from MCC/PIV');
colorbar;

fig4 = figure('Color','w','Position',[80 80 1500 850]);
tiledlayout(2,4,'TileSpacing','compact','Padding','compact');

nexttile; plot_piv_field(P, P.valid_fraction, [0 1]); title('Valid fraction'); colorbar;
nexttile; plot_piv_field(P, P.rmax, [0 1]); title('Maximum correlation'); colorbar;
nexttile; plot_piv_field(P, P.peak_gap, [0 0.12]); title('Independent peak gap'); colorbar;
nexttile; plot_piv_field(P, P.texture, []); title('Texture std'); colorbar;
nexttile; plot_piv_field(P, P.speed, [0 MAX_SPEED_MS]); title('Speed before final QC [m s^{-1}]'); colorbar;
nexttile; plot_piv_field(P, double(basic_qc), [0 1]); title('Basic QC'); colorbar;
nexttile; plot_piv_field(P, double(neighbor_qc), [0 1]); title('Neighbor QC'); colorbar;
nexttile; plot_piv_field(P, double(final_qc), [0 1]); title('Final QC'); colorbar;

sgtitle(sprintf('run06c: MCC/PIV quality-control diagnostics (%s mode)', PIV_MODE));

fig5 = figure('Color','w','Position',[80 80 1200 850]);
plot_map(lon, lat, SST1, [0 30]); hold on;
plot_vectors(P, final_qc & zoom_mask, dlon, dlat, PLOT_ARROW_SCALE, 'k');
xlim(ROI_LON); ylim(ROI_LAT);
text(ROI_LON(1)+0.05, ROI_LAT(2)-0.12, ...
    sprintf('ROI zoom: final vectors=%d', nnz(final_qc & zoom_mask)), ...
    'FontSize',15,'FontWeight','bold','BackgroundColor','w','Margin',6);
title('run06c: zoom around the 148E, 38.5N SST eddy');
colorbar;

%% ------------------------------------------------------------------------
% 6. Save figures and result files
% -------------------------------------------------------------------------
fig_bases = { ...
    sprintf('run06c_SST_%s_fig01_input_and_matching_images', MODE_TAG), ...
    sprintf('run06c_SST_%s_fig02_vectors_before_after_QC', MODE_TAG), ...
    sprintf('run06c_SST_%s_fig03_final_vectors_on_SST', MODE_TAG), ...
    sprintf('run06c_SST_%s_fig04_QC_diagnostics', MODE_TAG), ...
    sprintf('run06c_SST_%s_fig05_ROI_zoom_148E_385N', MODE_TAG)};
figure_files = save_figure_set([fig1 fig2 fig3 fig4 fig5], fig_bases, FIG_DIR);

settings = struct();
settings.piv_mode = PIV_MODE;
settings.mode_tag = MODE_TAG;
settings.window_size = WINDOW_SIZE;
settings.step = STEP;
settings.max_shift = MAX_SHIFT;
settings.min_valid_fraction = MIN_VALID_FRACTION;
settings.min_correlation = MIN_CORRELATION;
settings.min_peak_gap = MIN_PEAK_GAP;
settings.peak_exclude_radius = PEAK_EXCLUDE_RADIUS;
settings.min_texture_std = MIN_TEXTURE_STD;
settings.max_speed_ms = MAX_SPEED_MS;
settings.neighbor_diff_grid = NEIGH_DIFF_GRID;
settings.plot_arrow_scale = PLOT_ARROW_SCALE;
settings.preprocess_mode = PREPROCESS_MODE;
settings.highpass_removed = false;
settings.window_mean_removed_in_correlation = true;
settings.local_detrend_mode = LOCAL_DETREND_MODE;
settings.roi_lon = ROI_LON;
settings.roi_lat = ROI_LAT;
settings.roi_pad_deg = ROI_PAD_DEG;
settings.do_roi_crop = DO_ROI_CROP;

piv = build_piv_struct(RUN_ID, PRODUCT_NAME, file1, file2, P, ...
    basic_qc, neighbor_qc, final_qc, dt_hours, dt_seconds, dlon, dlat, settings);

source_name = PRODUCT_NAME;
lon_piv = piv.lon;
lat_piv = piv.lat;
dx_grid = piv.dx_grid;
dy_grid = piv.dy_grid;
U_piv = piv.U;
V_piv = piv.V;
speed_piv = piv.speed;
mask_basic = piv.mask_basic;
mask_neighbor = piv.mask_neighbor;
mask_final = piv.mask_final;
rmax = piv.qc.rmax;
valid_fraction = piv.qc.valid_fraction;
texture = piv.qc.texture;
peak_gap = piv.qc.peak_gap;
peak2 = piv.qc.peak2;
peak_ratio = piv.qc.peak_ratio;
edge_flag = piv.qc.edge_flag;

result_file = fullfile(DATA_DIR, sprintf('run06c_SST_%s_PIV_result_full.mat', MODE_TAG));
vector_file = fullfile(DATA_DIR, sprintf('run06c_SST_%s_PIV_vectors.mat', MODE_TAG));
csv_file    = fullfile(DATA_DIR, sprintf('run06c_SST_%s_PIV_vectors.csv', MODE_TAG));

save(result_file, 'piv','P','basic_qc','neighbor_qc','final_qc','zoom_mask', ...
    'lon','lat','SST1','SST2','I1','I2','settings','figure_files', ...
    'WINDOW_SIZE','STEP','MAX_SHIFT','MIN_VALID_FRACTION','MIN_CORRELATION', ...
    'MIN_PEAK_GAP','PEAK_EXCLUDE_RADIUS','MIN_TEXTURE_STD','MAX_SPEED_MS','NEIGH_DIFF_GRID','PLOT_ARROW_SCALE', ...
    'PREPROCESS_MODE','LOCAL_DETREND_MODE','dt_hours','dt_seconds','dlon','dlat','file1','file2','PIV_MODE','MODE_TAG','ROI_LON','ROI_LAT');

save(vector_file, 'piv','source_name','lon_piv','lat_piv','dx_grid','dy_grid', ...
    'U_piv','V_piv','speed_piv','mask_basic','mask_neighbor','mask_final', ...
    'rmax','valid_fraction','texture','peak_gap','peak2','peak_ratio','edge_flag', ...
    'dt_hours','dt_seconds','dlon','dlat','file1','file2','settings','PIV_MODE','MODE_TAG');

write_piv_csv(piv, csv_file);

if SAVE_LEGACY_FILENAMES
    legacy_result = fullfile(DATA_DIR, 'run06c_SST_PIV_result_full.mat');
    legacy_vector = fullfile(DATA_DIR, 'run06c_SST_PIV_vectors.mat');
    legacy_csv    = fullfile(DATA_DIR, 'run06c_SST_PIV_vectors.csv');
    copyfile(result_file, legacy_result);
    copyfile(vector_file, legacy_vector);
    copyfile(csv_file, legacy_csv);
end

fprintf('Saved full result   : %s\n', result_file);
fprintf('Saved vector product: %s\n', vector_file);
fprintf('Saved vector CSV    : %s\n', csv_file);
fprintf('Saved figures in    : %s\n', FIG_DIR);

function piv = build_piv_struct(run_id, source_name, file1, file2, P, basic_qc, neighbor_qc, final_qc, dt_hours, dt_seconds, dlon, dlat, settings)
    piv = struct();
    piv.run_id = run_id;
    piv.source_name = source_name;
    piv.input_file1 = file1;
    piv.input_file2 = file2;
    piv.dt_hours = dt_hours;
    piv.dt_seconds = dt_seconds;

    piv.lon = P.lon;
    piv.lat = P.lat;
    piv.ix = P.ix;
    piv.iy = P.iy;
    piv.dx_grid = P.dx;
    piv.dy_grid = P.dy;
    piv.U = P.U;
    piv.V = P.V;
    piv.speed = P.speed;

    piv.mask_basic = basic_qc;
    piv.mask_neighbor = neighbor_qc;
    piv.mask_final = final_qc;

    piv.qc = struct();
    piv.qc.rmax = P.rmax;
    piv.qc.valid_fraction = P.valid_fraction;
    piv.qc.texture = P.texture;
    piv.qc.peak_gap = P.peak_gap;
    piv.qc.peak2 = P.peak2;
    piv.qc.peak_ratio = P.peak_ratio;
    piv.qc.edge_flag = P.edge_flag;

    piv.grid = struct();
    piv.grid.dlon = dlon;
    piv.grid.dlat = dlat;
    piv.grid.units = 'degree';
    piv.grid.velocity_units = 'm s^-1';
    piv.grid.displacement_units = 'grid cells';

    piv.settings = settings;
end

function write_piv_csv(piv, csv_file)
    T = table( ...
        piv.lon(:), piv.lat(:), piv.ix(:), piv.iy(:), ...
        piv.dx_grid(:), piv.dy_grid(:), piv.U(:), piv.V(:), piv.speed(:), ...
        piv.mask_basic(:), piv.mask_neighbor(:), piv.mask_final(:), ...
        piv.qc.rmax(:), piv.qc.valid_fraction(:), piv.qc.texture(:), ...
        piv.qc.peak_gap(:), piv.qc.peak2(:), piv.qc.peak_ratio(:), piv.qc.edge_flag(:), ...
        'VariableNames', {'lon','lat','ix','iy','dx_grid','dy_grid','U_ms','V_ms','speed_ms', ...
        'basic_qc','neighbor_qc','final_qc','rmax','valid_fraction','texture','peak_gap','peak2','peak_ratio','edge_flag'});
    writetable(T, csv_file);
end


%% ========================================================================
% Local functions
% ========================================================================

function ensure_dir(dirname)
    if ~exist(dirname, 'dir')
        mkdir(dirname);
    end
end

function figure_files = save_figure_set(figHandles, baseNames, figDir)
    if numel(figHandles) ~= numel(baseNames)
        error('figHandles and baseNames must have the same length.');
    end
    ensure_dir(figDir);
    figure_files = struct('fig', cell(numel(figHandles),1), 'png', cell(numel(figHandles),1));
    for k = 1:numel(figHandles)
        fig_file = fullfile(figDir, [baseNames{k} '.fig']);
        png_file = fullfile(figDir, [baseNames{k} '.png']);
        if ~ishandle(figHandles(k)) || ~isvalid(figHandles(k))
            warning('Figure handle %d is invalid. Skipping figure save.', k);
            figure_files(k).fig = '';
            figure_files(k).png = '';
            continue;
        end
        set(figHandles(k), 'PaperPositionMode', 'auto');
        savefig(figHandles(k), fig_file);
        print(figHandles(k), png_file, '-dpng', '-r200');
        figure_files(k).fig = fig_file;
        figure_files(k).png = png_file;
        fprintf('Saved figure        : %s\n', png_file);
    end
end

function [latOut, Aout] = ensure_lat_ascending(latIn, Ain)
    latOut = latIn(:);
    Aout = Ain;
    if numel(latOut) >= 2 && latOut(1) > latOut(end)
        latOut = flipud(latOut);
        Aout = flipud(Aout);
    end
end

function [lon2, lat2, A2] = crop_domain(lon, lat, A, roiLon, roiLat, padDeg)
    if nargin < 6
        padDeg = 0;
    end
    xok = lon >= roiLon(1)-padDeg & lon <= roiLon(2)+padDeg;
    yok = lat >= roiLat(1)-padDeg & lat <= roiLat(2)+padDeg;
    if nnz(xok) < 10 || nnz(yok) < 10
        error('ROI crop is too small or outside the image. Check ROI_LON and ROI_LAT.');
    end
    lon2 = lon(xok);
    lat2 = lat(yok);
    A2 = A(yok, xok);
end

function B = robust_normalize(A, clipValue)
    x = A(isfinite(A));
    if isempty(x)
        B = A;
        return;
    end
    medx = median(x);
    madx = median(abs(x - medx));
    scale = 1.4826 * madx;
    if ~isfinite(scale) || scale <= 0
        scale = std(x);
    end
    if ~isfinite(scale) || scale <= 0
        scale = 1;
    end
    B = (A - medx) ./ scale;
    B(B >  clipValue) =  clipValue;
    B(B < -clipValue) = -clipValue;
end

function P = mcc_piv_core(I1, I2, lon, lat, win, step, maxShift, minValid, peakExcludeRadius)
    [ny, nx] = size(I1);
    halfWin = floor(win/2);
    margin = halfWin + maxShift + 2;

    iy_list = margin:step:(ny-margin+1);
    ix_list = margin:step:(nx-margin+1);
    [IX, IY] = meshgrid(ix_list, iy_list);

    ncy = numel(iy_list);
    ncx = numel(ix_list);

    dx = nan(ncy, ncx);
    dy = nan(ncy, ncx);
    rmax = nan(ncy, ncx);
    valid_fraction = nan(ncy, ncx);
    texture = nan(ncy, ncx);
    peak_gap = nan(ncy, ncx);
    peak2 = nan(ncy, ncx);
    peak_ratio = nan(ncy, ncx);
    edge_flag = false(ncy, ncx);

    shifts = -maxShift:maxShift;

    for jj = 1:ncy
        iy = iy_list(jj);
        y1 = (iy-halfWin):(iy+halfWin);

        for ii = 1:ncx
            ix = ix_list(ii);
            x1 = (ix-halfWin):(ix+halfWin);
            A = I1(y1, x1);
            texture(jj,ii) = nanstd2(A);

            % Skip windows with very little valid structure in image 1.
            if nnz(isfinite(A)) / numel(A) < minValid
                continue;
            end

            Cmat = nan(numel(shifts), numel(shifts));
            VFmat = nan(numel(shifts), numel(shifts));

            for sy = 1:numel(shifts)
                ddy = shifts(sy);
                y2 = y1 + ddy;
                for sx = 1:numel(shifts)
                    ddx = shifts(sx);
                    x2 = x1 + ddx;
                    B = I2(y2, x2);
                    [cc, vf] = window_corr(A, B);
                    Cmat(sy,sx) = cc;
                    VFmat(sy,sx) = vf;
                end
            end

            [bestR, idx] = max_finite(Cmat);
            if ~isfinite(bestR)
                continue;
            end

            [syBest, sxBest] = ind2sub(size(Cmat), idx);
            dyInt = shifts(syBest);
            dxInt = shifts(sxBest);

            rmax(jj,ii) = bestR;
            valid_fraction(jj,ii) = VFmat(syBest, sxBest);
            edge_flag(jj,ii) = abs(dxInt) == maxShift || abs(dyInt) == maxShift;

            % Independent second peak.  This avoids treating the shoulder of
            % the maximum as a separate candidate.
            secondR = second_peak_excluding_radius(Cmat, syBest, sxBest, peakExcludeRadius);
            peak2(jj,ii) = secondR;
            if isfinite(secondR)
                peak_gap(jj,ii) = bestR - secondR;
                if abs(secondR) > eps
                    peak_ratio(jj,ii) = bestR ./ secondR;
                end
            end

            if valid_fraction(jj,ii) < minValid
                continue;
            end

            % Sub-pixel correction using a 3-point parabolic fit.
            dxSub = dxInt;
            dySub = dyInt;
            if sxBest > 1 && sxBest < size(Cmat,2)
                dxSub = dxInt + parabola_subpixel(Cmat(syBest,sxBest-1), Cmat(syBest,sxBest), Cmat(syBest,sxBest+1));
            end
            if syBest > 1 && syBest < size(Cmat,1)
                dySub = dyInt + parabola_subpixel(Cmat(syBest-1,sxBest), Cmat(syBest,sxBest), Cmat(syBest+1,sxBest));
            end

            dx(jj,ii) = dxSub;
            dy(jj,ii) = dySub;
        end
    end

    P.ix = IX;
    P.iy = IY;
    P.lon = lon(IX);
    P.lat = lat(IY);
    P.dx = dx;
    P.dy = dy;
    P.rmax = rmax;
    P.valid_fraction = valid_fraction;
    P.texture = texture;
    P.peak_gap = peak_gap;
    P.peak2 = peak2;
    P.peak_ratio = peak_ratio;
    P.edge_flag = edge_flag;
end

function [bestR, idx] = max_finite(A)
    v = A(:);
    ok = isfinite(v);
    if ~any(ok)
        bestR = NaN;
        idx = NaN;
        return;
    end
    ii = find(ok);
    [bestR, kk] = max(v(ok));
    idx = ii(kk);
end

function secondR = second_peak_excluding_radius(Cmat, syBest, sxBest, rad)
    C2 = Cmat;
    ny = size(C2,1);
    nx = size(C2,2);
    y1 = max(1, syBest-rad);
    y2 = min(ny, syBest+rad);
    x1 = max(1, sxBest-rad);
    x2 = min(nx, sxBest+rad);
    C2(y1:y2, x1:x2) = NaN;
    [secondR, ~] = max_finite(C2);
end

function [cc, vf] = window_corr(A, B)
    m = isfinite(________) & isfinite(________);
    vf = nnz(m) / numel(A);
    if nnz(m) < 12
        cc = NaN;
        return;
    end
    a = A(m);
    b = B(m);
    a = a - ________(a);
    b = b - ________(b);
    den = sqrt(sum(a.^2) * sum(b.^2));
    if den <= 0 || ~isfinite(den)
        cc = NaN;
    else
        cc = sum(________ .* ________) ./ den;
    end
end

function s = parabola_subpixel(c1, c2, c3)
    if ~isfinite(c1) || ~isfinite(c2) || ~isfinite(c3)
        s = 0;
        return;
    end
    den = c1 - 2*c2 + c3;
    if abs(den) < 1e-12
        s = 0;
    else
        s = 0.5 * (c1 - c3) / den;
        if abs(s) > 1
            s = 0;
        end
    end
end

function [U, V, speed] = grid_displacement_to_velocity(dx, dy, latC, dlon, dlat, dt)
    meters_per_deg_lat = 111000;
    meters_per_deg_lon = 111000 .* cosd(latC);
    U = dx .* dlon .* meters_per_deg_lon ./ dt;
    V = dy .* dlat .* meters_per_deg_lat ./ dt;
    speed = sqrt(U.^2 + V.^2);
end

function keep2 = neighbor_consistency_qc(dx, dy, keep, maxDiffGrid)
    keep2 = keep;
    [ny, nx] = size(dx);
    for j = 1:ny
        for i = 1:nx
            if ~keep(j,i)
                keep2(j,i) = false;
                continue;
            end
            j1 = max(1,j-1); j2 = min(ny,j+1);
            i1 = max(1,i-1); i2 = min(nx,i+1);
            neigh = keep(j1:j2, i1:i2);
            dxn = dx(j1:j2, i1:i2);
            dyn = dy(j1:j2, i1:i2);
            neigh((j-j1+1), (i-i1+1)) = false;
            vv = neigh & isfinite(dxn) & isfinite(dyn);
            if nnz(vv) < 2
                % Keep isolated coastal vectors if they passed all basic QC.
                keep2(j,i) = keep(j,i);
            else
                mdx = median(dxn(vv));
                mdy = median(dyn(vv));
                dd = sqrt((dx(j,i)-mdx).^2 + (dy(j,i)-mdy).^2);
                keep2(j,i) = keep(j,i) && dd <= maxDiffGrid;
            end
        end
    end
end

function plot_map(lon, lat, A, climVals)
    imagesc(lon, lat, A);
    set(gca,'YDir','normal');
    axis tight;
    xlabel('Longitude');
    ylabel('Latitude');
    if ~isempty(climVals)
        caxis(climVals);
    end
end

function plot_vectors(P, mask, dlon, dlat, scaleFactor, colorSpec)
    if nargin < 6
        colorSpec = 'k';
    end
    qx = P.dx .* dlon .* scaleFactor;
    qy = P.dy .* dlat .* scaleFactor;
    quiver(P.lon(mask), P.lat(mask), qx(mask), qy(mask), 0, colorSpec, 'LineWidth', 1.0, 'MaxHeadSize', 1.8);
end

function plot_piv_field(P, F, climVals)
    imagesc(P.lon(1,:), P.lat(:,1), F);
    set(gca,'YDir','normal');
    axis tight;
    xlabel('Longitude');
    ylabel('Latitude');
    if ~isempty(climVals)
        caxis(climVals);
    end
end

function inside = roi_mask(P, roiLon, roiLat)
    inside = P.lon >= roiLon(1) & P.lon <= roiLon(2) & P.lat >= roiLat(1) & P.lat <= roiLat(2);
end

function s = nanstd2(A)
    x = A(isfinite(A));
    if numel(x) <= 1
        s = NaN;
    else
        s = std(x);
    end
end

function m = nanmedian2(A)
    x = A(isfinite(A));
    if isempty(x)
        m = NaN;
    else
        m = median(x);
    end
end

function v = nanmin2(A)
    x = A(isfinite(A));
    if isempty(x)
        v = NaN;
    else
        v = min(x);
    end
end

function v = nanmax2(A)
    x = A(isfinite(A));
    if isempty(x)
        v = NaN;
    else
        v = max(x);
    end
end

function [lon2, lat2, B] = block_nanmean_2d(lon, lat, A, factor)
    if factor == 1
        lon2 = lon;
        lat2 = lat;
        B = A;
        return;
    end
    ny = floor(size(A,1) / factor) * factor;
    nx = floor(size(A,2) / factor) * factor;
    A = A(1:ny, 1:nx);
    lon = lon(1:nx);
    lat = lat(1:ny);

    B = nan(ny/factor, nx/factor);
    for jj = 1:(ny/factor)
        rr = (jj-1)*factor + (1:factor);
        for ii = 1:(nx/factor)
            cc = (ii-1)*factor + (1:factor);
            block = A(rr, cc);
            x = block(isfinite(block));
            if ~isempty(x)
                B(jj,ii) = mean(x);
            end
        end
    end
    lon2 = mean(reshape(lon, factor, []), 1);
    lat2 = mean(reshape(lat, factor, []), 1).';
end
