Tutorials about Geometric Transformations, Feature Detection, and Color Processing

Piotr Wendykier, Matthias Odisio

Wolfram Research, Inc.

Abstract:
These tutorials focus on three disciplines of image processing: geometric transformations, feature detection, and color processing. We not only introduce the Mathematica commands in each of these disciplines, but also provide the motivation, the background, and the concepts behind the different image processing techniques. The aim of the tutorial series is to convey a deeper understanding of some built-in image processing functions so that you can utilize these in your day-to-day image processing workflow more effectively.

Topics of Discussion

Geometric Transformations

Padding, cropping, rotating, and resizing

Linear transforms

Nonlinear transforms

Feature Detection

Points

Contours

Faces, barcodes, text

Color Processing

Color primitives

Color spaces

Histogram processing

Geometric Transformations

Mathematica supports complete matrix and functional representations of geometric transformations in any number of dimensions since version 6. In this tutorial we will explain how to apply various transformations to Image and Image3D objects.

Animate[Graphics[{Texture[ImageRotate[pattern=,t,Floor[1/Sqrt[2] ImageDimensions[pattern]]]],GeometricTransformation[Polygon[{{0,0},{1,0},{1/2,Sqrt[3]/2}},VertexTextureCoordinates->{{1/2,Sqrt[3]/2},{0,0},{1,0}}],Join[Table[ReflectionTransform[{Cos[(2 π t)/3+π/2],Sin[(2 π t)/3+π/2]}],{t,0,2}],Table[RotationTransform[(2 π t)/3],{t,0,3}]]]},ImageSize->Medium],{t,0,2 π},SaveDefinitions->True]

Geometric Transformations: Padding, Cropping, Rotating, Resizing

Image Toolbar

Built-in Functions

ImagePad

ImageCrop, ImageTrim, ImageTake

ImageReflect, ImageRotate

ImageResize

ImageCrop with a single argument removes borders of uniform color. Internally it uses BorderDimensions.

ImageCrop[i=]
BorderDimensions[i]

ImageTrim gives the smallest subimage that includes all of the specified points.  Here we trim areas that contain no keypoints.

i=;
xy=ImageKeypoints[i,"KeypointStrength"->10^-5];
HighlightImage[i,xy,"HighlightColor"->Red]
ImageTrim[i,xy,20]

Reflection of a 3D image :

ImageReflect[,#]&/@{Top->Bottom, Top->Right}

Geometric Transformations: Linear Transforms

The purpose of geometric transformations is to deform an image by altering its geometry. A geometric operation transforms a given image I to a new image I’ by modifying the coordinates of image pixels,

I(x,y)→I'(x',y')

i.e., the value of the image function I originally located at (x,y) moves to the position (x’,y’) in the new image I’.  
To model this process, we first need a mapping function

T:(ℝ^2)→(ℝ^2)

that specifies for original 2D coordinate point x=(x,y) the corresponding target point x’=(x’,y’) in the new image I’,

x'=T(x)

A transformed coordinate produced by the mapping function T will, in general, no longer fall into a discrete raster point. The solution to this problem is to compute intermediate pixel values for the transformed image by a process called interpolation.

Affine Mappings

Translation by a vector ((Subscript[d, x]),(Subscript[d, y]))

(Subscript[T, x]):x'=x+(Subscript[d, x])
(Subscript[T, y]):y'=y+(Subscript[d, y])

Scaling by the factors (Subscript[s, x]) and  (Subscript[s, y])

(Subscript[T, x]):x'=(Subscript[s, x])·x
(Subscript[T, y]):y'=(Subscript[s, y ])·y

Shearing by the factors (Subscript[b, x]) and  (Subscript[b, y])

(Subscript[T, x]):x'=x+(Subscript[b, x])·y
(Subscript[T, y]):y'=y+(Subscript[b, y])·x

Rotation by an angle α

(Subscript[T, x]):x'=x·cos(α)-y·sin(α)
(Subscript[T, y]):y'=x·sin(α)+y·cos(α)

The above operations constitute the class of affine transformations. To simplify the concatenation of mappings, it is advantageous to specify all operations in the form of matrix-vector multiplications. Notice that the translation cannot be formulated as matrix-vector multiplication. This difficulty can be resolved with homogeneous coordinates

x=(({{x}, {y}}))           converts  to          (OverHat[x])=(({{hx}, {hy}, {h}}))        with arbitrary h≠0

We can now express any combination of affine transformations as a matrix-vector multiplication

(({{x'}, {y'}, {1}})) =(({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[b, 1])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[b, 2])}, {0, 0, 1}}))·(({{x}, {y}, {1}}))

