Applications in Computer Vision
using Mathematica 10

Matthias Odisio, Wolfram Research

Join the Conversation #WolframTechConf

Outline

Computer vision in Mathematica 10

Texture Classification

Image Straightening

Depth Map Estimation from Stereo

Computer Vision in Mathematica10

Image Processing

Live image acquisition

Fundamental building blocks: 
Point-based operations (histogram, color), linear and non-linear filters, wavelets, morphology, spatial transformations

Feature detection, matching, and tracking:
Corners, keypoints, edges, lines, ridges

Segmentation: 
Clustering, watershed transform, region growing, graph-based segmentation (Foresting transform), mean shift

Registration: 
Feature-based alignment, dense estimation of geometrical transformation

Computational photography:
Tone mapping, exposure bracketing, focus bracketing, inpainting, background removal

Object detection and recognition: 
Faces, text, barcodes, general identification

Machine Learning

Clustering, nearest neighbors

Classification and Prediction:
Logistic regression, neural network, support vector machine, random forest, ...

Identification based on deep convolutional networks

Dimensionality reduction

Supporting Areas

Geometry: Regions, 2D and 3D geometrical transformations

Statistics: Distributions

Numerical and Symbolic Linear algebra

Optimization

Texture Classification using Haralick’s Features

A set of features defined from the cooccurrence matrix

The cooccurrence matrix represents the probabilities of observing specific pairs of values for neighboring pixels.

For example, in the image below, you can observe that bright pixels have left-right neighbors of typically similar brightness, that dark pixels have typically dark neighbors, and that mid-tone pixels have neighbors of varied intensity.

ImageCooccurrence[,32,{{1,0,1}}]//MatrixPlot

To achieve invariance to image rotation, you can average cooccurrence matrices computed using the eight spatial neighbor configuration (top, top-right, right, etc.).

