Correlation based similarity measures-Summary

11 04 2010

Correlation based matching typically produces dense depth maps by calculating the disparity at each pixel within a neighborhood. This is achieved by taking a square window of certain size around the pixel of interest in the reference image and finding the homologous pixel within the window in the target image, while moving along the corresponding scanline. The goal is to find the corresponding (correlated) pixel within a certain disparity range d (d E [0,….,dmax]) that minimizes the associated error and maximizes the similarity.

In brief, the matching process involves computation of the similarity measure for each disparity value, followed by an aggregation and optimization step. Since these steps consume a lot of processing power, there are significant speed-performance advantages to be had in optimizing the matching algorithm.

The images can be matched by taking either left image as the reference (left-to-right matching, also known as direct matching) or right image as the reference (right-to-left matching, also known as reverse matching) [2].

Similarity Measure Formula
Sum of Absolute Differences (SAD) Sum of Absolute Differences
Zero-mean Sum of Absolute Differences (ZSAD)
Locally scaled Sum of Absolute Differences (LSAD)
Sum of Squared Differences (SSD) Sum of Squared Differences
Zero-mean Sum of Squared Differences (ZSSD)
Locally scaled Sum of Squared Differences (LSSD)
Normalized Cross Correlation (NCC) Normalized Cross Correlation
Zero-mean Normalized Cross Correlation (ZNCC)
Sum of Hamming Distances (SHD) Sum of Hamming Distances

Sum of Absolute Differences (SAD) is one of the simplest of the similarity measures which is calculated by subtracting pixels within a square neighborhood between the reference image I1 and the target image I2 followed by the aggregation of absolute differences within the square window, and optimization with the winner-take-all (WTA) strategy [1]. If the left and right images exactly match, the resultant will be zero.

In Sum of Squared Differences (SSD), the differences are squared and aggregated within a square window and later optimized by WTA strategy. This measure has a higher computational complexity compared to SAD algorithm as it involves numerous multiplication operations.

Normalized Cross Correlation is even more complex to both SAD and SSD algorithms as it involves numerous multiplication, division and square root operations.

Sum of Hamming Distances is normally employed for matching census-transformed images (can be used on images that have not been census transformed) by computing bitwise-XOR of the values in left and right images, within a square window. This step is usually followed by a bit-counting operation which results in the final Hamming distance score.

Example: Tsukuba

TsukubaLeftLeft Image TsukubaRightRight Image
SAD Disparity Map ZSAD Disparity Map LSAD Disparity Map
SSD Disparity Map ZSSD Disparity Map LSSD Disparity Map
NCC Disparity Map ZNCC Disparity Map TsukubaSHD9x9DispRange=0-16ColorSHD Disparity Map
TsukubaGroundTruthColor

Ground Truth Disparity Map

Disparity Range: 0-16 ,Window Size: 9×9, The hotter the color, the closer it is; the cooler the color the farther it is.

Generic MATLAB code:

% *************************************************************************
% Title: Function-Compute Correlation between two images using various 
% similarity measures with Left Image as reference.
% Author: Siddhant Ahuja
% Created: March 2010
% Copyright Siddhant Ahuja, 2010
% Inputs: 
% 1. Left Image (var: rightImage), 
% 2. Right Image (var: leftImage),
% 3. Correlation Window Size (var: corrWindowSize), 
% 4. Minimum Disparity in X-direction (var: dMin), 
% 5. Maximum Disparity in X-direction (var: dMax),
% 6. Method used for calculating the correlation scores (var: method)
% Valid values include: 'SAD', 'LSAD', 'ZSAD', 'SSD', 'LSSD', ZSSD', 'NCC',
% 'ZNCC'
% Outputs: 
% 1. Disparity Map (var: dispMap), 
% 2. Time taken (var: timeTaken)
% Example Usage of Function: [dispMap, timeTaken]=denseMatch('tsukuba_left.tiff', 'tsukuba_right.tiff', 9, 0, 16, 'ZNCC');
% *************************************************************************
function [dispMap, timeTaken]=denseMatch(rightImage, leftImage, corrWindowSize, dMin, dMax, method)
% Grab the image information (metadata) of left image using the function imfinfo
leftImageInfo=imfinfo(leftImage);
% Grab the image information (metadata) of right image using the function imfinfo
rightImageInfo=imfinfo(rightImage);
% Since Dense Matching is applied on a grayscale image, determine if the
% input left image is already in grayscale or color
if(getfield(leftImageInfo,'ColorType')=='truecolor')
% Read an image using imread function, convert from RGB color space to
% grayscale using rgb2gray function and assign it to variable leftImage
    leftImage=rgb2gray(imread(leftImage));
