Algorithme d'Inpainting en utilisant la méthode du gradient conjugué

Bonjour à tous,
Je travaille actuellement sur un algorithme d'inpainting: l'entrée est une image, j'applique un masque (disposition aléatoire de pixels noirs), puis l'algorithme tente de récupérer l'image originale en utilisant la méthode de gradient conjugué. Je l'ai testé pour les images en niveaux de gris, il fonctionne sans problème. Je suis maintenant en train de l'implémenter pour les images colorées. Je me sers de Matlab. Une image colorée est reconnue par Matlab comme un tableau de dimension m x n x 3 (m est la hauteur, et n est la largeur). Voici le script que j'ai écrit:
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 clear all; tic; %%%%%%%%%%%%%%Lit l'image d'entrée I = imread('Google.png'); %ici, une icône 256x256, en uint8 I = double(I); %On la convertit en double [ N1,N2 ] = size(I(:,:,1)); %On stocke les dimensions de l'image %%%%%%%%%%%%%%%%%%Pour créer le masque Masque = zeros(size(I));%Mêmes dimensions que l'image d'entrée Masque(:,:,1) = max(min(round(0.6*rand(size(I(:,:,1)))),1),0);%Disposition aléatoire de 0 & 1 comme coefficients Masque(:,:,2) = Masque(:,:,1); %Je pose comme contrainte d'avoir le même masque pour les 3 matrices stockées dans I Masque(:,:,3) = Masque(:,:,1); Masque = double(Masque); Ic = I.*Masque; %On crée l'image corrompue en multipliant par le masque %%%%%%%%%%%%%%%%%%%% Méthode du gradient conjugué Ifix = zeros(size(I)); %Pour stocker les résultats for i=1:3 tic; disp('Étape : ');disp(i); res = 1; %idée de "précision" pour l'algorithme b = Ic(:,:,i); xk = b; %xk est une suite récurrente, avec comme 1er élément l'image corrompue, et qui va converger vers le résultat : l'image "réparée" dk = b - matrice_A(xk,Masque(:,:,i)); while res >10^(-2), %Tant que la précision est trop grand, on effectue les étapes du gradient conjugué rk = b - matrice_A(xk,Masque(:,:,i)); Ark = matrice_A(rk,Masque(:,:,i)); alphak = sum(sum(rk.*rk))/sum(sum(Ark.*rk)); xk = xk + alphak*dk; temp = rk; rk = b - matrice_A(xk,Masque(:,:,i)); betak = sum(sum(rk.*rk))/sum(sum(temp.*temp)); dk = rk + betak*dk; res = dk(:)'*dk(:)/(N1*N2); end Ifix(:,:,i) = xk; disp('Durée de l étape :'); toc; end Ifix = uint8(Ifix);%Je convertis en uint8 pour l'afficher correctement avec la fonction "image" image(Ifix); toc;
La fonction "matrice_A" est la suivante :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 function [ J ] = matrice_A(I,Masque) %En gros, pour multiplier I par Masque, en respectant des contraintes liées à une interpolation spline [ N1,N2 ] = size(I); J = zeros(size(I)); %sortie for i=1:N1, bool_iplus = i<N1; iplus = min(i+1,N1); bool_imoins = i>1; imoins = max(i-1,1); for j=1:N2, bool_jplus = j<N2; jplus = min(j+1,N2); bool_jmoins = j>1; jmoins = max(j-1,1); if(Masque(i,j)>0) J(i,j) = 1/36*(16*I(i,j) + 4*(bool_iplus*I(iplus,j) + bool_imoins*I(imoins,j) + ... bool_jplus*I(i,jplus) + bool_jmoins*I(i,jmoins)) + ... (bool_iplus*bool_jplus*I(iplus,jplus) + bool_imoins*bool_jplus*I(imoins,jplus) + ... bool_imoins*bool_jmoins*I(imoins,jmoins) + bool_iplus*bool_jmoins*I(iplus,jmoins))); else J(i,j) = -(N1*N2)*(-8*I(i,j) +1*(I(iplus,j) + I(imoins,j) + I(i,jplus) + I(i,jmoins)) + ... (I(iplus,jplus) + I(imoins,jplus) + I(imoins,jmoins) + I(iplus,jmoins))); end end end
Ce programme marche, par exemple avec cette image en entrée : Pièce jointe 207647
J'obtiens cette image en résultat : Pièce jointe 207648
Et pour cette entrée : Pièce jointe 207649
On obtient : Pièce jointe 207650

Je pense que vous serez d'accord pour dire que le résultat pourrait être meilleur. Cependant, je ne sais pas quoi changer pour améliorer le rendu. Des idées ?
Aussi, le programme a tourné pendant 60 minutes pour ce résultat, je pense également qu'il pourrait être optimisé (surtout pour de plus grosses images, ici elles sont de taille assez faible).
Avez vous des propositions ? Des choses que je pourrait faire différemment ?
Merci beaucoup !