L'algorithme de Brox est une méthode variationnelle classique pour l'estimation du flux optique dense, introduite par Thomas Brox en 2004. Cette approche repose sur la minimisation d'une fonction d'énergie globale combinant la constance de l'intensité, la constance du gradient et une contrainte de lissage spatial.

classdef BroxSolver
    % BroxSolver : Implémentation de l'algorithme de flux optique variationnel
    
    properties
        poids_lissage   % Paramètre Alpha (lissage)
        poids_gradient  % Paramètre Gamma (gradient)
        poids_couplage  % Paramètre Omega (couplage)
        iterations      % Nombre maximum d'itérations
        epsilon         % Constante de régularisation pour éviter la division par zéro
        niveaux_pyr     % Profondeur de la pyramide d'images
    end
    
    methods
        function obj = BroxSolver(alpha, gamma, omega, iter, pyr)
            % Constructeur de l'estimateur
            if nargin < 1, alpha = 0.001; end
            if nargin < 2, gamma = 0.8;   end
            if nargin < 3, omega = 0.7;   end
            if nargin < 4, iter = 150;    end
            if nargin < 5, pyr = 5;       end
            
            obj.poids_lissage = alpha;
            obj.poids_gradient = gamma;
            obj.poids_couplage = omega;
            obj.iterations = iter;
            obj.niveaux_pyr = pyr;
            obj.epsilon = 1e-9;
        end
        
        function [champ_u, champ_v] = calculer(obj, im1, im2)
            % Normalisation et conversion en niveaux de gris
            if size(im1, 3) == 3, im1 = rgb2gray(im1); end
            if size(im2, 3) == 3, im2 = rgb2gray(im2); end
            
            src = im2double(im1) * 255;
            tgt = im2double(im2) * 255;
            
            % Construction des pyramides gaussiennes
            pyr1 = obj.creerPyramide(src);
            pyr2 = obj.creerPyramide(tgt);
            
            % Initialisation du flux
            [h, w] = size(pyr1{obj.niveaux_pyr});
            u = zeros(h, w);
            v = zeros(h, w);
            
            % Traitement multi-échelle (du grossier au fin)
            for n = obj.niveaux_pyr:-1:1
                curr_im1 = pyr1{n};
                curr_im2 = pyr2{n};
                [H, W] = size(curr_im1);
                
                % Redimensionnement du flux précédent si nécessaire
                if n < obj.niveaux_pyr
                    u = imresize(u, [H, W], 'bilinear') * (W / size(u, 2));
                    v = imresize(v, [H, W], 'bilinear') * (H / size(v, 1));
                end
                
                % Raffinement itératif par warping
                for k = 1:obj.iterations
                    % Gauchissement (Warping) de l'image cible
                    im2_warped = obj.warp(curr_im2, u, v);
                    
                    % Calcul des dérivées spatiales et temporelles
                    [Ix, Iy] = gradient(curr_im1);
                    It = im2_warped - curr_im1;
                    
                    % Estimation des gradients du flux
                    ux = gradient(u); uy = gradient(u, 1);
                    vx = gradient(v); vy = gradient(v, 1);
                    
                    % Calcul de l'énergie de lissage
                    lissage_u = obj.poids_lissage^2 * (ux.^2 + uy.^2);
                    lissage_v = obj.poids_lissage^2 * (vx.^2 + vy.^2);
                    
                    % Terme de gradient et données
                    denom = sqrt(Ix.^2 + Iy.^2 + obj.epsilon);
                    du = -(Ix .* It) ./ (denom + obj.epsilon);
                    dv = -(Iy .* It) ./ (denom + obj.epsilon);
                    
                    % Mise à jour avec régularisation
                    u = u + obj.filtrer(du, 2);
                    v = v + obj.filtrer(dv, 2);
                end
            end
            champ_u = u;
            champ_v = v;
        end
        
        function pyr = creerPyramide(obj, img)
            pyr = cell(1, obj.niveaux_pyr);
            pyr{1} = img;
            for i = 2:obj.niveaux_pyr
                pyr{i} = imresize(pyr{i-1}, 0.5, 'bicubic');
            end
        end
        
        function warped = warp(~, img, u, v)
            [h, w] = size(img);
            [x, y] = meshgrid(1:w, 1:h);
            warped = interp2(x, y, img, x + u, y + v, 'bilinear', 0);
        end
        
        function img_f = filtrer(~, img, sigma)
            f = fspecial('gaussian', ceil(sigma*3)*2+1, sigma);
            img_f = imfilter(img, f, 'replicate');
        end
        
        function afficherFlux(~, u, v)
            mag = sqrt(u.^2 + v.^2);
            ang = atan2(v, u);
            hsv = zeros([size(u), 3]);
            hsv(:,:,1) = (ang + pi) / (2*pi);
            hsv(:,:,2) = 1;
            hsv(:,:,3) = min(mag ./ max(mag(:)), 1);
            imshow(hsv2rgb(hsv));
        end
    end
end

Exemple d'utilisation pratique

Pour utiliser cette classe, il suffit d'instancier l'objet BroxSolver et de lui passer deux images consécutives d'une séquence vidéo.

% Chargement des données
imgA = imread('sequence_01.jpg');
imgB = imread('sequence_02.jpg');

% Initialisation du solveur avec des paramètres personnalisés
% alpha=0.002, gamma=0.75, omega=0.8, iterations=100, niveaux=4
solveur = BroxSolver(0.002, 0.75, 0.8, 100, 4);

% Calcul du flux optique
[u, v] = solveur.calculer(imgA, imgB);

% Visualisation
figure;
solveur.afficherFlux(u, v);
title('Champ de vecteurs - Algorithme de Brox');

Fondements théoriques

Modèle d'énergie

L'algorithme cherche à minimiser une fonctionnelle d'énergie composée de trois termes principaux :

• Terme de données : Il suppose que l'intensité des pixels reste constante entre deux images, mais ajoute également une contrainte sur la constance du gradient d'image pour mieux gérer les changements d'illumination.
• Terme de lissage : Il impose une régularité spatiale au champ de vecteurs (u, v) pour garantir que les vecteurs voisins ont des directions et des amplitudes similaires.
• Terme de couplage : Permet d'unifier les composantes horizontales et verticales du mouvement.

Stratégie Multi-échelle

Le flux optique classique échoue souvent en cas de grands déplacements (supérieurs à 1 pixel). L'implémentation utilise une pyramide gaussienne :

1. L'image est réduite plusieurs fois pour capturer les mouvements globaux à basse résolution.
2. Le flux estimé à une échelle grossière est dilaté et sert de point de départ pour l'échelle supérieure.
3. Le processus de "warping" (gauchissement) compense le mouvement déjà estimé avant de calculer le résidu à chaque itération.

Tableau des paramètres et impacts

Domaines d'application

• Analyse du mouvement : Suivi de trajectoires de véhicules ou de piétons.
• Stabilisation vidéo : Estimation du bougé de la caméra pour correction logicielle.
• Compression : Estimation de vecteurs de mouvement pour les codecs vidéo.
• Imagerie médicale : Suivi de la déformation d'organes entre plusieurs scans.