else if(getfield(leftImageInfo,'ColorType')=='grayscale')
% If the image is already in grayscale, then just read it.        
        leftImage=imread(leftImage);
    else
        error('The Color Type of Left Image is not acceptable. Acceptable color types are truecolor or grayscale.');
    end
end
% Since Dense Matching is applied on a grayscale image, determine if the
% input right image is already in grayscale or color
if(getfield(rightImageInfo,'ColorType')=='truecolor')
% Read an image using imread function, convert from RGB color space to
% grayscale using rgb2gray function and assign it to variable rightImage
    rightImage=rgb2gray(imread(rightImage));
else if(getfield(rightImageInfo,'ColorType')=='grayscale')
% If the image is already in grayscale, then just read it.        
        rightImage=imread(rightImage);
    else
        error('The Color Type of Right Image is not acceptable. Acceptable color types are truecolor or grayscale.');
    end
end
% Find the size (columns and rows) of the left image and assign the rows to
% variable nrLeft, and columns to variable ncLeft
[nrLeft,ncLeft] = size(leftImage);
% Find the size (columns and rows) of the right image and assign the rows to
% variable nrRight, and columns to variable ncRight
[nrRight,ncRight] = size(rightImage);
% Check to see if both the left and right images have same number of rows
% and columns
if(nrLeft==nrRight && ncLeft==ncRight)
else
    error('Both left and right images should have the same number of rows and columns');
end
% Convert the left and right images from uint8 to double
leftImage=im2double(leftImage);
rightImage=im2double(rightImage);
% Check the size of window to see if it is an odd number.
if (mod(corrWindowSize,2)==0)
    error('The window size must be an odd number.');
end
% Check whether minimum disparity is less than the maximum disparity.
if (dMin>dMax)
    error('Minimum Disparity must be less than the Maximum disparity.');
end
% Create an image of size nrLeft and ncLeft, fill it with zeros and assign
% it to variable dispMap
dispMap=zeros(nrLeft, ncLeft);
% Find out how many rows and columns are to the left/right/up/down of the
% central pixel based on the window size
win=(corrWindowSize-1)/2;
% The objective of CC, NCC and ZNCC is to maxmize the
% correlation score, whereas other methods try to minimize
% it.
maximize = 0;
if strcmp(method,'NCC') || strcmp(method,'ZNCC')
    maximize = 1;
end
tic; % Initialize the timer to calculate the time consumed.
for(i=1+win:1:nrLeft-win)
    % For every row in Left Image
    for(j=1+win:1:ncLeft-win-dMax)
        % For every column in Left Image
        % Initialize the temporary variable to hold the previous
        % correlation score
        if(maximize)
            prevcorrScore = 0.0;
        else
            prevcorrScore = 65532;
        end
        % Initialize the temporary variable to store the best matched
        % disparity score
        bestMatchSoFar = dMin;
        for(d=dMin:dMax)
            % For every disparity value in x-direction
            % Construct a region with window around central/selected pixel in left image
            regionLeft=leftImage(i-win : i+win, j-win : j+win);
            % Construct a region with window around central/selected pixel in right image
            regionRight=rightImage(i-win : i+win, j+d-win : j+d+win);
            % Calculate the local mean in left region
            meanLeft = mean2(regionLeft);
            % Calculate the local mean in right region
            meanRight = mean2(regionRight);
            % Initialize the variable to store temporarily the correlation
            % scores
            tempCorrScore = zeros(size(regionLeft));
            % Calculate the correlation score
            if strcmp(method,'SAD')
                tempCorrScore = abs(regionLeft - regionRight);
            elseif strcmp(method,'ZSAD')
                tempCorrScore = abs(regionLeft - meanLeft - regionRight + meanRight);
            elseif strcmp(method,'LSAD')
                tempCorrScore = abs(regionLeft - meanLeft/meanRight*regionRight);
            elseif strcmp(method,'SSD')
                tempCorrScore = (regionLeft - regionRight).^2;
            elseif strcmp(method,'ZSSD')
                tempCorrScore = (regionLeft - meanLeft - regionRight + meanRight).^2;          
            elseif strcmp(method,'LSSD')
                tempCorrScore = (regionLeft - meanLeft/meanRight*regionRight).^2;
            elseif strcmp(method,'NCC')
                % Calculate the term in the denominator (var: den)
                den = sqrt(sum(sum(regionLeft.^2))*sum(sum(regionRight.^2)));
                tempCorrScore = regionLeft.*regionRight/den;
            elseif strcmp(method,'ZNCC')
                % Calculate the term in the denominator (var: den)
                den = sqrt(sum(sum((regionLeft - meanLeft).^2))*sum(sum((regionRight - meanRight).^2)));
                tempCorrScore = (regionLeft - meanLeft).*(regionRight - meanRight)/den;
            end
            % Compute the final score by summing the values in tempCorrScore,
            % and store it in a temporary variable signifying the distance
            % (var: corrScore)
            corrScore=sum(sum(tempCorrScore));
            if(maximize)
                if(corrScore>prevcorrScore)
                    % If the current disparity value is greater than
                    % previous one, then swap them
                    prevcorrScore=corrScore;
                    bestMatchSoFar=d;
                end
            else
                if (prevcorrScore > corrScore)
                    % If the current disparity value is less than
                    % previous one, then swap them
                    prevcorrScore = corrScore;
                    bestMatchSoFar = d;
                end
            end
        end
        % Store the final matched value in variable dispMap
        dispMap(i,j) = bestMatchSoFar;
    end
