Working with Bumblebee2 stereo camera and FireflyMV monocular camera in Ubuntu 11.10 + ROS Electric

I had a hard time trying to get the stereo vision camera  – Bumblebee2 and monocular camera – FireFlyMV from Point Grey, working with Ubuntu 11.10 and ROS Electric. Finally ironed out all the bugs, and wrote this tutorial which takes you step by step to getting it running on my laptop.

I have tested it with the following environment setup:

• Origin PC Eon 15s with Firewire 1394b port
• Point Grey’s Bumblebee2 BB2-03S2C (640×480 at 48 fps and Raw16 pixel format)
• Point Grey’s FireflyMV FMVU-13S2C (1328×1048 at 23FPS, 640×480 at 60 FPS)
• Ubuntu 11.10 (Oneiric)
• ROS Electric

1. If you already have a Firewire 1394a/b port in your laptop, then you can ignore this step and move on. If you are using an express card, then you might have to go through the installation and check if firewire functionality works. Make sure that the kernel of the Ubuntu 11.10 OS recognises the express card. I tested this part on a Toshiba Tecra R700 laptop with Bytec Express card BT-ECU2FW which has 2 Firewire and 1 USB port. The express card and firewire functionality did not work out of the box. I followed the instructions from : https://help.ubuntu.com/community/ExpressCard to install my card.

    

   Setup the kernel option for PCI Express Hotplug driver, by opening a Terminal window and typing:

   $ sudo modprobe pciehp pciehp_force=1

   Then, load the ACPI PCI hotplug driver:

   $ sudo modprobe acpiphp

   Finally, setup the drivers to be loaded automatically on startup:

   $ gksudo gedit /etc/modules

   This will open up a file. Add the following line to the end of it:

   acpiphp

   Save and exit.
   Restart the computer.
   Once in the OS, remove the express card and plug it back in again. Open a terminal window, and type:

   $ dmesg

   You should see something similar to the following, as reported by “firwire_ohci” and “firewire_core”:

   firewire_ohci 0000:03:00.0: enabling device (0000 -> 0003)
   firewire_ohci 0000:03:00.0: PCI INT A -> GSI 18 (level, low) -> IRQ 18
   firewire_ohci 0000:03:00.0: setting latency timer to 64
   firewire_ohci: Added fw-ohci device 0000:03:00.0, OHCI v1.10, 4 IR + 8 IT contexts, quirks 0x11
   firewire_core: created device fw0: GUID 001106660000015c, S400
   

   If the firewire_core reports that a device was created, then everything is good to go. Else, your express card was not recognized, and so you may need to troubleshoot this by searching the forums online.

2. Check for 32/64-bit processor: Next, you need to figure out whether your computer supports 32-bit or 64-bit applications. This is needed for the next step where we will download a specific archive file depending on the type of support. To find out what your system capability is, in a terminal window type the following command:

   $ grep --color=always -iw lm /proc/cpuinfo

   If you see lm highlighted in red, then your processor is capable of 64-bit computing.

3. Download and install Libdc1394 library:
   Install all the packages that will be needed for the following step:

   sudo apt-get install libraw1394-11 libgtk2.0-0 libgtkmm-2.4-dev libglademm-2.4-dev libgtkglextmm-x11-1.2-dev libusb-1.0-0 libusb-dev libusb-1.0-0-dev libraw1394-dev git autoconf libtool pkg-config libsdl1.2-dev

   Now we will grab the source code from the online repository:

   cd

   git clone git://libdc1394.git.sourceforge.net/gitroot/libdc1394/libdc1394/

   cd ~/libdc1394/libdc1394

   autoreconf -i -s

   ./configure

   You should see the following output at the end:

   Configuration (libdc1394):
   
       Linux support:                      Enabled
         Legacy drivers (libraw1394):      Enabled
       Mac OS X support:                   Disabled (Mac OS not detected)
       Windows support:                    Disabled (Windows not detected)
       IIDC-over-USB support:              Enabled
       Build examples:                     Yes
         Build SDL/OpenGL examples:        Yes
         Build XV examples:                No
         Build V4L examples:               No
   

   If you do not see the IIDC-over-USB support as enabled, then you would have to re-install libusb packages. FireflyMV cameras implement IIDC 1394-based Digital Camera Specification, so they would not appear in the guvcview utility. However, if the libdc1394 is properly set up with IIDC over USB enabled, they would show up in a program called coriander. We will install coriander later.

   sudo make

   sudo make install

    

   If for some odd reason the above does not work for you, then you might have to try the harder and NOT PREFERRED way as outlined below. For most users, you can skip the following instructions and move on to the next step. YOu can grab the libdc1394v2 package from online which is needed for working with the sensor. In a new terminal window, type the following in your home directory:

   $ svn co http://cornell-ros-pkg.googlecode.com/svn/trunk/bumblebee2/camera_drivers_experimental/libdc1394v2/

   Once the download finishes, type the following:

   $ cd libdc1394v2/

   $ make

   You will get the following error, and compilation will terminate:

   dc1394_vloopback.c:49:28: fatal error: linux/videodev.h: No such file or directory

   To solve this problem, we will open up the file dc1394_vloopback.c and change the naming of the include file. Type the following in the terminal:

   $ gedit ~/libdc1394v2/build/libdc1394v2-svn/examples/dc1394_vloopback.c

   Go to line number 49, and replace

   linux/videodev.h 

   with

   libv4l1-videodev.h 

   Save and exit.

   $ cd ~/libdc1394v2/build/libdc1394v2-svn

   $ ./configure

   $ sudo make

   Everything should compile just fine now.

   $ sudo make install