This 2D coordinate transformation is called an affine mapping, where (Subscript[b, 1]) and (Subscript[b, 2]) specify the translation and (Subscript[a, 11]), (Subscript[a, 12]), (Subscript[a, 21]), (Subscript[a, 22])  aggregate the scaling, shearing, and rotation terms. An affine mapping transforms straight lines to straight lines, triangles to triangles, and rectangles to parallelograms. The distance ratio between points on a straight line remains unchanged after affine mapping.

AffineTransform[{{{(Subscript[a, 1,1]),(Subscript[a, 1,2])},{(Subscript[a, 2,1]),(Subscript[a, 2,2])}},{(Subscript[b, 1]),(Subscript[b, 2])}}]

ImagePerspectiveTransformation can be used to apply affine transformations to Image and Image3D objects.

t=AffineTransform[{{{1,0.2},{0.2,0.5}},{.1,.1}}]
res=ImagePerspectiveTransformation[,t]

The inverse transformation is computed using InverseFunction.

invt = InverseFunction[t]
ImagePerspectiveTransformation[res,invt]

The default value of PlotRange option: {{0,1},{0,height/width}} does not give us the original image back. To include all of the transformed pixels, we need to use PlotRange->All.

{res=ImagePerspectiveTransformation[,t, PlotRange->All],ImagePerspectiveTransformation[res,invt, PlotRange->All]}

Projective Mappings

In contrast to the affine transformation, which provides a mapping between arbitrary triangles, the projective transformation is a linear mapping between arbitrary quadrilaterals. Any combination of projective transformations can be expressed as the following matrix-vector multiplication:

(({{x'}, {y'}, {1}})) =(({{(Subscript[a, 11]), (Subscript[a, 12]), (Subscript[b, 1])}, {(Subscript[a, 21]), (Subscript[a, 22]), (Subscript[b, 2])}, {(Subscript[a, 31]), (Subscript[a, 32]), 1}}))·(({{x}, {y}, {1}}))

Straight lines are preserved as well under this transformations. The projective transformation is the most general transformation that maps any Nth-order algebraic curve to another Nth-order algebraic curve.

ImagePerspectiveTransformation[,LinearFractionalTransform[{{0.9,0.3,0.},{0.3,0.4,-0.1},{.5,0.1,1.}}], PlotRange->All]

Finding Geometric Transforms

Using FindGeometricTransform one can find the mapping between two geometries. That mapping can later be used to register two images. FindGeometricTransform  returns also an estimate of the average alignment error.

{e1,t1}=FindGeometricTransform[i1=,i2=]
ImagePerspectiveTransformation[i2,t1,DataRange->Full]

The same result can be achieved by calling ImageAlign.

ImageAlign[i1,i2]

FindGeometricTransform also accepts sets of corresponding positions which you can determine using ImageCorrespondingPoints.