end
% Stop the timer to calculate the time consumed.
timeTaken=toc;
end

For other MATLAB Codes, please visit:

1. Sum of Absolute Differences
2. Sum of Squared Differences
3. Normalized Cross Correlation
4. Sum of Hamming Distances

References

1. T. Kanade, H. Kano, and S. Kimura, “Development of a video-rate stereo machine,” in Image UnderstandingWorkshop, Monterey,CA, 1994, p. 549–557.
2. D. Scharstein and R. Szeliski, “A taxonomy and evaluation of dense two-frame stereo correspondence algorithms,” International Journal of Computer Vision, vol. 47(1/2/3), pp. 7-42, Apr. 2002.
3. S. Chambon and A. Crouzil. évaluation et comparaison de mesures de corrélationrobustes aux occultations. Rapport de recherche 2002-34-R, IRIT, Université PaulSabatier, Toulouse, France, December 2002.





Correlation based similarity measure-Normalized Cross Correlation (NCC)

20 05 2009

Example 1: Tsukuba

TsukubaLeftLeft Image TsukubaRightRight Image
TsukubaNCC9x9DispRange=0-16ColorNCC Disparity Map
Disparity Range: 0-16
Window Size: 9×9
The hotter the color, the closer it is; the cooler the color the farther it is.

Example 2: Stereogram

StereogramLeftLeft Image StereogramRightRight Image
StereogramNCC9x9DispRange=0-16ColorNCC Disparity Map
Disparity Range: 0-16
Window Size: 9×9
The hotter the color, the closer it is; the cooler the color the farther it is.

MATLAB Code-NCC Right to Left Matching

% *************************************************************************
% Title: Function-Compute Correlation between two images using the 
% similarity measure of Normalized Cross Correlation (NCC) with Right Image 
% as reference.
% Author: Siddhant Ahuja
% Created: May 2008
% Copyright Siddhant Ahuja, 2008
% Inputs: Left Image (var: leftImage), Right Image (var: rightImage),
% Window Size (var: windowSize), Minimum Disparity (dispMin), Maximum
% Disparity (dispMax)
% Outputs: Disparity Map (var: dispMap), Time taken (var: timeTaken)
% Example Usage of Function: [dispMap, timeTaken]=funcNCCR2L('StereogramLeft.jpg', 'StereogramRight.jpg', 9, 0, 16);
% *************************************************************************
function [dispMap, timeTaken]=funcNCCR2L(leftImage, rightImage, windowSize, dispMin, dispMax)
try 
    % Grab the image information (metadata) of left image using the function imfinfo
    leftImageInfo=imfinfo(leftImage);
    % Since NCCR2L is applied on a grayscale image, determine if the
    % input left image is already in grayscale or color
    if(getfield(leftImageInfo,'ColorType')=='truecolor')
    % Read an image using imread function, convert from RGB color space to
    % grayscale using rgb2gray function and assign it to variable leftImage
        leftImage=rgb2gray(imread(leftImage));
        % Convert the image from uint8 to double
        leftImage=double(leftImage);
    else if(getfield(leftImageInfo,'ColorType')=='grayscale')
    % If the image is already in grayscale, then just read it.        
            leftImage=imread(leftImage);
            % Convert the image from uint8 to double
            leftImage=double(leftImage);
        else
            error('The Color Type of Left Image is not acceptable. Acceptable color types are truecolor or grayscale.');
        end
    end
