Covariance Descriptor Fundamentals

Covariance descriptors map visual regions into a compact statistical representation. By capturing first- and second-order pixel statistics, this method encodes spatial layout, intensity distribution, and edge orientation within a single symmetric matrix. The resulting descriptor is invariant to affine transformations of the feature space and demonstrates strong resilience to illumination shifts and partial occlusion.

Feature Space Formulation

Each pixel p within a candidate window is projected into a multi-channel vector combining positional coordinates, radiometric values, and directional gradients:

f(p) = [x, y, I(x,y), |∇<sub>x</sub>|, |∇<sub>y</sub>|]<sup>T</sup>

• x, y: Normalized pixel coordinates relative to the region origin
• I(x,y): Grayscale intensity sampled at the target coordinate
• |∇x|, |∇y|: Magnitudes obtained via Sobel kernel convolution

Matrix Aggregation & Similarity Metric

Given N active pixels in a target zone, the covariance matrix is derived as:

C = (1/(N-1)) * Σ(f_i - μ)(f_i - μ)<sup>T</sup>

where μ represents the empirical mean of the feature ensemble. Because covariance matrices reside on a Riemannian manifold, standard Euclidean subtraction yields misleading results. Instead, we employ a log-Euclidean distance grounded in generalized eigenvalue decomposition:

D(C₁, C₂) = √( Σ ln²(λᵢ) )

The scalars λᵢ are extracted from the pencil problem C₁v = λC₂v. Positive, finite eigenvalues are retained to guarantee numerical stability before applying the logarithmic transformation.

MATLAB Implementation Architecture

The detection workflow integrates proposal generation, descriptor extraction, Riemannian matching, and spatial filtering. The following modules demonstrate a streamlined pipeline optimized for clarity and execution speed.

Pipeline Orchestrator

%% Core Detection Loop
clear; clc; close all;

scene = imread('test_frame.jpg');
is_rgb = size(scene, 3) == 3;
luma = rgb2gray(scene);

figure('Name', 'Covariance Vehicle Tracker', 'NumberTitle', 'off');
subplot(2,2,[1 3]); imshow(scene); title('Original Frame');

% Initialize reference models
car_proto = load('ref_car_desc.mat').matrix;
truck_proto = load('ref_truck_desc.mat').matrix;

ifisempty(car_proto)
    fprintf('Constructing default prototypes from frame center...\n');
    car_proto = build_template(luma, [0.2, 0.2, 0.6, 0.6]);
    truck_proto = build_template(luma, [0.4, 0.1, 0.5, 0.5]);
end

% Generate spatial hypotheses
proposals = assemble_multi_scale_grid(luma);
matched_results = [];
similarity_cutoff = 0.55;

for k = 1:numel(proposals)
    cand = proposals{k};
    raw_roi = luma(cand.r:cand.r+cand.h-1, cand.c:cand.c+cand.w-1);
    
    descriptor = compute_region_stat(raw_roi);
    d_car = evaluate_manifold_dist(car_proto, descriptor);
    d_trk = evaluate_manifold_dist(truck_proto, descriptor);
    
    best_sim = min(d_car, d_trk);
    if best_sim < similarity_cutoff
        category = (d_car < d_trk) ? 'Car' : 'Truck';
        confidence = 1.0 - best_sim;
        matched_results(end+1,:) = [cand.c, cand.r, cand.w, cand.h, confidence, strcmp(category,'Car')];
    end
end

% Enforce spatial exclusivity
final_output = execute_nms(matched_results, 0.35);

% Render overlay
overlay = copy( scene );
for i = 1:size(final_output,1)
    rect = round(final_output(i,1:4));
    cls = final_output(i,6) > 0.5 ? 'Car' : 'Truck';
    color = final_output(i,6) > 0.5 ? [0 1 0] : [1 0 0];
    overlay = insertShape(overlay, 'Rectangle', rect, 'LineWidth', 2, 'Color', color);
    lbl = sprintf('%s: %.2f', cls, final_output(i,5));
    overlay = insertText(overlay, [rect(1)+5 rect(2)-5], lbl, ...
        'FontSize', 11, 'TextColor', 'w', 'BoxOpacity', 0.6);
end