4. Installing Coriander to check camera output:

   cd

   git clone git://coriander.git.sourceforge.net/gitroot/coriander/coriander/ 

   sudo apt-get install ftplib3 ftplib-dev libxv-dev x11proto-video-dev libgnomeui-dev

   ./configure

   sudo make

   sudo make install

   sudo coriander

   This will popup the coriander program, and if you have Bumblebee2 and/or FireflyMV connected then you can view the output. You should see the following window:
   
   Select the Services Tab from above, and make sure Format is Format 7, Mode_3 for Bumblebee2. Under Options, within the same tab, select Nearest as the coding scheme, BGGR Bayer pattern is selected by default, set the stereo coding scheme to St. Interlaced as shown below:
   
   Click on the Display button and you should see another window pop up with the stereo images, one on top of another

5. Installation of FlyCapture program from PointGrey: Create your own user account at PointGrey’s website, login, and from this link: http://ptgrey.com/support/downloads/download.asp, grab the appropriate 64 bit FlyCapture2 v2.3 Beta 8, archive file from their website, found under the heading of Software and Drivers for Imaging Products (Linux operating systems).
   Extract the tar file into the folder where you downloaded it.
   Open a terminal and change directory to the extracted archive folder from within the terminal:

   $ cd ~/Downloads/flycapture2-2.3.1.8-amd64/

   Then, type the following to install the program:

   $ sudo sh install_flycapture.sh

   The installation should finish successfully. Restart the computer.
   After reboot, you should now see a menu item for Point Grey Research’s FlyCap 2 program.
   Power up your BUmblebee2 sensor and connect the firewire cable.
   Launch the Flycap2 application, and you should now see the following Camera Selection window:
   
   Click on the Configure Selected button, and you should now see the following window appear:
   
   From the left hand menu, select Custom Video Modes and on the right hand pane, make sure Mode: 3 is selected. Then click on Apply on the bottom left. You can now close this window.