catch
    % if it is not an image but a variable
    leftImage=leftImage;
end
try
    % Grab the image information (metadata) of right image using the function imfinfo
    rightImageInfo=imfinfo(rightImage);
    % Since NCCR2L is applied on a grayscale image, determine if the
    % input right image is already in grayscale or color
    if(getfield(rightImageInfo,'ColorType')=='truecolor')
    % Read an image using imread function, convert from RGB color space to
    % grayscale using rgb2gray function and assign it to variable rightImage
        rightImage=rgb2gray(imread(rightImage));
        % Convert the image from uint8 to double
        rightImage=double(rightImage);
    else if(getfield(rightImageInfo,'ColorType')=='grayscale')
    % If the image is already in grayscale, then just read it.        
            rightImage=imread(rightImage);
            % Convert the image from uint8 to double
            rightImage=double(rightImage);
        else
            error('The Color Type of Right Image is not acceptable. Acceptable color types are truecolor or grayscale.');
        end
    end
catch
    % if it is not an image but a variable
    rightImage=rightImage;
end
% Find the size (columns and rows) of the left image and assign the rows to
% variable nrLeft, and columns to variable ncLeft
[nrLeft,ncLeft] = size(leftImage);
% Find the size (columns and rows) of the right image and assign the rows to
% variable nrRight, and columns to variable ncRight
[nrRight,ncRight] = size(rightImage);
% Check to see if both the left and right images have same number of rows
% and columns
if(nrLeft==nrRight && ncLeft==ncRight)
else
    error('Both left and right images should have the same number of rows and columns');
end
% Check the size of window to see if it is an odd number.
if (mod(windowSize,2)==0)
    error('The window size must be an odd number.');
end
% Check whether minimum disparity is less than the maximum disparity.
if (dispMin>dispMax)
    error('Minimum Disparity must be less than the Maximum disparity.');
end
% Create an image of size nrLeft and ncLeft, fill it with zeros and assign
% it to variable dispMap
dispMap=zeros(nrLeft, ncLeft);
% Find out how many rows and columns are to the left/right/up/down of the
% central pixel based on the window size
win=(windowSize-1)/2;
tic; % Initialize the timer to calculate the time consumed.
for(i=1+win:1:nrLeft-win)
    for(j=1+win:1:ncLeft-win-dispMax)
        prevNCC = 0.0;
        bestMatchSoFar = dispMin;
        for(dispRange=dispMin:1:dispMax)
            ncc=0.0;
            nccNumerator=0.0;
            nccDenominator=0.0;
            nccDenominatorRightWindow=0.0;
            nccDenominatorLeftWindow=0.0;
            for(a=-win:1:win)
                for(b=-win:1:win)
                   nccNumerator=nccNumerator+(rightImage(i+a,j+b)*leftImage(i+a,j+b+dispRange));
                   nccDenominatorRightWindow=nccDenominatorRightWindow+(rightImage(i+a,j+b)*rightImage(i+a,j+b));
                   nccDenominatorLeftWindow=nccDenominatorLeftWindow+(leftImage(i+a,j+b+dispRange)*leftImage(i+a,j+b+dispRange));
                end
            end
            nccDenominator=sqrt(nccDenominatorRightWindow*nccDenominatorLeftWindow);
            ncc=nccNumerator/nccDenominator;
            if (prevNCC < ncc)
                prevNCC = ncc;
                bestMatchSoFar = dispRange;
            end
        end
        dispMap(i,j) = bestMatchSoFar;
    end
end
% Stop the timer to calculate the time consumed.
timeTaken=toc;








Follow

Get every new post delivered to your Inbox.

Join 80 other followers