Skip to content

B. Linear Filters

Jump to bottom
lcdurann edited this page Jul 13, 2020 · 9 revisions

Two types of linear filters were considered on the project, namely the MACE filters and the HBCOM filters.

MACE Filters

Linear filters are computed by forming linear combinations of training images in frequency space. Although this sounds simple, the characteristics of the correlation output can be controlled by choosing the superposition coefficients so that a particular property is optimized. That is the particular case of linear filters like the MACE. The MACE Filter minimizes the Average Correlation Energy of the correlation output (hence the name). But also, this particular filter is designed so that the correlation output has a fixed peak size for the training images used. So, the MACE filter is a constrained linear filter.

Mathematical definition

When training images are viewed as column vectors xi, and the MACE correlation filter, as a column vector h, the correlation output of the i-th training image is simply gi = Xi · h. Where the dot product is defined in a complex space, and Xi is a diagonal matrix with xi along its diagonal. The Average Correlation Energy is simply

N Eave = Σi Ei = Σi gi · gi = Σi h · Di · h

Where Di is a diagonal matrix that contains the power spectrum of the i-th image. If the matrix D is defined as one that contains the average power spectrum of the training set, then:

Eave = Σi h · D · h

The MACE filter is chosen to minimize this quantity constrained by the requirement that the peak intensities for the training set, X · h = u. Where X is a matrix whose columns are the images inside the training set, and u is a predefined vector, usually with all components equal to one. By means of the Lagrange multipliers method it can be shown that the optimization process leads to:

h = D-1<\b> · X<\b> · (X<\b> · D-1<\b> · X<\b>)-1 u<\b>

Computation algorithm

The algorithm for the computation of a MACE filter without registration is very straightforward:

  1. Compute the Average Density Spectrum and the X matrix: This step can be implemented thusly in MATLAB.
%% Generating matrixes for filter computation
for i = 3:numel(images)
    % Read image from directory
    filename = [basedir '/' dirname '/' images(i).name];
    im = fftshift(fft2(rgb2gray(imread(filename)),m,n));
    % Update matrix of images
    X(:,i-2) = im(:);
    % Update Power Spectra Matrix
    P = P + im(:) .* conj(im(:)) * 1.0/(numel(images)-2);
end
  1. Compute the product D-1 · X: Since images have resolutions of 1280x720, the average power spectrum matrix is of large size, and can easily cause overflow errors. Hence, in MATLAB, the product can be computed as follows:
%% Computation of filter
dims = size(X);
iP = 1 ./ P;
% Compute DX & iDX matrix
for i = 1:dims(2)
    iDX(:,i) = iP .* X(:,i);
end
  1. Perform the matrix product computation for the filter: This corresponds to computing the triple matrix product in the mathematical definition of the MACE filter. In MATLAB:
XDX = (conj(X') * iDX);
filter = (iDX / XDX) * u;
filter = reshape(filter,orsize(1),orsize(2));

For a more refined construction of the filter some sort of registration of the training images is necessary. This is an iterative process in which the images that enter in the filter are chosen with the criteria that a normalization metric of the correlation output, when applied to a candidate training image should be the worst over the training set. In other words, the filter is updated to contain the worst performing image in the training set. Usually, some sort of alignment of the output peaks is necessary. The registration process is somewhat lengthy and the reader is referred to the function MACE_Filter.m for more details.

Expected Performance

HCOM Filters

The composite filter, represented by HCOM is also based on a linear combination involving a subset from the subject's image database, such subset is known as the Training set or the Reference set. This type of linear filter is well consider for the intended application since it allows to merge several versions of reference images and it also can be useful to solve sensitivity issues thanks to rotations or scaling operations [1].

Mathematical definition

The way this filter can be constructed is by a simple standard linear combination:

HCOM= Σi ai Ri

Where the weights ai are specifically chosen to optimize some cost function which depends on the desire application, whereas the term Ri represents the Fourier spectrum of the i-th training image used to construct the filter [1].

As was mentioned before, the proper election of the linear combination coefficients is a crucial step to enhance the filter behaviour, ensuring its correct performance. Thus, one of these selection criteria is based on the Peak-to-Correlation Energy value or PCE, which is defined as the energy of the peak on the plane, normalized to the total energy of the plane.

Now, when constructing the HCOM filter by means of the PCE criteria is important to remark that apart from the reference set used on the linear combination, one must select another image from the database called Filter image whose purpose consists on serve as the sample from which the linear combination coefficients are to be found. Hence, based on the study done by de la Tocnaye, Quémener and Pétillot [2], some minor changes were carried out so the process to select the weights is done as follows:

  1. Compute individually all the Fourier transforms of the training set images, that is, calculate Ri for all values of i and also calculate the spectrum T for a suitable filter image, which for this specific process is going to be treated as a test image.
  2. Perform the product conj(T)*Ri for all values of i and calculate the inverse Fourier transform in each case. This would be defined as a Partial correlation or PC, since it represents the correlation between the target and a single reference. Therefore, if there are N samples on the training set, one gets virtually N different Partial Correlation Planes.
  3. Now, define a total correlation plane as the inverse transform of the convolution between the target image and a spectrum compose of all the training samples. Thus, this plane is represented as:

C= FT-1{conj(T)*Σi Ri}

Note that C represents the total correlation plane for a case in which the convolution is computed between the test sample and a linear combination of the training images but with unit valued weights.
  1. Now, since each partial correlation gives a different peak, a PCE value for each of them can be defined as the energy of the respective peak, normalized to the energy of the plane C defined above:

PCEi=E_peaki/EC

Where E_peaki represents the peak energy on the i-th partial correlation and EC is the energy of the total correlation plane described before. The respective energies are calculated as:

E_peaki=|max(PCi)|2 ; EC=Σnm|c(n,m)|2

Where max(PCi) represents the peak of the i-th partial correlation and c(n,m) is the value assign to the point n,m on the plane C and correspondingly |c(n,m)|2 is the energy of the point in question.

  1. Finally the weights are found as ai=[PCEi]-b . where b is is a constant value between 0.5 and 1.5. For this project the chosen value was b=0.8.

Therefore, the HCOM filter is constructed by means of a sample of training images and a filter image which is treated as a target in order to find the optimal values for the linear combination coefficients. Thus, the filter is finally design and store for the respective subject under which it was made, this means that each individual is assigned their own filter, built following the previous steps.

HBCOM filters

A binarized composite filter can also be designed and is defined as follows [2]:

HBCOM= sgn[Re{Σi ai Ri}]

Where Re{x} represents the real part of x and sgn is the sign function which is defined as

sgn(x)= 1 if x>0; 0 if x=0; -1 if x<0.

The method used to calculate the weights remains the same and the process to design the filter by computational methods is described in the next section

Computation algorithm

Before designing the HBCOM filter, is worth mentioning that in order to get its best performance, the images must go through a pre-processing stage, where the entire set of samples (which includes the training set and the filter image) must follow certain specified steps for their acquisition and refinement. The filter is constructed in MATLAB, by means of the following algorithm:

  1. Choose the filter image from database:
    %% Read Reference Image
    % The system path of the image is generated using directory object
    
    refimagpath = [base '/*sample' num2str(refimag) '.png'];
    RefImagLoc = dir(refimagpath);
    fullname = [curr_loc dataFolder '/' RefImagLoc.name];
    sample = imread(fullname);
  1. Calculate the test image FT and initialize other variables
%% Define variables for filter computation
    fftsamp = fft2(sample);             % Template for coeff optimization
    filter = zeros(size(fftsamp));       % Filter
    totalxcorr = zeros(size(fftsamp));
  1. Select size of training set and perform process to find weights
for j = 1:subspace_size
            % read image from training set
            filename = [base '/' Data(j).name];
            im = imread(filename);
            usedImages{end+1} = filename;
            disp(['Used image: ' Data(j).name]);
            % Compute xcorr with template
            fftim = fft2(im);
            corrplane = abs(fftshift(ifft2(...
                conj(fftsamp) .* fftim...
                )));
            % cumulative xcorr
            totalxcorr = totalxcorr + corrplane;
            % Plane energy
            PE = max(corrplane(:))^2;
            % Update filter
            filter = filter + (PE)^(-0.8) * fftim;
        end
        % Total energy
        etot = sum(abs(totalxcorr).^2,'all');
        % Binarize + real
        filter = sign(real(1.0/etot^(-0.8) * filter));

As was mentioned before, the filter is calculated and store for future implementations using a specific filter image and a chosen number of training samples. However, it is important to pick out properly the image in which the filter bases its construction on, for more details about the HBCOM filter construction, the reader is refereed to the function HBCOM_Filter.m.

Expected Performance

Once the filter is completed and ready to use, it is expected to be capable of accurately verifying a test sample belonging to the subject's database but ofuscarse, different from the training set and the filter image. Moreover, the quality of the filter performance is also determined by its discrimination capacity, that is, whether or not the filter can distinguish a false input sample. Hence, a good behaviour of the HBCOM filter should give results as the one exhibit bellow.

Based on the image shown above, is clear that the HBCOM filters are able to properly discriminate between a false sample and a real one, additionally, the correlation plane for the true class gives a truly sharp and bright peak with almost no noise, which is a sign of reliability on its performance. Additionally, is worth mentioning that these filter enhances its result on the correlation plane in terms of peak intensity and brightness, when a binirazation is done during the preprocessing stage, however this discussion is expanded on the last two sections Methodology and Simulations

The fundamental references of the project are:

Clone this wiki locally