6. Install the Bumblebee2 ROS Package: Now we will download the Bumblebee2 ROS package into the home folder, in order to work with the sensor and publish the left and right images as topics. In a new terminal window, type:

   $ svn co http://cornell-ros-pkg.googlecode.com/svn/trunk/bumblebee2/

   Once the download has finished, type the following in the terminal:

   $ rm -rvf ~/bumblebee2/camera_drivers_experimental/libdc1394v2/

   We are basically removing the libdc1394v2 folder as we do not need it in that folder anymore (we have already compiled and installed it separately before). Now, we will remove the sr4k folder, which is also not required:

   $ rm -rvf ~/bumblebee2/camera_drivers_experimental/sr4k/

   Now we will add the locations of libdc1394v2 and bumblebee2 to the ROS_PACKAGE_PATH, so that ROS can find them for the next step of compilation:

   $ cd

   $ sudo gedit .bashrc

   Enter your password when prompted, and add the following two lines to the end of the file:

   export ROS_PACKAGE_PATH=~/libdc1394v2:$ROS_PACKAGE_PATH
   export ROS_PACKAGE_PATH=~/bumblebee2:$ROS_PACKAGE_PATH

   Save and exit.
   In a new terminal window, type:

   $ cd ~/bumblebee2/camera_drivers_experimental/

   $ rosmake --rosdep-install

   The build will succeed, however when you try to launch the stereo camera you will get the following error:

   libdc1394 error: Format_7 Error_flag_1 is set: in dc1394_format7_set_roi (format7.c, line 864): Handshaking failed after setting size and position
   
   [ WARN] [1334707170.281146370]: Failed to set format7_3 parameters for Bumblebee2
   [FATAL] [1334707170.288842875]: [camera] exception opening device: [Camera1394::open]: Failed to open device!
   Segmentation fault, stopping camera driver.
   

   To solve it, do the following:

   $ gedit ~/bumblebee2/camera_drivers_experimental/camera1394/src/dev_camera1394.cpp

   Do a search for DC1394_USE_MAX_AVAIL, 0, 0, and following that, replace 1024 with 640, and 768 with 480. Now, search for unsigned int cols = 1024 and replace 1024 with 640, and in the same line 768 with 480. Save and exit.
   Type the following in the terminal window:

   $ cd ~/bumblebee2/camera_drivers_experimental/

   $ rosmake --rosdep-install

   The build will succeed again.
   Next we will edit the calibration files.

   $ gedit ~/bumblebee2/camera_calibration/bumblebee_left.yaml

   Edit the image_width: to be 640 and image_height to be 480. Save and exit.

   $ gedit ~/bumblebee2/camera_calibration/bumblebee_right.yaml

   Edit the image_width: to be 640 and image_height to be 480. Save and exit.

   $ gedit ~/bumblebee2/camera_drivers_experimental/camera1394/launch/bumblebee2.yaml 

   Edit fps to be 48, bayer_pattern to be BGGR, and ensure bayer_method is Nearest. Next you will have to set the calibration file locations. This will be something like:

   settings_url: file:///home/your_user_name/bumblebee2/camera_calibration/bumblebee_right.yaml
   settings_url2: file:///home/your_user_name/bumblebee2/camera_calibration/bumblebee_left.yaml

   Make sure to fill in your user name in the 2 places above. Save and exit.
   Next, we will edit the launch file:

   $ gedit ~/bumblebee2/camera_drivers_experimental/camera1394/launch/startcamera.launch 

   Find and replace all test_camera with bumblebee2.
   Let us launch the stereo camera package now:

   $ roslaunch camera1394 startcamera.launch 

   If everything goes well, you will not see any errors. In a new terminal window, you can view the published topics:

   $ rostopic list

   You should see the following:

   /bumblebee2/left/camera_info
   /bumblebee2/left/image_raw
   /bumblebee2/left/image_raw/compressed
   /bumblebee2/left/image_raw/compressed/parameter_descriptions
   /bumblebee2/left/image_raw/compressed/parameter_updates
   /bumblebee2/left/image_raw/theora
   /bumblebee2/left/image_raw/theora/parameter_descriptions
   /bumblebee2/left/image_raw/theora/parameter_updates
   /bumblebee2/parameter_descriptions
   /bumblebee2/parameter_updates
   /bumblebee2/right/camera_info
   /bumblebee2/right/image_raw
   /bumblebee2/right/image_raw/compressed
   /bumblebee2/right/image_raw/compressed/parameter_descriptions
   /bumblebee2/right/image_raw/compressed/parameter_updates
   /bumblebee2/right/image_raw/theora
   /bumblebee2/right/image_raw/theora/parameter_descriptions
   /bumblebee2/right/image_raw/theora/parameter_updates
   /rosout
   /rosout_agg
   

   In order to view the images, type the following:

   $ rosrun image_view image_view image:=/bumblebee2/left/image_raw

   In a new terminal window, type the following to view the right camera image:

   $ rosrun image_view image_view image:=/bumblebee2/right/image_raw

   You should now be able to see both the images.

Correlation based similarity measures-Summary

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)
Zero-mean Sum of Absolute Differences (ZSAD)
Locally scaled Sum of Absolute Differences (LSAD)
Sum of Squared Differences (SSD)
Zero-mean Sum of Squared Differences (ZSSD)
Locally scaled Sum of Squared Differences (LSSD)
Normalized Cross Correlation (NCC)
Zero-mean Normalized Cross Correlation (ZNCC)
Sum of Hamming Distances (SHD)

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