subplot(2,2,2); imshow(overlay); title('Detection Overlay');
subplot(2,2,4); histogram(matched_results(:,5), 'FaceColor',[0.2 0.5 0.8]);
title('Score Distribution'); xlabel('Certainty'); ylabel('Count');
fprintf('Pipeline complete. Total targets: %d\n', size(final_output,1));

Descriptor Extractor

function stat_matrix = compute_region_stat(gray_patch)
    [H, W] = size(gray_patch);
    [Gx, Gy] = gradient(double(gray_patch));
    mag_x = abs(Gx); mag_y = abs(Gy);
    
    [X_coord, Y_coord] = meshgrid(1:W, 1:H);
    feat_stack = [X_coord(:), Y_coord(:), gray_patch(:), mag_x(:), mag_y(:)];
    
    stat_matrix = cov(feat_stack, 1);
    
    % Tikhonov regularization for positive definiteness
    stat_matrix = stat_matrix + eye(5) * 1e-6;
end

Riemannian Distance Calculator

function diff_val = evaluate_manifold_dist(M_ref, M_target)
    try
        [~, E] = eig(M_ref, M_target);
        vals = diag(E);
        
        valid_idx = vals > 0 & isfinite(vals);
        log_specs = log(vals(valid_idx));
        
        if ~isempty(log_specs)
            diff_val = sqrt(sum(log_specs.^2));
        else
            diff_val = Inf;
        end
    catch
        diff_val = norm(M_ref - M_target, 'fro');
    end
end

Proposal Generator

function hyp_list = assemble_multi_scale_grid(img_gray)
    [R, C] = size(img_gray);
    scales = [0.6, 0.8, 1.0, 1.25];
    base_dims = [40, 64, 96];
    
    hyp_list = {};
    idx = 1;
    for sc = scales
        for dim = base_dims
            win = round(dim * sc);
            if win > min(R,C)/1.5 || win < 24, continue; end
            
            stride = max(4, round(win/3));
            for r = 1:stride:(R-win)
                for c = 1:stride:(C-win)
                    hyp_list{idx} = struct('r', r, 'c', c, 'h', win, 'w', win);
                    idx = idx + 1;
                end
            end
        end
    end
    disp(sprintf('Hypotheses generated: %d', length(hyp_list)));
end

Non-Maximum Suppression Routine

function kept = execute_nms(scores_thresh, overlap_limit)
    if isempty(scores_thresh), kept = []; return; end
    
    [~, rank] = sort(scores_thresh(:,5), 'descend');
    sorted = scores_thresh(rank,:);
    kept = [];
    
    while ~isempty(sorted)
        top = sorted(1,:);
        kept = [kept; top'];
        if numel(sorted)==1, break; end
        
        overlaps = zeros(size(sorted,1)-1, 1);
        for j = 2:size(sorted,1)
            b1 = top; b2 = sorted(j,:);
            
            xi1 = max(b1(1), b2(1)); yi1 = max(b1(2), b2(2));
            xi2 = min(b1(1)+b1(3), b2(1)+b2(3)); yi2 = min(b1(2)+b1(4), b2(2)+b2(4));
            
            if xi2>xi1 && yi2>yi1
                inter = (xi2-xi1)*(yi2-yi1);
                union = prod(b1(3:4)) + prod(b2(3:4)) - inter;
                overlaps(j-1) = inter/union;
            end
        end
        
        keeper = find(overlaps < overlap_limit) + 1;
        sorted = sorted([1; keeper], :);
    end
end

Prototype Aggregator

function model_mat = build_template(reference_img, norm_coords)
    [H, W] = size(reference_img);
    cx = floor(norm_coords(1)*W) + 1;
    cy = floor(norm_coords(2)*H) + 1;
    cw = ceil(norm_coords(3)*W);
    ch = ceil(norm_coords(4)*H);
    
    patch = reference_img(max(1,cy):min(H,cy+ch-1), ...
                          max(1,cx):min(W,cx+cw-1));
    model_mat = compute_region_stat(patch);
end

Deployment Configuraton

Computational Profile

Implementation considerations favor hardware-accelerated linear algebra routines for real-time throughput. When deployed alongside deep neural architectures, covariance descriptors excel as lightweight verification filters or embedded system primitives where memory footprints must remain constrained.