How to calculate the Data-Term for Patch Based Inpainting

Hello everyone;

I've been working on inpainting algorithms (without AI or DNN) for a few months now. The aim is to later modify these inpainting algorithms and test them against each other. But with patch-based inpainting I now have the technical problem that I simply don't understand how to determine the data term.

I understand that the data term (Dp) contains the structure values and is therefore of crucial importance for inpainting, as the structures that run into the inpainting region are treated with priority. For this purpose, gradient values or isophote values should be determined directly at the edge pixel and a priority should be determined from these.

I fully understood the purpose of the Dp and how it should be used afterwards. Unfortunately, I don't understand how I can actually calculate the Dp and this is what I'm hoping for concrete help with.

This is my code so far. It may works if the Data Term is not very important fpr the inpainting task as the data term is just constant zero so far to test if the surrounding function is working. The paper is used for creating this is found in the first line. As I comment my code in German only I added translations in '()' behind them.

%Basisquelle: https://www.degruyter.com/document/doi/10.1515/jisys-2013-0031/html?lang=de

%To Parameters:
%pic        the original (uint8 coloured-RGB) Picture where the inpainting Region can be found. 
%           The Inpainting Region has only full black RGB-Pixels. I am testing with small Picture 
%           of size 100x100 Pixels now, but I want to test with 1k and 2k Pictures later.
%mask       Logical 2d Array with same size as pic. 1 is part of Inpainting region
%patchsize  Opt. Param - Odd Int that says how big a patch is
%searchdist Opt. Param - Used in Step 3 for search areal for possible best patch
%alpha      The alpha used in the Paper for P-Calculating
%beta       The beta used in the Paper for P-Calculating
%gamma      The [small] omega used in the Paper for RC-Calculating