Left Image

Right 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	SHD Disparity Map
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.

How to select a proper Disparity Value for Stereo-Vision

I have been asked many times how to determine the proper disparity range? If you are taking the images (left, right and ground truth disparity maps) from Middlebury Stereo’s website, then the answer is pretty easy. The minimum disparity value is 0. The maximum disparity value is the maximum value of the pixel in the ground truth map, divided by the scale factor. If you have a stereo vision set-up taking real-world imagery, then you have to do a bit of math to calculate the values.

Here is one of the equations:

In the above equation,

r is the range of the object that you are trying to determine,
b is the baseline of the two cameras, i.e. the distance between the centers of the cameras,
f is the focal length of the image sensor,
x is the pixel size of the image sensor,
N is the maximum disparity value.

For example, lets say that the image sensor has the pixel size of 17um, focal length of 2.8mm, baseline of 28mm, and maximum disparity value is from 5-35. Range vs. disparity values can be plotted in a graphical form as below:

As can be seen above, for the characteristics of the stereo vision system selected, the maximum range of the object that you are trying to detect really depends on the range of disparity values.  For disparity values between 5-35, the detectable range of objects is roughly 15-100 cm.

However, there is another trade-off that people forget about, primarily between the uncertainty in detecting objects at a particular range and the actual range itself. Here is the equation which relates both these variables:

In the above equation, there are two new variables: Δr and ΔN.

Δr is the uncertainty in detecting the object at a certain range,
ΔN is the change in disparity value.

Taking the previous example, lets say that the image sensor has the pixel size of 17um, focal length of 2.8mm, baseline of 28mm, maximum disparity values between 5-35, and the range of 15-100cm, the uncertainty in detecting object at the desired range vs. range values can be plotted in a graphical form as below:

From the above, and for the characteristics of the stereo vision system selected, for a range of 50cm, the uncertainty in the range is 5cm, i.e. if an object has been detected at 50cm, the object could be anywhere from 47.5cm to 52.5cm with a delta of 5cm.

Thus you can see that there are some trade-offs that you should consider while developing a stereo-vision system, and selecting a proper disparity value.

For more information, please refer to:

[1] Khaleghi B., Ahuja S., Wu Q.M.J., “An improved real-time miniaturized embedded stereo vision system (MESVS-II),” IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops, 2008 (CVPRW’08), pp.1-8, 23-28 June 2008.
[2] Khaleghi B., Ahuja S., Wu Q.M.J., “A New Miniaturized Embedded Stereo-Vision System (MESVS-I)”, Canadian Conference on Computer and Robot Vision, 2008 (CRV ‘08), pp.26-33, 28-30 May 2008

Left-Right Consistency (LRC) Check

Left-Right Consistency (LRC) check is performed to get rid of the half-occluded (objects scene in one image, and not in other) pixels in the final disparity map. This is computed by taking the computed disparity value in one image, and re-projecting it in the other image. If the difference in the values is less than a given threshold, then the pixels are half-occluded.

Example

Aloe-Left Image	Aloe-Right Image
Aloe-Ground Truth Left to Right Disparity Map	Aloe-Ground Truth Right to Left Disparity Map
Aloe-Left to Right SAD Disparity Map Window Size: 15×15 Disparity Range: 0-70	Aloe-Right to Left SAD Disparity Map Window Size: 15×15 Disparity Range: 0-70
Aloe-Ground Truth Occlusion Map	Aloe-LRC Map (threshold: 2)

MATLAB Code