cooccurrenceMatrix[image_,n_:256]:=With[{kernels={{{0,0,1}},{{1,0,0}},{{1},{0},{0}},{{0},{0},{1}},{{1,0,0},{0,0,0},{0,0,0}},{{0,0,0},{0,0,0},{0,0,1}},
{{0,0,1},{0,0,0},{0,0,0}},{{0,0,0},{0,0,0},{1,0,0}}}},
Mean[ImageCooccurrence[image,n,#]&/@kernels]
];

cooccurrenceMatrix[,32]//MatrixPlot

Haralick has proposed 14 features to describe textures. Each of these features can be computed from the cooccurrence matrix.
Because the 14 features are usually correlated, it is classic to only use the following subset of 5 features:
Homogeneity, Contrast, Correlation, Entropy, and Local Homogeneity.

homogeneity[cm_]:=Total[cm^2,∞];
contrast[cm_]:=Sum[(n-1)^2*Total[Join[Diagonal[cm,n],Diagonal[cm,-n]]],{n,Length[cm]}];
correlation[cm_]:=Module[{miux,mt,m,miuy,stdx,stdy,in=Range[Length@cm]},{m,mt}=Total/@{cm,Transpose[cm]};
{miuy,miux}=Dot[in,#]&/@{m,mt};
{stdy,stdx}={((in-miuy)^2).m,((in-miux)^2).mt};
(Total[KroneckerProduct[in,in]*cm,2]-miux*miuy)/Sqrt[stdx*stdy]];
entropy[cm_]:=With[{m=DeleteCases[Flatten@cm,0|0.]},-Total[m*Log[m]]];
localhomogeneity[cm_]:=Sum[cm[[i,j]]/(1+(i-j)^2),{i,Length@cm},{j,Length@cm}];
textureFeatures[image_/;ImageChannels@image>1]:=Module[
{channels,result,c,s},
channels=ColorSeparate[RemoveAlphaChannel@image,"LAB"];
channels[[2]]=ImageAdjust[channels[[2]],{0,0,1},{-0.8,0.94}];
channels[[3]]=ImageAdjust[channels[[3]],{0,0,1},{-1.13,0.94}];
result=Join[
Table[haralick5[cooccurrenceMatrix[c]],{c,channels}],
Table[haralick5[Mean[{cooccurrenceMatrix[s],cooccurrenceMatrix[Reverse@s]}]],{s,Subsets[channels,{2}]}]
];
Flatten@result
];

haralick5[cm_]:={homogeneity[cm],contrast[cm],correlation[cm],entropy[cm],localhomogeneity[cm]};
textureFeatures[image_]:=haralick5[cooccurrenceMatrix[image]];

textureFeatures[]

What about color images? You can generalize the concept of cooccurence to pairs of images by tallying the occurrences of pixel values pairs where pixels are taken from the two images.  This allows to compute the 5 Haralick features for each of the cross-channel subsets. For example, if you work in the “LAB” color space, you can compute the 5 Haralick features for the pairs: {L, L}, {L, A}, {L, B}, {A, A}, {A, B}, {B, B}.

Following this framework, a color image is described by a feature vector of 30 entries:

textureFeatures[]

Texture Classification using Haralick’s Features

Application: Wood and Bark classification from images using a portion of UIUC texture database

Svetlana Lazebnik, Cordelia Schmid, and Jean Ponce. A Sparse Texture Representation Using Local Affine Regions. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 27, no. 8, pp. 1265-1278, August 2005.

Consider 40 images for each of 6 types of wood and barks:

You will now build a classifier based on the 5 Haralick features for grayscale images.
You can separate this dataset into 6x20 images for training, and 6x20 images for testing.

(*
trainingset=Table[
Association[names[[i,2]]->Table[textureFeatures[Import@img],{img,trainingnames[[i]]}]],
{i,Length@names}];
*)

Get[FileNameJoin[{NotebookDirectory[],"UIUC-texture-training.mx"}]]

The task of the classifier is to correctly interpret a vector of 5 features and assign it to one of the 6 categories:

classifier=Classify[Association@trainingset,PerformanceGoal->"Quality"]

You can now assess the performance of your bark and wood classifier on the test dataset:

Get[FileNameJoin[{NotebookDirectory[],"UIUC-texture-testing.mx"}]]

cm=ClassifierMeasurements[classifier,Association@testingset]

Overall accuracy:

cm["Accuracy"]

You can see from the confusion matrix that the classifier is struggling with the type of wood T04:

cm["ConfusionMatrixPlot"]

A Study of Perspective for Image Straightening

How can you make these images more Euclidean?

A Study of Perspective for Image Straightening

Toy example using synthetic data

image=;

Line detection

Our approach relies on the detection of lines in the image.

First, create an image where lines are the prominent features:

contours=ImageAdjust@GradientFilter[image,2]

Find the straight lines:

lines=ImageLines[contours,0.05,0.05]

Show[image,Graphics[{Thick,Red,Line/@lines}]]

Line intersections

Due to the perspective, parallel lines on the cube actually intersect at the so-called vanishing points.
Vanishing points enable the computation of geometric properties of the image.

One could compute the line intersections with RegionIntersection, but RegionIntersection is too safe for our purpose:

RegionIntersection[Line[{{10,10},{100,10}}],Line[{{20,20},{100,20}}]]

LineIntersectionsPolar[{{xa1_,ya1_},{xa2_,ya2_}},{{xb1_,yb1_},{xb2_,yb2_}}]:=
With[
{denominator = -(xb1-xb2) (ya1-ya2)+(xa1-xa2) (yb1-yb2),
X=xa1 (xb2 (-ya2+yb1)+xb1 (ya2-yb2))+xa2 (xb2 (ya1-yb1)+xb1 (-ya1+yb2)),
Y=xb2 (ya1-ya2) yb1+xa1 ya2 yb1-xb1 ya1 yb2-xa1 ya2 yb2+xb1 ya2 yb2+xa2 ya1 (-yb1+yb2)},
{If[denominator==0,∞,Norm[{X,Y}]/Abs[denominator]],ArcTan[X,Y]}
];

For numerical robustness, it helps to work in a coordinate system with its origin at the image center instead of the bottom-left.

clines=Map[(#-ImageDimensions[image]/2.)&,lines,{-2}]

intersections=Apply[LineIntersectionsPolar,Subsets[clines,{2}],{1}]

Line directions

You can find these intersections that are way outside the image:

directions=Mod[Cases[intersections,{_?(#>Max[ImageDimensions[image]]&),_}][[All,2]]/Degree,180]

There may not be always 3 directions, as an image could exhibit only 1 or 2 vanishing points.

Select the directions that most likely represent the vertical axis:

Select[directions,(60<#<120)&]

We take the median as the most likely candidate for a vertical direction in the image and rotate the image be the difference to 90° to turn it upright.

φ=90-Median[Select[directions,(60<#<120)&]]

ImageRotate[image,φ Degree, Background-> White]

A Study of Perspective for Image Straightening

Two practical examples with human-made structures

Rubik’s cube

image=;

Find the contours in the cube:

contours=ImageAdjust@RidgeFilter[ColorNegate@ColorConvert[image,"Grayscale"],8]

Find the straight lines:

lines=ImageLines[contours,0.18,0.1];
Show[image,Graphics[{Thick,Orange,Line/@lines}]]

clines=Map[(#-ImageDimensions[image]/2)&,lines,{-2}];

Find the intersections:

intersections=Apply[LineIntersectionsPolar,Subsets[clines,{2}],{1}]

Select the intersections corresponding to the vanishing points and get their orientations:

directions=Mod[Cases[intersections,{_?(#>Max[ImageDimensions[image]]&),_}][[All,2]]/Degree,180]

Select the orientations most likely to correspond to the vertical:

Select[directions,(60<#<120)&]

φ=90-Median[Select[directions,(60<#<120)&]]

ImageRotate[image,φ Degree, Background-> White]

Urban environment

image=;

contours=ImageAdjust@GradientFilter[image,1]

lines=ImageLines[contours,0.12,0.05];
Show[image,Graphics[{Thick,Orange,Line/@lines}]]

clines=Map[(#-ImageDimensions[image]/2)&,lines,{-2}];

intersections=Apply[LineIntersectionsPolar,Subsets[clines,{2}],{1}]

directions=Mod[Cases[intersections,{_?(#>Max[ImageDimensions[image]]&),_}][[All,2]]/Degree,180]

Select[directions,(60<#<120)&]

φ=90-Median[Select[directions,(60<#<120)&]]

ImageRotate[image,φ Degree, Background-> White]

Depth Map Estimation from Stereo

Our assumptions about the stereo rig acquiring pictures from the left and the right of a scene, here the 3D point X:

Import["http://www.ni.com/cms/images/devzone/tut/image4573967107314428008.png"]

Depth Map Estimation from Stereo

Images courtesy of
D.Scharstein and R.Szeliski. High-accuracy stereo depth maps using structured light. In IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2003), volume 1, pages 195 - 202, Madison, WI, June 2003.

{left,right}={,};

TabView[{left,right}]

Depth Map Estimation from Stereo: The Fundamental Matrix

Definition

The fundamental matrix expresses the geometrical relationship between stereo images, mapping a point in one image to a line in the other image.
For corresponding points (Subscript[p, 1]) and (Subscript[p, 2]), the fundamental matrix F satisfies (({{(Subscript[p, 2])}, {1}})).F.(({{(Subscript[p, 1])}, {1}}))=0.

Finding sets of corresponding points

We detect SURF keypoints in each image, extract their descriptors of the local appearance at the detected locations, and then match each of these descriptors across the two sets on the left and on the right. This does not involve any geometrical constraint:

{pleft,pright}=ImageCorrespondingPoints[left,right];

The first 100 points in correspondence:

MapThread[Show[#1,Graphics[{Yellow,MapIndexed[Inset[#2[[1]],#1]& ,#2]}],ImageSize->350]&,{{left,right},{pleft[[;;99]],pright[[;;99]]}}]

Depth Map Estimation from Stereo: The Fundamental Matrix

Estimation of the Fundamental matrix

The fundamental matrix verifies (({{(Subscript[p, 2])}, {1}})).F.(({{(Subscript[p, 1])}, {1}}))=0 for each pair of corresponding points. 
We use this formula to cast the problem as a linear algebra task consisting in determining the 8 free parameters of F.

There are various algorithms yielding an estimate for F. We will use RANSAC using the 7-point method:

F=Image`RegistrationDump`FindFundamentalMatrix[pleft,pright,Method->"Ransac"];

MatrixForm[F]

Are we done yet? No. We need dense, ie. for each pixel, information about the disparity between the stereo images.

Depth Map Estimation from Stereo: Rectification

The fundamental matrix is enough for us to derive dense correspondences by scanning lines across the stereo images. 
But that would be too computationally expensive. Instead rectification will warp the images so that corresponding points are on the same raster horizontal line.

Import["http://www.ni.com/cms/images/devzone/tut/lrlwhxvh4348425072667506924.jpg"]

The code is a bit complex. Essentially, we are looking for an homography that allows to make the two images parallel by sending the epipoles to ∞:

ImageStereoRectify[{left_,right_},{pleft_,pright_},F_,center_:Automatic]:=Module[{res,x0,epipoleleft,epipoleright,T,R,G,Hright,M,H0,sol,Ha,Hleft,prighthat,xrighthat,plefthat,a,b,c},
x0=If[center===Automatic,
ImageDimensions[right]/2.,
N@center
];
epipoleleft=First@NullSpace[F];
If[Last@epipoleleft≠0.,epipoleleft/=Last@epipoleleft];
epipoleright=First@NullSpace[Transpose@F];
If[Last@epipoleright≠0.,epipoleright/=Last@epipoleright];
(* If an epipole is in the image,we can not proceeed, since epipoles are to be sent to Infinity *)If[(0≤First@epipoleleft≤First@ImageDimensions@left&&0≤Last@epipoleleft≤Last@ImageDimensions@left)||(0≤First@epipoleright≤First@ImageDimensions@right&&0≤Last@epipoleright≤Last@ImageDimensions@right),Throw@$Failed];

T=TranslationTransform[-x0][[1]];
x0=T.epipoleright;
R=RotationTransform[-ArcTan[x0[[1]],x0[[2]]]][[1]];
x0=R.x0;
G={{1.,0.,0.},{0.,1.,0.},{-1./x0[[1]],0.,1.}};
Hright=G.R.T;
(* Bring epipole at infinity {1,0,0} *)
M=Last@decomposeFundamentalMatrix[F];
H0=Hright.M;
prighthat=Hright.#&/@toP2[pright];
xrighthat=fromP2[prighthat][[All,1]];
plefthat=H0.#&/@toP2[pleft];

Ha={{a,b,c},{0.,1.,0.},{0.,0.,1.}};
sol=NMinimize[Total[(fromP2[Ha.#&/@plefthat][[All,1]]-xrighthat)^2],{a,b,c}];
Ha=Ha/.sol[[2]];
Hleft=Ha.H0;
res={ImagePerspectiveTransformation[left,Inverse@Hleft,PlotRange->All,Padding->"Fixed"],ImagePerspectiveTransformation[right,Inverse@Hright,PlotRange->All,Padding->"Fixed"]};
res
];
decomposeFundamentalMatrix[F_]:=Module[{u,d,v,S,M},
{u,d,v}=SingularValueDecomposition[F];
S=u.{{0,1,0},{-1,0,0},{0,0,0}}.Transpose@u;
M=u.{{0,-1,0},{1,0,0},{0,0,1}}.ReplacePart[d,{3,3}->(d[[1,1]]+d[[2,2]])/2].Transpose[v];
{S,M}];

toP2[x_]:=Join[x,ConstantArray[1.,{Length@x,1}],2];
fromP2[x_]:=Most[#]/Last[#]&/@x;

{rectleft,rectright}=ImageStereoRectify[{left,right},{pleft,pright},F];

d=Min/@Transpose[ImageDimensions/@{rectleft,rectright}];
{rectleft,rectright}=ImageTake[#,{1,d[[2]]},{First@ImageDimensions@#-d[[1]]+1,-1}]&/@{rectleft,rectright};

TabView[{rectleft,rectright}]

Corresponding points are aligned horizontally indeed:

Show[a=ImageAssemble[{rectleft,rectright}],Graphics[{Yellow,Thick,Line[{{0,#},{First@ImageDimensions@a,#}}]&/@Range[20,400,50]}]]

Depth Map Estimation from Stereo: Disparity

Concept of dense matching for disparity estimation

For each point in, say, the left image, we can scan the corresponding line in the right image and find corresponding pixel for which the local patch bests matches the local patch around the point in the left image.

For example, for this point in the left image, we are looking at the best local match in the band from the right image:

index={50,353};

patchradius=4;
searchradius=20;

patchleft=ImageTake[rectleft,index[[1]]+patchradius*{-1,1},index[[2]]+patchradius*{-1,1}]

bandright=Image[ImageTake[rectright,index[[1]]+patchradius*{-1,1},index[[2]]+searchradius*{-1,1}],Magnification->8]

xcor=ImageCorrelate[ColorConvert[bandright,"Grayscale"],ColorConvert[patchleft,"Grayscale"],EuclideanDistance, Padding->None]

δ=Flatten@ImageData@xcor;

ListLinePlot@δ

This is our estimate of the disparity using the left image as reference:

d=(Length[δ]+1)/2-First@Ordering[δ,1]

Depth Map Estimation from Stereo: Disparity

Estimation of the disparity (displacement, and best match value)

We can estimate the disparity for all the pixels at once, if we are willing to accept bad results near the boundary of the image.

This code estimates the disparity considering one image as fixed and the other one being shifted:

disparity[imgfixed_,imgmoving_,patchradius_,searchradius_,gradientweight_:0.]:=
Module[
{imgMinDist,imgDisparity,shifted,tx,ssd,better,gfixed,gshifted},
imgMinDist=ConstantImage[$MaxMachineNumber,ImageDimensions@imgfixed];
imgDisparity=ConstantImage[0,ImageDimensions@imgfixed];
If[gradientweight>0,
gfixed=GradientFilter[imgfixed,2]
];

Do[
shifted=ImagePerspectiveTransformation[imgmoving,TranslationTransform[{tx,0}],ImageDimensions@imgfixed,DataRange->Full,Padding->"Fixed"];
ssd=MeanFilter[ColorDistance[shifted,imgfixed],patchradius];
If[gradientweight>0,
ssd=ImageAdd[ssd,ImageMultiply[MeanFilter[ColorDistance[GradientFilter[shifted,2],gfixed],patchradius],gradientweight]]
];
better=Binarize[ImageSubtract[imgMinDist,ssd],0];
imgMinDist=ImageAdd[ImageMultiply[imgMinDist,ColorNegate@better],ImageMultiply[ssd,better]];
imgDisparity=ImageAdd[ImageMultiply[imgDisparity,ColorNegate@better],ImageMultiply[ConstantImage[tx,ImageDimensions@imgfixed],better]],
{tx,-searchradius,searchradius}];
{imgMinDist,imgDisparity}
]

Let’s compute the left disparity and the right disparity:

{mindistLeft,disparityLeft}=disparity[ColorConvert[rectleft,"Grayscale"],ColorConvert[rectright,"Grayscale"],2,10,5];

{mindistRight,disparityRight}=disparity[ColorConvert[rectright,"Grayscale"],ColorConvert[rectleft,"Grayscale"],2,10,5];

There can be cases where the algorithm will fail, for example in case of occlusions. 
We remove false detections by enforcing similar results in the left and right disparity maps:

maskdisp=Binarize[ImageDifference[disparityLeft,ImageMultiply[disparityRight,-1]],{0,2}]

Finally we consolidate the result by merging the two disparity maps by selecting the best match for each pixel:

maskrightisbetter=Binarize[ImageSubtract[mindistLeft,mindistRight],0];

consolidatedMinDist=ImageMultiply[
ImageAdd[ImageMultiply[mindistRight,maskrightisbetter],ImageMultiply[mindistLeft,ColorNegate@maskrightisbetter]],
maskdisp];

consolidatedDisparity=ImageMultiply[
ImageAdd[ImageMultiply[ImageMultiply[disparityRight,-1],maskrightisbetter],ImageMultiply[disparityLeft,ColorNegate@maskrightisbetter]],
maskdisp];

Voila! We are not doing too well on the foreground:

c=Colorize[ImageAdjust@consolidatedDisparity,ColorFunction->ColorData["TemperatureMap"]]

Going further

General Resources

Computer Vision [guide/ComputerVision]

Image Processing [guide/ImageProcessing]

Machine Learning [guide/MachineLearning]

Questions?