pts=ImageCorrespondingPoints[i3=,i4=];
MapThread[HighlightImage[#1,#2,Method->{"DiskMarkers",4}]&,{{i3,i4},pts}]

A basic image stitching.

{e2,t2}=FindGeometricTransform[pts[[1]],pts[[2]]];
{w,h}=ImageDimensions[i4];
tmp=ImagePerspectiveTransformation[i4,t2,DataRange->Full,PlotRange->{{0,First@t2[{w,0}]},{0,h}}];
ImageCompose[tmp,{i3,.7},Round@({w,h}/2)]

Both FindGeometricTransform and ImageAlign have the TransformationClass option that allows to specify a geometrical relation between images.  Here we list possible values of the TransformationClass option in order of increasing complexity:

({{"Translation", translation only}, {"Rigid", translation and rotation}, {"Similarity", translation, rotation, and scaling}, {"Affine", linear transformation and translation}, {"Perspective", linear fractional transformation}})

Geometric Transformations: Nonlinear Transforms

A nonlinear mapping cannot be expressed as a simple matrix-vector multiplication in homogeneous coordinates. In general, two types of mappings are possible: source-to-target and target-to-source. In the source-to-target approach, for every pixel position (x,y) in the source image, the corresponding transformed position in the target image is computed

(x',y')=T(x,y)

The problem with this method is that, depending on the geometric transformation T, some elements in the target image may never receive a source pixel value. Conversely, a single pixel in the target image may be hit by multiple source pixels and thus image content may be lost. Due to all of these complications, source-to-target mapping is not the method of choice.

In the target-to-source case, for every discrete pixel position (x', y') in the target image, we compute the corresponding continuous point

(x,y)=(T^( -1))(x',y')

in the source image plane using the inverse geometric transformation (T^( -1)). The major advantage of the target-to-source method is that all pixels in the target image are computed and filled exactly once such that no holes or multiple hits can occur. This mapping, however, requires the inverse geometric transformation (T^( -1)) to be available.

Nonlinear Forward Transform

ImageForwardTransformation implements the source-to-target mapping. This example shows that there are holes in the resulting image.

ImageForwardTransformation[,#^2&,"Interpolated"->False,PlotRange->All]

Nonlinear Inverse Transform

ImageTransformation implements the target-to-source mapping. ImageTransformation computes the inverse mapping (T^( -1)) by calling InverseFunction. Here we show the same transformation (#^2&) computed by the target-to-source mapping.

ImageTransformation[,Sqrt,PlotRange->All, Padding->"Fixed"]

The technique of generating various artistic image effects via nonlinear deformations is called image warping. This is a twirl transformation.

f[pt_]:=With[{r=Norm[pt],a=ArcTan@@pt},
{0.6 r (Sin[(a+12. r)]),0.6 r (Cos[(a+12. r)])}];
ImageTransformation[,f,DataRange->{{-1,1},{-1,1}}]

Here is another transformation called fisheye effect.

f[pt_]:=With[{s={0.5,0.5}},Module[{r,a},
r=((Norm[pt-s])^2)/Norm[s];a=ArcTan@@(pt-s);
s+r{Cos[a],Sin[a]}]]
ImageTransformation[,f]

Feature Detection

Mathematica provides immediate functions for detecting and extracting features in images.  Specific geometrical features such as edges and corners, as well as general keypoints, can be used to register and compare images, or to infer meaningful properties of the 3D world.

Feature Detection: Points

Corners

A corner is defined as an angled bright shape like this one: . Such pattern is does not suffer from the aperture problem, thus it allows to localize a corner accurately and makes it suitable for tracking.

CornerFilter computes a cornerness map at a given spatial scale:

CornerFilter[,3]//ImageAdjust

Corners can be detected with ImageCorners by finding the local maxima of the brightest regions in the cornernerss map. Courtesy of Tom Sherlock:

stars=;

pts=ImageCorners[stars,5];
Short[pts]

That was a lot of corners! We can further control the detection by retaining only the strongest points that are spatially apart:

With[{mindistance=20,n=50},
pts=ImageCorners[stars,5,Automatic,mindistance,MaxFeatures->n]];

HighlightImage[stars,pts,Method->{"CrossMarkers",5}]

Corners can be tracked across images using the function ImageFeatureTrack:

pts=With[{mindistancebetweencorners=10,scale=4},
ImageFeatureTrack[list,ImageCorners[list[[1]],scale,Automatic,mindistancebetweencorners,MaxFeatures->50],"MaxDistance"->mindistancebetweencorners]
];
pts=DeleteCases[pts,_Missing,2,Heads->False];

t=Table[HighlightImage[list[[i]],pts[[i]],Method->{"CrossMarkers",5},"HighlightColor"->Green],{i,Length@list}];

ListAnimate[t,SaveDefinitions->True]

Points of interest

Keypoints are similar to corners in that they are detected based on local spatial geometry, but the detection also considers cornerness across scale, that is by looking at local maxima of cornerness in a 3D space (x, y, σ). Local properties of the neighborhood around each keypoint of an image can be extracted in a fashion robust to change in scale, translation, rotation, and additive luminance.

The function ImageKeypoints detects keypoints using the SURF method.

img=;
points=ImageKeypoints[img,{"Position","Scale","Orientation","ContrastSign"},MaxFeatures->100];
Show[img,Graphics[Table[{{If[p[[4]]==1,Yellow,Red],Circle[p[[1]],p[[2]]*2.5],Line[{p[[1]],p[[1]]+p[[2]]*2.5*{Cos[p[[3]]],Sin[p[[3]]]}}]}},{p,points}]]]

ImageKeypoints can also return a description of the landmark regions of an image as a list of vectors of 64 reals, where each vector is a robust description of the local appearance around each detected keypoint.

desc=ImageKeypoints[img,"Descriptor"];
Dimensions[desc]

Such descriptors can be used for object detection across images, assuming objects do have many points of interest. The function ImageCorrespondingPoints uses the keypoints’ descriptors to determine the correspondences:

pts=ImageCorrespondingPoints[i3=,i4=];

Show[ImageAssemble[{i3,i4}],Graphics[{Red,PointSize[.02],MapThread[If[#2===Missing[],{Red,Point[#1]},Arrow[{#1,#2+{ImageDimensions[i3][[1]],0}}]]&,pts]}]]

Feature Detection: Contours

Edges

In the function EdgeDetect edges are defined according to the traditional approached developed by Canny. An edge is a set of points located at the boundary between image regions, and they are computed by linking pixels where the spatial gradient is high and locally maximal. Edges are detected at a specified scale using spatial derivating kernels optimized according to various methods. In practice, an edge can have an arbitrary shape.

img=;

GradientFilter[img,8]//ImageAdjust

GradientFilter[img,8,"NonMaxSuppression"->True]//ImageAdjust

This one-liner is a compound of the building blocks composing EdgeDetect:

GradientFilter[img,8,"NonMaxSuppression"->True]//MorphologicalBinarize

Manipulate[
EdgeDetect[,radius,Method->m],
{radius,2,20,4},{m,{"Canny","ShenCastan","Sobel"}},SaveDefinitions->True]

With 3D image of a filled cube, what does EdgeDetect detect? Faces!

i=;
Image3D[#,ClipRange->{{15,30},All,10+{15,30}}]&/@{i,EdgeDetect[i]}

Lines

if edges can have arbitrary shapes, lines are straight, geometrical. The function ImageLines detects bright lines on dark background. Typically best results are obtained on binary images, but ImageLines also operates on grayscale images.

cube=;

ecube=EdgeDetect[cube,3,Method->"ShenCastan"]

ImageLines can use two methods: “Hough”, based on the voting mechanism of the Hough transform, and “Ransac”, based on line fitting and outlier detection.

ImageLines has 2 parameters that typically need to be adjusted based on the application: one to control the strength of the detected lines, and one to specify how distinct the lines should be from each other.

Manipulate[
(linescube=ImageLines[ecube,thr,dist];
Show[cube,Graphics[{Thick,White,Line/@linescube}]]),
{{thr,.12},0,1},{{dist,0.04},0,1},SaveDefinitions->True]

Another typical application is the computation of vanishing points for perspective camera calibration. Here the method “Ransac” is used, and only the strongest 6 lines are returned:

img=;
lines=ImageLines[EdgeDetect[img,4,.1],0,0,Method->"Ransac",MaxFeatures->6];
Show[img,Graphics[{Thick,Red,Line/@lines}]]

Feature Detection: Towards High-Level Tasks & Recognition

Face Detection

The function FindFaces uses a classic Viola-Jones face classifier that detects faces as square bounding boxes. In its current version, FindFaces works best with upright faces, and typically yields fast and effective results:

With[{img=RandomImage[]},HighlightImage[img,Graphics[{EdgeForm[{Red,Thick}],FaceForm[],Rectangle@@@FindFaces[img]}]]]

Barcode Recognition

A large number of 1D and 2D barcodes can be recognized with the function BarcodeRecognize, including UPC, Code39, Code128, QR, and PDF417.

image=;

BarcodeRecognize[image]

BarcodeRecognize can provide other varied information about the barcodes in the image:

bc=BarcodeRecognize[image,{"BoundingBox","Data","Format"}]

Show[image,Graphics[{EdgeForm[{AbsoluteThickness[4],Red}], FaceForm[],Rectangle@@bc[[1]], Inset[Style[bc[[2]]<>"\n"<>bc[[3]],Bold,Red,FontSize->20], RegionCentroid[Rectangle@@bc[[1]]]]}]]

Text Recognition

The function TextRecognize performs OCR, and operates with white text on dark background, and vice versa. OCR is available in the languages English, French, German, Italian, Portuguese, Russian, and Spanish. In presence of noise, change in illumination, small characters, gutter deformation, bleeding from the verso, etc., pre-processing typically improves the result:

image=;

TextRecognize[image]

b=LocalAdaptiveBinarize[image,40,{0.85,0.07}]

TextRecognize[b]

Color Processing

Colors are involved in every aspect of our lives. They play an important role in everyday activities such as television, photography and printing. Mathematica has a comprehensive support for device-dependent and device-independent color spaces suitable for color processing and analysis.

flags=;
colors=DominantColors[flags,Automatic,{"Coverage","Color"},Masking->AlphaChannel[flags]];
PieChart[colors[[All,1]],ColorFunction->(#/.Rule@@@colors&),ColorFunctionScaling->False]

Color Processing: Color Primitives and Color Spaces

Mathematica provides a systemwide support for the following color primitives: GrayLevel,  RGBColor, CMYKColor, Hue, XYZColor, LABColor, LCHColor and LUVColor.  These primitives correspond to the following color spaces that characterize every Image and Image3D object: “Grayscale”, “RGB”, “CMYK”, “HSB”, “XYZ”, “LAB”, “LCH”, and “LUV”.

Device-Dependent Color Spaces

The RGB color space is well suited for use in programming, as it is simple to manipulate and maps directly to the typical display hardware. However, the measured distance within RGB color space does not proportionally correspond to our perception of colors. Doubling the value of the red component does not necessarily result in a color that is twice as red.

ConstantImage/@{RGBColor[1/2,0,0],RGBColor[1,0,0]}

In addition, brightness changes are also perceived as nonlinear. Since any coordinate movement modifies color tone, saturation and brightness all at once, color selection in the RGB color space is difficult. In other color spaces, such as HSB, the saturation is represented by one component and can be modified independently.

Manipulate[ImageApply[{hue*#[[1]],saturation*#[[2]],brightness*#[[3]]}&,],{{hue,1},0,1},{{saturation,1},0,1},{{brightness,1},0,1}]

RGB is an additive color system, which means that all colors start with black and are created by adding the primary colors. In contrast, the CMYK color space used for printing is a subtractive color scheme, where each printed color reduces the intensity of the reflected light at that location. Color printing requires at least three primary colors: cyan (C), magenta (M) and yellow (Y). Setting C=M=Y=0 results in the color white (no ink), while setting C=M=Y=1 (complete saturation of all three inks) results in the color black. In practice, the color produced by fully saturating all three inks is not physically a true black. Therefore, the three primary colors are usually supplemented with a black ink (K). Here is an example of simple conversion between RGB and CMYK color spaces.

ColorConvert[RGBColor[0,1,0], CMYKColor]//InputForm

However, this conversion is not sufficient in print production. Special transfer functions stored in color profiles provided with every printer need to be applied to tune the results. Mathematica fully supports ICC v2 and v4 color profiles. This example shows how a CMYK image will look in a print.

rgb = Import["sRGB_v4_ICC_preference.icc"];
ColorConvert[,rgb["PerceptualProfile"]]

The relationship between two colors is deeply influenced by the color model that is used to describe them. The function Blend performs a color mixing in the colors native space thus preserving its natural metric and its characteristics.

Grid[Partition[Labeled[ArrayPlot[{Table[Blend[ColorConvert[{Blue,Red},#],n],{n,0,1,.01}]},AspectRatio->.25,Frame->None],#]&/@
{Hue,RGBColor,LCHColor,LABColor,LUVColor,XYZColor},3],Spacings->{0,1}]

Blend can be also used to change the saturation of the color.

Manipulate[Blend[{i,g},c],{{c,0.5},1,0},Initialization:>(i=;g=ColorConvert[i,"Gray"])]

Device-Independent Color Spaces

All color spaces described above are device-dependent - they are limited by the gamut, i.e. the range of colors that a particular device is capable of. Color systems that describe colors in a device-independent fashion are called colorimetric. Mathematica supports four device-independent color spaces developed by the CIE (International Commission on Illumination): XYZ, L*a*b* LC*h* and L*u*v*. Here we describe in some details only two of them: L*a*b* and L*C*h*. The L*a*b* color model was developed with the goal of linearizing the representation with respect to human color perception. It is a widely used color model for high-quality photographic applications. The dimensions in this color space are the luminosity L* and the two color components a* and b*, which specify the color hue and saturation along the green-red and blue-yellow axes, respectively.

ChromaticityPlot3D["RGB","Lab", BoxRatios->Automatic]

The coordinates in the L*C*h* color model are derived from L*a*b* using a cylindrical transformation that makes them easier to comprehend. The L* component represents luminosity (the same as in L*a*b*), the C* represents chroma or saturation and the h* component represents hue. In this example we perform color adaptation in the L*C*h* space.

img=ColorConvert[,LCHColor];
masks={,,};
ImageCompose[img,MapThread[{alpha,chroma,Δh}|->ImageApply[#{1,chroma,1,1}+{0,0,Δh,0}&,SetAlphaChannel[img,alpha]],{masks,{1.5,1,.2},{.75,0,.4}}]]

The CIE system specifies a number of standard illuminants for a variety of light sources. Each illuminant is described by a spectral radiant power distribution and the correlated color temperature. In practice, color temperature is only meaningful for light sources that correspond closely to the radiation of some black body, i.e., those on a line from reddish/orange via yellow and more or less white to blueish white. Two most important illuminants are D50 and D65. D50 emulates the spectrum of direct sunlight with the color temperature of 5000° K. D65 has a color temperature 6500° K and is designed to emulate the average daylight. The standard illuminants, besides describing the ambient viewing light, also define the reference white points in various color spaces in the CIE color system.

ChromaticityPlot["RGB",WhitePoint->{"D50","D55","D65","D75"}]

Here is an example of how to change the color temperature of an image. Constant k defines the orange direction in the L*a*b* color space.

With[{k=Tan[2 Pi ColorConvert[Orange,LCHColor][[3]]]},Manipulate[ImageApply[# +{0,α,k α}&,ColorConvert[,LABColor]],{{α,0},-.5,.5}]]

Color Processing: Histogram Processing

Histograms are used to depict image statistics in an easily interpreted visual format. By visual inspection of a histogram we can easily conclude if the image is properly exposed. That is why all digital cameras provide histogram overlay on the viewfinder. Poor exposure results in a loss of information that cannot be recover later using image processing.

ImageHistogram[,Appearance->"Separated",FrameTicks->True,ImageSize->Large]

This example illustrates how varying the contrast of an image affects its histogram .

Grid[{imgs ={, ,},ImageHistogram/@imgs, {"Low contrast", "Normal contrast", "High contrast"}}]

This example shows what happens with the histogram when image contains only a small subset of possible pixel values.

Grid[{imgs ={,,},ImageHistogram/@imgs, {"High dynamic range", "Low dynamic range", "Very low dynamic range"}}]

Automatic Contrast Adjustment

ImageAdjust allows to automatically modify the pixel values such that the available range of values is fully covered.

Grid[{imgs ={img=,ImageAdjust[img]},ImageHistogram/@imgs,{"Before adjustment", "After adjustment"}}]

Histogram Equalization

A frequent task is to adjust two different images in such a way that their resulting intensity distributions are similar. The goal of histogram equalization is to find and apply a point operation such that the histogram of the modified image approximates a uniform distribution.

Grid[{imgs ={img=,HistogramTransform[img], HistogramTransform[img,NormalDistribution[.5,.2]]},ImageHistogram/@imgs,{"Before equalization","After equalization (uniform)", "After equalization (normal)"}}]

HistogramTransform can also be used to transfer colors between images. It performs operations in the color space of the input images.

{HistogramTransform[img=,ref=],HistogramTransform[ColorConvert[img,"LAB"],ColorConvert[ref,"LAB"]]}

Questions?