% *************************************************************************
% Title: Function-Left/Right Consistency Check
% 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), Threshold for the check (var: thresh) typically 2.0
% Outputs: Disparity Map (var: dispMap), Time taken (var: timeTaken)
% Example Usage of Function: [dispMapLRC, timeTaken]=funcLRCCheck('TsukubaLeft.jpg', 'TsukubaRight.jpg', 9, 0, 16,2);
% *************************************************************************
function [dispMapLRC, timeTaken]=funcLRCCheck(leftImage, rightImage, windowSize, dispMin, dispMax, thresh)
% Initiate the Timer to calculate the time consumed.
tic;
% Perform SAD Correlation based matching (Right to Left)
[dispMapR2L, timeTakenR2L]=funcSADR2L(leftImage, rightImage, windowSize, dispMin, dispMax);
% Perform SAD Correlation based matching (Left to Right)
[dispMapL2R, timeTakenL2R]=funcSADL2R(leftImage, rightImage, windowSize,dispMin , dispMax);
% Prepare matrix for subtraction and scale it for comparison
dispMapL2R=-dispMapL2R;
% Find the size (columns and rows) of the L2R Disparity map and assign the rows to
% variable nrLRCCheck, and columns to variable ncLRCCheck
[nrLRCCheck,ncLRCCheck] = size(dispMapL2R);
% Create an image of size nrLRCCheck and ncLRCCheck, fill it with zeros and assign
% it to variable dispMapLRC
dispMapLRC=zeros(nrLRCCheck,ncLRCCheck);
% 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;
for(i=1:1:nrLRCCheck)
    for(j=1:1:ncLRCCheck)
        xl=j;
        xr=xl+dispMapL2R(i,xl);
        if (xr>ncLRCCheck||xr<1)
            dispMapLRC(i,j) = 0; %% occluded pixel
        else            
            xlp=xr+dispMapR2L(i,xr);
            if (abs(xl-xlp)<thresh)
                dispMapLRC(i,j) = -dispMapL2R(i,j);  %% non-occluded pixel            
            else
                dispMapLRC(i,j) = 0; %% occluded pixel                        
            end
        end
    end
end
% Terminate the Timer to calculate the time consumed.
timeTaken=toc;

Compute Variance Map of an Image

Variance map of an image is calculated by taking a square window of a set size around a center pixel, and calculating the variance of the values of the pixels.

The variance within the window can be calculated from the following equation:

Example

Monopoly-Left Image

Monopoly-Left Image-Variance Map (Pixels marked as black are low-texture regions)

MATLAB Code

% *************************************************************************
% Title: Function-Compute Variance map of the image
% Author: Siddhant Ahuja
% Created: May 2008
% Copyright Siddhant Ahuja, 2008
% Inputs: Input Image (var: inputImage), Window Size (var: windowSize),
% Threshold (var: thresh) Typical value is 140
% Outputs: Variance Map (var: varianceImg) , Time taken (var: timeTaken)
% Example Usage of Function: [varianceImg, timeTaken]=funcVarianceMap('MonopolyLeftColor.png', 5, 9);
% *************************************************************************
function [varianceImg, timeTaken]=funcVarianceMap(inputImage, windowSize, thresh)
try 
    % Grab the image information (metadata) of input image using the function imfinfo
    inputImageInfo=imfinfo(inputImage);
    if(getfield(inputImageInfo,'ColorType')=='truecolor')
    % Read an image using imread function, convert from RGB color space to
    % grayscale using rgb2gray function and assign it to variable inputImage
        inputImage=rgb2gray(imread(inputImage));
        % Convert the image from uint8 to double
        inputImage=double(inputImage);
    else if(getfield(inputImageInfo,'ColorType')=='grayscale')
    % If the image is already in grayscale, then just read it.        
            inputImage=imread(inputImage);
            % Convert the image from uint8 to double
            inputImage=double(inputImage);
        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
    inputImage=inputImage;
end
% Find the size (columns and rows) of the input image and assign the rows to
% variable nr, and columns to variable nc
[nr,nc] = size(inputImage);
% 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
% Create an image of size nr and nc, fill it with zeros and assign
% it to variable meanImg
meanImg=zeros(nr, nc);
% Create an image of size nr and nc, fill it with zeros and assign
% it to variable varianceImg
varianceImg=zeros(nr, nc);
% 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.
% Compute a map of mean values
for(i=1+win:1:nr-win)
    for(j=1+win:1:nc-win)
        sum=0.0;
        for(a=-win:1:win)
            for(b=-win:1:win)
                sum=sum+inputImage(i+a,j+b);
            end
        end
        meanImg(i,j)=sum/(windowSize*windowSize);
    end
end
% Compute a map of variance values
for(i=1+win:1:nr-win)
    for(j=1+win:1:nc-win)
        sum=0.0;
        for(a=-win:1:win)
            for(b=-win:1:win)
                sum=sum+((inputImage(i+a,j+b)-meanImg(i,j))^2);
            end
        end         
        var=sum/((windowSize*windowSize)-1);
        % Apply threshold to produce a binarized variance map
        if (var > thresh)
            varianceImg(i,j) = 255;
        else
            varianceImg(i,j) = 0;
        end
    end