%[klein] omega wurde zu gamma umgeschrieben damit es weniger
%Verwechslungsgefahr mit [groß]Omega gibt.
%(the small omega was changed to gamma to make it easier to read and differ from big Omega in this code)
function [out] = patchBasedInpainting(pic,mask,patchsize,searchdist,alpha,beta,gamma)

    if ~exist("alpha","var")
        alpha = 0.2;
    end
    if ~exist("beta","var")
        beta = 0.8;
    end
    if ~exist("gamma","var")
        gamma = 0.7;
    end

    if gamma > 1 || 0 > gamma
        error("Gamma must be between zero and one!");
    end
    if (alpha+beta) > 1 || 0 > (alpha+beta)
        error("Sum of alpha and beta must be between zero and one!");
    end

    if ~exist("patchsize","var")
        patchsize = 9;
    else
        if isEven(patchsize) || (patchsize < 3)
            error("The Patchsize is not allowed to be odd!");
        end
    end

    ps = floor(patchsize/2);

    if ~exist("searchdist","var")
        searchdist = 33;
    end

    searchdist = searchdist + ps;
    %Searchdist beschreibt im welchem Umfeld um den höchst priorisierten
    %Pixel herum nach dem besten Patch gesucht werden soll. Dabei gibt
    %search dist bei der eingabe einen Wert vor, der dann noch anschließend
    %um ps erweitert wird.
    %(Searchdist tells in what area around a highst priority Pixel should
    %searched for the best patch. searchdist is give by caller function and
    %extended by ps)

    [y,x,N] = size(pic);

    %Phi fertiges Bildgebiet
    %(Phi is finished Picture Area)
    %Omega unfertiges Bildgebiet
    %(Omega is unfinished Picture Area - aka Inpainting Region)
    %deltaOmega Grenze zwischen Phi und Omega in Omega
    %(deltaOmega is the border between Phi and Omega and is part of Omega)

    Phi = pic .* uint8(~mask); %1 in Maske heißt Teil der Inpainting region (A One in the mask meas be part of inpaintign region)
    Omega = mask;
    Cmap = double(~Omega); % Map für das Vertrauen eiens Pixels (1 ist max) (Map for confidence for each Pixel - One is max)
    iterationcounter = 0;

    barmax = sum(Omega,"all");
    bartext = "Patchbased Inpainting in Iteration: ";
    bar = waitbar(0,append(bartext,string(0)));
    while 1 %Abbruch im Schritt 1 enthalten (Abort in Step 1)
        %Schritt 1: Bestimme deltaOmega (Step 1 - find deltaOmega)
        deltaOmega = zeros(sum(uint64(Omega),'all')*2,2);
        at = 1;

        if(isempty(deltaOmega)) %Falls die Restmaske leer ist, muss man auch nicht mehr zählen.
            break;              %(Break id there are no Pixel in Omega left)
        end

        for i = 1 + ps : y - ps
            for j = 1 + ps : x - ps
                if Omega(i,j) == 1
                    if sum(Omega(i-1:i+1,j-1:j+1),"all") < 9 %Überprüfe ob in der 9 Nachbarschaft des Pixels die Summe <9 ist
                        deltaOmega(at,:) = [i,j];       %(Check if in 9 Neigborhood sum is less to 9)
                        at = at + 1;
                    end
                end
            end
        end
        deltaOmega = deltaOmega(1:at-1,:);
        to = size(deltaOmega,1);

        %Schritt 2: Bestimmte die höchste Priorität
        %(Step 2: Find highest Priority Pixel)
        maxP = -inf;
        pixel_x = 0;
        pixel_y = 0;
        for at = 1 : to
            %Bestimme C bzw. RC
            %(Calculate C)
            p = deltaOmega(at,:);
            py = p(1,1);
            px = p(1,2);
            %[py,px] = deltaOmega(at,:);
            %patch = ~Omega(py-ps:py+ps,px-ps:px+ps);
            patch = Cmap(py-ps:py+ps,px-ps:px+ps);
            Cp = sum(patch,"all") / patchsize^2;

            RCp = (1-gamma) * Cp + gamma;

            %Bestimme D
            %(Calculate D)
            
            %(This is where I need to add Code. Somehow I need to extract a
            %Value between 0 and 1 by using gradients.
            Dp = 0;

            %Bestimme P
            %(Calculate P)
            Pp = alpha * RCp + beta * Dp;

            if Pp > maxP
                maxP = Pp;
                cForMaxP = Cp;
                pixel_y = py;
                pixel_x = px;
            end
        end
        if N == 3
            Phi_p = Phi(pixel_y-ps:pixel_y+ps,pixel_x-ps:pixel_x+ps,:);
        else
            Phi_p = Phi(pixel_y-ps:pixel_y+ps,pixel_x-ps:pixel_x+ps);
        end
        Omega_p = Omega(pixel_y-ps:pixel_y+ps,pixel_x-ps:pixel_x+ps);
        
        %Schritt 3: Bestimme den besten Patch aus Phi zum Problem
        %(Step 3: Find the best Patch inside Phi for the Pixel)
        best_value = inf;
        best_y = 0;
        best_x = 0;
        for i = pixel_y - searchdist : pixel_y + searchdist
            for j = pixel_x - searchdist : pixel_x + searchdist
                if (j+ps > x || i-ps < 1 ||j-ps < 1 || i+ps > y)
                    continue;%Die Maske darf din Bild Rand nicht überschreiten!
                            %(Mask is not allowed to exit Image Border!)
                end
                Q_mask = Omega(i-ps:i+ps,j-ps:j+ps); % Patch muss vollständig in Phi liegen!
                if mean2(Q_mask) == 0                % (Patch must be 100% inside Phi!)
                    if N == 3
                        Phi_q = Phi(i-ps:i+ps,j-ps:j+ps,:); 
                    else
                        Phi_q = Phi(i-ps:i+ps,i-ps:i+ps);
                    end

                    %err=SSE(Phi_p,Phi_q,Omega_p); %Der Einfluss von Omega ist viel zu groß, als das man hier die std. Funktion nutzen könnte.
                    %(Influence of Omega is to big to use the st. MATLAB Function)
                    err = immse(Phi_q,Phi_p);
                    if err < best_value
                        best_y = i;
                        best_x = j;
                        best_value = err;
                        bestPhiQ = Phi_q; % Nur für Debugzwecke
                                          %(Just for Debugging)
                    end
                end
            end
        end

        %Schritt 4 Fülle Omega mit den Pixeln Phi_q
        %(Step 4: Fill Omega with Pixels taken from Phi_q)
        for i = 0 : ps*2
            for j = 0 : ps*2
                if Omega(pixel_y-ps+i,pixel_x-ps+j) == 1 %Aktualisiere nur Pixel in Omega
                                                         %(Only update Pixel
                                                         %inside Omega)
                    Omega(pixel_y-ps+i,pixel_x-ps+j) = 0; %Nun Teil von Phi
                                                          %(Now a part of
                                                          %Phi)
                    if N == 3
                        Phi(pixel_y-ps+i,pixel_x-ps+j,:) = Phi(best_y-ps+i,best_x-ps+j,:);
                    else
                        Phi(pixel_y-ps+i,pixel_x-ps+j) = Phi(best_y-ps+i,best_x-ps+j);
                    end
                    Cmap(pixel_y-ps+i,pixel_x-ps+j) = cForMaxP; %Setzte Vertrauenswerte
                    %(update confidence)
                end
            end
        end


        %Schleifenende
        %Backpart of Looping
        iterationcounter = iterationcounter + 1;
        remainingPixels = sum(Omega,"all");
        %disp(string(remainingPixels));
        waitbar((barmax-remainingPixels)/barmax,bar,append(bartext,string(iterationcounter)));
    end
    close(bar);
    out = Phi;
end

%Optimierungs Ideen:
%Mittelwertsfilter Faltung auf der Maske um C zu bestimmen
%Um Omega zu finden sich den Rahmen merken - dies sollte realtiv viele
%Nullvergleiche vermeiden.

%(Idea for efficieny improfing:
%Mean Convolution for calculating C Values
%Save a Inner bounding area for Omega - so upper parts are not checked anymore)
%Finding ways to avoid zero checkings

Answer 1

[Shubh]

Hi Robin,

The data term (Dp) in patch-based inpainting algorithms plays a critical role, especially in how the inpainting process prioritizes different regions of the image. Generally, it's calculated based on image features like gradients or isophotes (lines of constant intensity) at the edges of the inpainting region. The goal is to give higher priority to those areas where the image structure (like edges or textures) is more prominent, facilitating a more coherent reconstruction of the missing parts.

Here's a way to calculate the data term (Dp) using gradient information:

1. Calculate Image Gradients: Compute the gradients of the image in both the x and y directions. This can be done using gradient operators like Sobel or Prewitt.
2. Determine Gradient at the Inpainting Boundary: For each pixel on the boundary of the inpainting region (part of deltaOmega in your code), calculate the gradient magnitude.
3. Normalize the Gradient Values: To get the data term in the range [0, 1], normalize the gradient magnitudes. One common approach is to divide by the maximum gradient magnitude found at the inpainting boundary.
4. Assign the Data Term: Use these normalized gradient magnitudes as your data term (Dp) for each boundary pixel.

Here's an example of how you might modify your MATLAB code to include this:

% ... [rest of your existing code] ...

% Step 2: Calculate Gradients of the Image
% You might need to convert the image to grayscale if it's RGB
grayImage = rgb2gray(pic);  
[Gx, Gy] = imgradientxy(grayImage, 'sobel');

% Step 2 Continued: Find highest Priority Pixel
for at = 1 : to
    % ... [existing code for calculating C and RC] ...

    % Calculate D
    % Determine the gradient magnitude at the boundary pixel
    py = p(1,1);
    px = p(1,2);
    gradMagnitude = sqrt(Gx(py, px)^2 + Gy(py, px)^2);

    % Normalize the gradient magnitude to get the data term
    % Assuming maxGradMagnitude is the maximum gradient magnitude found at the inpainting boundary
    Dp = gradMagnitude / maxGradMagnitude;  % Normalize to [0, 1]

    % Calculate P
    Pp = alpha * RCp + beta * Dp;
    % ... [rest of your existing priority calculation code] ...
end

% ... [rest of your existing code] ...

In this code snippet:

• 'imgradientxy' is used to calculate the x and y gradients of the image. You may need to pre-process your image into a suitable format (like grayscale) for this operation.
• 'sqrt(Gx(py, px)^2 + Gy(py, px)^2)' computes the gradient magnitude at each boundary pixel.
• 'gradMagnitude / maxGradMagnitude' normalizes this magnitude to a [0, 1] range.

Note: You'll need to determine 'maxGradMagnitude' beforehand by evaluating the gradient magnitudes at all boundary pixels and finding the maximum value. This pre-processing step ensures that your data term is normalized correctly.

This approach will give you a data term that reflects the strength of the image structure at the inpainting boundary, aiding in prioritizing areas with more significant structural information.

Hope this helps!

[Robin Breuer]

Well first of all thank you very much. I am now going to adopt this into my code to check this. It's helping anyways but there is still one thing I wonder about.

Is it really a good idea to perform normal filtering or normal gradient calculation on the image? I ask in particular because I'm still struggling with when, where or how I can use the pixels from the inpainting region. At least in my head, practically every pixel in the inpainting region is at an illegitimate value - so NaN, so to speak. If I carry out Sobel filtering, the theoretically non-existent pixels from the inpainting region are included in the gradient calculation in the edge area and I just wonder whether this is even allowed. I was already able to find out through testing that in the subsequent step where the squared error between the patch to be filled and the best patch is calculated, I should not include the inpainting region in the calculation.

In any case, I will now implement the gradient solution first. I would also like to try the isophote solution - do you have a basic idea for that? What I had seen before was to calculate the isophotes (e.g. using this methodology: https://matlab.mathworks.com/open/fileexchange/v1?id=48811) and then use a mean value filter from the isophote values for angle and magnitude to calculate a value, which poses three questions to me:

1: Are the values from the inpainting region allowed to be included in the mean value filtering?

2: Can I actually combine magnitude and angle into a meaningful value?

3: In my main source paper the image is rotated 90° for this calculation - why?

Either way, thank you very much so far.

Edit:

I tried out to use the normal. magnitudes and indeed I got better results as before but not as good as I hoped for.