end
% Stop the timer to calculate the time consumed.
timeTaken=toc;

Finding Low-texture or textureless regions in images

Stereo vision algorithms typically compute erroneous results in regions where there is a little or no texture in the scene. They are defined in [1] as regions where the squared horizontal intensity gradient averaged over a square window of a given size is below a given threshold.

Example

Bowling1-Left Image

Bowling1-Left Image-Textureless Map (Pixels marked as white are low-texture regions)

MATLAB Code

% *************************************************************************
% Title: Function-Find Textureless regions of an image
% Notes: Textureless regions are defined as regions where the squared horizontal
% intensity gradient averaged over a square window of a given size 
% (windowSize) is below a given threshold (thresh);
% Author: Siddhant Ahuja
% Created: May 2008
% Copyright Siddhant Ahuja, 2008
% Inputs: Input Image (var: inputImage), Window Size (var: windowSize),
% Threshold (var: thresh) Typical value is 4
% Outputs: Textureless Map (var: texturelessImg) , Time taken (var: timeTaken)
% Example Usage of Function: [texturelessImg, timeTaken]=funcTexturelessRegions('imL.png', 9, 4);
% *************************************************************************
function [texturelessImg, timeTaken]=funcTexturelessRegions(inputImage, windowSize, thresh)
% Read the input image
try
    % Read an image using imread function
    inputImage=imread(inputImage);
    % grab the number of rows, columns, and channels
    [nr, nc, nChannels]=size(inputImage);
     % Grab the image information (metadata) of input image using the function imfinfo
    inputImageInfo=imfinfo(inputImage);
    % Determine if input left image is already in grayscale or color
    if(getfield(inputImageInfo,'ColorType')=='truecolor')
        colored=1;
    else if(getfield(inputImageInfo,'ColorType')=='grayscale')
        colored=0;
        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
    % grab the number of channels
    [nr, nc, nChannels]=size(inputImage);
    if(nChannels)>1
        colored=1;
    else
        colored=0;
    end
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
% Create an image of size nr and nc, fill it with zeros and assign
% it to variable texturelessImg
texturelessImg=zeros(nr, nc);
% Create an image of size nr and nc, fill it with zeros and assign
% it to variable sqGradImg
sqGradImg=zeros(nr,nc);
% 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;
inputImage=double(inputImage);
tic; % Initialize the timer to calculate the time consumed.
% Produce Squared Horizontal Gradient image sqGradImg
for (i=1:1:nr)
    for (j=1:1:nc-1)
        sum = 0.0;        
        for (k=1:1:nChannels)
            diff =  inputImage(i,j,k) - inputImage(i,j+1,k);
            sum = sum + (diff*diff); 
        end
        sum = sum / nChannels;
        sqGradImg(i,j+1) = sum;
        if (j==1)
            sqGradImg(i,j) = sum;
        end
        if (sum > sqGradImg(i,j))
            sqGradImg(i,j) = sum;
        end
    end
end
% Compute average within predefined box window of size windowSize x
% windowSize
for (i=1+win:nr-win)
    for (j=1+win:nc-win)       
        % go over the square window
        sum = 0.0;
        avg = 0.0;
        for (a=-win:1:win)
            for (b=-win:1:win)                
                sum = sum + sqGradImg(i+a,j+b);
            end
        end           
        % Compute the average
        avg = sum / (windowSize*windowSize);
        % Apply threshold
        if (avg < (thresh*thresh))
            texturelessImg(i,j) = 255; % mark detected textureless pixel as white
        end
    end
end
% Stop the timer to calculate the time consumed.
timeTaken=toc;

References

[1] 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.

Finding Depth Discontinuous regions in stereo-images

Stereo vision algorithms typically compute erroneous results in regions where there is a sudden change in the depth between objects in the scene. They are defined in [1] as regions where neighboring disparities differ by more than a certain gap, dilated by a window of a given width.

Example

Flowerpots-Left Image	Flowerpots-Right Image
Flowerpots-Ground Truth Disparity Map=Left-to-Right Image	Flowerpots-Ground Truth Disparity Map-Right-to-Left Image
Flowerpots-Depth-Discontinuous Regions (Pixels marked as white are depth-discontinuous)

MATLAB Code

% *************************************************************************
% Title: Function-Find Discontinuous regions of an image
% Notes: 
% 1. According to paper [A Taxonomy and Evaluation of Dense Two-Frame Stereo Correspondence Algorithms], Depth Discontinuous regions are defined as pixels whose neighboring
% disparities differ by more than a certain gap, dilated
% by a window of width windowSize.
% 2. According to paper [ An Experimental Comparison of Stereo Algorithms], A pixel is a depth 
% discontinuity if any of its (4-connected) neighbors has a disparity that differs by more than 1 from 
% its disparity. Neighboring pixels that are part of a sloped surface can easily differ by 1 pixel, but
% should not be counted as discontinuities.
% Author: Siddhant Ahuja
% Created: May 2008
% Copyright Siddhant Ahuja, 2008
% Inputs: Input Image-Depth Map (var: inputImage), Window Size (var: windowSize),
% Threshold (var: thresh) Typical value is 2
% Outputs: Discontinuous Map (var: discontinuousImg) , Time taken (var: timeTaken)
% Example Usage of Function: [discontinuousImg, timeTaken]=funcDiscontinuousRegions('disp_l_r.png', 9, 2);
% *************************************************************************
function [discontinuousImg, timeTaken]=funcDiscontinuousRegions(inputImage, windowSize, thresh)
% Read the input image
try
    % Read an image using imread function
    inputImage=imread(inputImage);
    % grab the number of rows, columns, and channels
    [nr, nc, nChannels]=size(inputImage);
     % Grab the image information (metadata) of input image using the function imfinfo
    inputImageInfo=imfinfo(inputImage);
    % Determine if input left image is already in grayscale or color
    if(getfield(inputImageInfo,'ColorType')=='truecolor')
        inputImage=rgb2gray(inputImage);
    else if(getfield(inputImageInfo,'ColorType')=='grayscale')
        inputImage=inputImage;
        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
    % grab the number of channels
    [nr, nc, nChannels]=size(inputImage);
    if(nChannels)>1
        inputImage=rgb2gray(inputImage);
    else
        inputImage=inputImage;
    end
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
% Create an image of size nr and nc, fill it with zeros and assign
% it to variable discontinuousImg
discontinuousImg=zeros(nr, nc);
% Create an image of size nr and nc, fill it with zeros and assign
% it to variable edgeImg
edgeImg=zeros(nr, nc);
% 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;
inputImage=double(inputImage);
tic; % Initialize the timer to calculate the time consumed.
% Find variation/edges in disparity values in horizontal and vertical directions
for (i=1:1:nr-1)
    for (j=1:1:nc-1)
        % Traverse in horizontal direction
        if (abs(inputImage(i,j)-inputImage(i,j+1)) > thresh )
            edgeImg(i,j) = 255;
            edgeImg(i,j+1) = 255;
        end
        % Traverse in vertical direction
        if (abs(inputImage(i,j)-inputImage(i+1,j)) > thresh )
            edgeImg(i,j) = 255;
            edgeImg(i+1,j) = 255;
        end
    end
end
% Dilate within the window
for (i=1+win:nr-win)
    for (j=1+win:nc-win)       
        % Go over the square window
        sum = 0.0;
        for (a=-win:1:win)
            for (b=-win:1:win)                
                sum = sum + edgeImg(i+a,j+b);
            end
        end  
        % Apply Threshold
        if (sum>0)
            discontinuousImg(i,j) = 255;
        else
            discontinuousImg(i,j) = 0;
        end
    end
end
% Stop the timer to calculate the time consumed.
timeTaken=toc;

References

[1] 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.

Finding regions that are occluded in left and right images

Stereo vision algorithms typically compute erroneous results when the objects in the scene are fully occluded or half occluded. They are defined in [1] as regions where the left-to-right disparity lands at a location with a larger disparity i.e. regions that are visible in one image, and not in the other image.

Example

Aloe-Left Image	Aloe-Right Image
Aloe-Ground Truth Disparity Map (Left-to-Right)	Aloe-Ground Truth Disparity Map (Right-to-Left)
Aloe-Occluded Regions (Pixels marked as black are occluded)

MATLAB Code

% *************************************************************************
% Title: Function-Find Occluded regions of an image
% Notes: Occluded regions are defined as regions that are occluded in the
% matching image, i.e., where the forward-mapped disparity
% lands at a location with a larger (nearer) disparity.
% Author: Siddhant Ahuja
% Created: September 2008
% Copyright Siddhant Ahuja, 2008
% Inputs: Left to Right Disparity Image (var: dispMapL2R), Right to Left Disparity Image (var: dispMapR2L),
% Scale factor of ground truth map (var: scale) typically 4.0, Threshold for the check (var: thresh) typically 2.0
% Outputs: Disparity Map (var: dispMap), Time taken (var: timeTaken)
% Example Usage of Function: [occludedImg, timeTaken]=funcOccludedRegions('disp_l_r.png', 'disp_r_l.png', 4, 2);
% *************************************************************************
function [occludedImg, timeTaken]=funcOccludedRegions(dispMapL2R, dispMapR2L, scale, thresh)
% Initiate the Timer to calculate the time consumed.
tic;
dispMapL2R=imread(dispMapL2R);
dispMapR2L=imread(dispMapR2L);
dispMapL2R = floor(double (dispMapL2R)/scale);
dispMapR2L = floor(double (dispMapR2L)/scale);
% Prepare matrix for subtraction and scale it for comparison
dispMapL2R=-dispMapL2R;
% Find the size (columns and rows) of the L2R Disparity map and assign the rows to
% variable nrLRCCheck, and columns to variable ncLRCCheck
[nrLRCCheck,ncLRCCheck] = size(dispMapL2R);
% Create an image of size nrLRCCheck and ncLRCCheck, fill it with zeros and assign
% it to variable occludedImg
occludedImg=zeros(nrLRCCheck,ncLRCCheck);
for(i=1:1:nrLRCCheck)
    for(j=1:1:ncLRCCheck)
        xl=j;
        xr=xl+dispMapL2R(i,xl);
        if (xr>ncLRCCheck||xr<1)
            occludedImg(i,j) = 0; %% occluded pixel
        else            
            xlp=xr+dispMapR2L(i,xr);
            if (abs(xl-xlp)<thresh)
                occludedImg(i,j) = 255;  %% non-occluded pixel            
            else
                occludedImg(i,j) = 0; %% occluded pixel                        
            end
        end
    end
end
% Terminate the Timer to calculate the time consumed.
timeTaken=toc;

References

[1] 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.

Correlation based similarity measure-Sum of Hamming Distances (SHD)

Example 1: Tsukuba

Left Image	Right Image
SHD 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

Left Image	Right Image
SHD 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-SHD Right to Left Matching

% *************************************************************************
% Title: Function-Compute Correlation between two images using the 
% similarity measure of Sum of Hamming Distances (SHD) 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]=funcSHDR2L('TsukubaLeft.jpg', 'TsukubaRight.jpg', 9, 0, 16);
% *************************************************************************
function [dispMap, timeTaken]=funcSHDR2L(leftImage, rightImage, windowSize, dispMin, dispMax)
try 
    % Grab the image information (metadata) of left image using the function imfinfo
    leftImageInfo=imfinfo(leftImage);
    % Since SHDR2L 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
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 SHDR2L 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
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;
numberOfBits=8;
tic; % Initialize the timer to calculate the time consumed.
for(i=1+win:1:nrLeft-win)
    for(j=1+win:1:ncLeft-win-dispMax)
        min=0;
        position=0;
        rightWindow=rightImage(i-win:i+win, j-win:j+win);
        for(dispRange=dispMin:1:dispMax)
            if (j+win+dispRange <= ncLeft)
                leftWindow=leftImage(i-win:i+win, j-win+dispRange:j+win+dispRange);
                bloc3=bitxor(rightWindow,leftWindow);
                distance=uint8(zeros(windowSize,windowSize));
                for (k=1:1:numberOfBits)
                    distance=distance+bitget(bloc3,k);
                end
                dif=sum(sum(distance));
                if (dispRange==0)
                    min=dif;
                elseif (min>dif)
                    min=dif;
                    position=dispRange;
                end
            end 
        end
        dispMap(i,j) = position;
    end
end
% Stop the timer to calculate the time consumed.
timeTaken=toc;

Correlation based similarity measure-Normalized Cross Correlation (NCC)

Example 1: Tsukuba

Left Image	Right Image
NCC 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

Left Image	Right Image
NCC 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;