Image Processing and Analysis Quick Start

Shadi Ashnai
Wolfram Research
training@wolfram.com

Image Processing and Analysis Quick Start

Image Objects in Mathematica®

Importing & Creating

Images can be created from 2- or 3-dimensional arrays of data, imported from files stored on the disk, or dragged and dropped into a notebook.

Create a grayscale image from a 2D array:

Image[(({{0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0}, {0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0}, {0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0}, {1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1}})),ImageSize->100]

Create an RGB image from a 3D array:

Image[(({{(({{0.9}, {0.7}, {0.8}})), (({{0.1}, {0.6}, {0.9}})), (({{0.9}, {0.4}, {0.2}})), (({{0.5}, {0.3}, {0.9}}))}, {(({{0.5}, {0.1}, {0.7}})), (({{0.1}, {0.2}, {0.4}})), (({{0.5}, {0.7}, {0.6}})), (({{0.3}, {0.7}, {0.5}}))}, {(({{0.6}, {0.2}, {0.8}})), (({{0.9}, {0.1}, {0.9}})), (({{0.2}, {0.5}, {0.7}})), (({{0.2}, {0.2}, {0.3}}))}, {(({{0.9}, {0.3}, {0.7}})), (({{0.6}, {0.8}, {0.8}})), (({{1}, {0.4}, {0.1}})), (({{0.2}, {0.8}, {0.7}}))}})),ImageSize->100]

Import an image stored on the disk:

Import["ExampleData/peacock.tif"]

Image Dimensions

Use ImageDimensions to get the width and height of an image:

ImageDimensions[]

Note: Dimensions returns rows and columns of the data stored in an image and thus returns reversed dimensions compared to ImageDimensions:

Dimensions[ImageData[]]

Image Coordinates

Mathematica’s standard Image coordinate system (same as Graphics coordinate system):

Pixel indices (same as List coordinate system):

Image Channels

Extract the number of channels that exist in an image:

ImageChannels[]

ImageChannels[]

In a single-channel image, each pixel is stored as a scalar:

ImageData[,"Byte"][[1,1]]

In a multichannel image, each pixel is stored as a list of channel values:

ImageData[,"Byte"][[1,1]]

Image Color Space

The ColorSpace option defines how the values of an image are translated into a color or grayscale representation on the screen. Many color spaces, including grayscale, RGB, HSB, CMYK, and device-independent CIE color spaces Lab, Luv, and XYZ, are currently supported:

Options[]

ImageColorSpace[]

Rendering in different color spaces:

ColorConvert[,"RGB"]

ColorConvert[,"HSB"]

Processing in different color spaces (taking the channel mean):

ImageApply[Mean,ColorConvert[,"RGB"]]

ImageApply[Mean,ColorConvert[,"HSB"]]

Create a an image with an arbitrary, high number of channels (12 channels in this case):

Image[RandomReal[1,{5,5,12}],ImageSize->100]

Note that in this setting ColorSpace->Automatic is used:

Options[%,ColorSpace]

Image Type

Several encoding schemes for pixel values exist in order to compensate for efficient memory consumption versus computation precision. Supported data types and their corresponding storage properties are:

Use ImageType to extract the underlying data type:

ImageType[]

Memory consumption is much less when using only one byte per value:

ByteCount[]

ByteCount[Image[,"Real"]]

Computations, on the other hand, might be more accurate with larger data types. Note that values are truncated and clipped to the interval [0, 1] when using nonreal data types:

Through[{Min,Max}[ImageData[,Automatic]]]

ImageAdd[,0.5]

Through[{Min,Max}[ImageData[%,Automatic]]]

Note: since the values are clipped, even ImageAdjust cannot undo this clipping:

ImageAdjust[%%]

Perform a similar operation on an image of a real type:

ImageAdd[Image[,"Real"],0.5]

Through[{Min,Max}[ImageData[%]]]

ImageAdjust[%%]

Some image processing tends to preserve the data type, while some has to resort to a real type:

Blur[,0.7]

ImageType[%]

GaussianFilter[,0.7]

ImageType[%]

Interactive Image Assistance

Image Processing Capabilities

Pixel Operators

Local Filters

Morphological Operations

Segmentation

Feature Detection & Tracking

Registration

Capabilities—Pixel Operators

Operations including color conversion, level adjustment, brightening and darkening, intensity thresholding, and much more are all in the pixel operation category.

i=;

Result of different built-in operators:

table=Transpose@
{{"grayscale",ColorConvert[i,"Grayscale"]},
{"negation",ColorNegate[i]},
{"threshold",Binarize[i,.2]//ImageAdjust},
{"brightening",ImageAdjust[i,{0,.5}]}};
Row[MapThread[Panel[#2,#1,ImageSize->200]&,table]]

Define an operation that performs selective coloring:

ImageApply[If[1.5#[[1]]>#[[2]]+#[[3]],#,ConstantArray[Mean[#],3]]&,i]

Capabilities—Local Filtering

Many filters are predefined, and many common convolution-based filtering functions are built in. Mathematica's implementations make use of Intel's MKL (Math Kernel Library) and IPP (Integrated Performance Primitives) libraries to maximize performance on supported platforms.

i=;

Result of different built-in filters:

table=Transpose@
{{"Gaussian",GaussianFilter[i,5]},
{"bilateral",BilateralFilter[i,5,.1]},
{"gradient",GradientFilter[i,5]//ImageAdjust},
{"Gabor",GaborFilter[i,5,{1,1}]//ImageAdjust}};
Row[MapThread[Panel[#2,#1,ImageSize->200]&,table]]

Define a custom filter to create a jitter effect:

ImageFilter[RandomChoice[Flatten[#,1]]&,i,5,Interleaving->True]

Application: An Edge-Preserving Smoothing for a Better Segmentation

i=;

p=PeronaMalikFilter[i,20,0.02,2]

c = SelectComponents[Binarize[p], Large];
HighlightImage[i, c]

Capabilities—Morphological Operations

Many basic and advanced morphological operations, from erosion and dilation to skeletonization, are already available in Mathematica.

i=;

Result of some built-in morphological operators:

table=Transpose@
{{"Erosion",Erosion[i,5]},
{"Dilation",Dilation[i,5]},
{"Filling",FillingTransform[i]},
{"Thinning",Thinning[i]}};
Row[MapThread[Panel[#2,#1,ImageSize->200]&,table]]

Represent the connectivity:

MorphologicalGraph[i]

Application: Solve a Maze

Use watershed segmentation to solve a maze:

maze=;

WatershedComponents[maze]//Colorize

res=ColorNegate@Image[WatershedComponents[maze] ,"Bit"]

HighlightImage[maze,Dilation[res,1]]

Capabilities—Segmentation

Basic segmentation algorithms may be based on different properties of the object to be segmented, including its color, shape, position, and more. A perfect segmentation is often achieved either by combining multiple segmentation and other algorithms or by utilizing user interaction.

Binary Segmentation

Mathematica offers a large set of built-in segmentation algorithms based on intensity, connected intensity, color, region growing, active contours, and more.

Intensity thresholding:

Binarize[]

Segment by knowing the color of the background:

ChanVeseBinarize[,"TargetColor"->RGBColor[0,0,.9]]

Region growing:

RegionBinarize[,{{250,220}},.5]

Multi-Region Segmentation

Connected components:

MorphologicalComponents[,.5]//Colorize

Watershed segmentation:

WatershedComponents[FillingTransform[,.9]]//Colorize

Color segmentation:

ClusteringComponents[,3]//Colorize

Capabilities—Feature Detection & Recognition

Detection, recognition, classification, and tracking of features are commonly used in computer vision algorithms. There are many basic algorithms in this category built into Mathematica.

Points of Interest

Keypoints and image alignment:

i1=;i2=;

HighlightImage[#,ImageKeypoints[#,MaxFeatures->Min[ImageDimensions[#]]]]&/@{i1,i2}

ImageCorrespondingPoints[i1,i2]

ImageAlign[i1,i2]

Manipulate[ImageCompose[i1,{RemoveAlphaChannel[ImageAlign[i1,i2],White],a}],{a,0,1}]

Lines of Interest

Straight or curved lines:

image=;

EdgeDetect[image]

lines=ImageLines[EdgeDetect[image],MaxFeatures->5]

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

Objects

Detect faces:

i=;

Show[i,Graphics[{EdgeForm[{Thick,Green}],Transparent,Rectangle@@#}&/@FindFaces[i,{0,36}]]]

Recognition

OCR:

TextRecognize[]

Capabilities—Transformation & Registration

Geometric transforms are used for changing an image perspective, creating effects like looking through a magnifying glass, image alignment for 2D image stitching, or 3D image reconstruction.

i=;

Perspective transformation:

ImagePerspectiveTransformation[i,{{2,0,0},{1,2,0},{1.5,-1/100,1.1}},PlotRange->All,Padding->LightGray]

Custom transformation to create a fisheye effect:

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

ImageTransformation[i,(Subscript[f, {110,120}/200])]

Application: Image Stitching

images={,,,};

Clear[ImageStitch];
ImageStitch[right_,left_]:=Module[{canvas},
canvas=ImageCrop[right,2ImageDimensions[right],Right,Padding->Transparent];
ImageCompose[canvas,ImageAlign[canvas,left]]//ImageCrop
]

ImageStitch[images[[2]],images[[3]]]

p=Fold[ImageStitch,images[[1]],Rest[images]]

ImageRotate[p,-10Degree]

Capabilities—More Tools & Techniques

Measurements via ImageMeasurements:

image=;

Maximum value of each channel:

ImageMeasurements[image,"Max"]

Mean channel values:

ImageMeasurements[image,"Mean"]

Mean pixel intensity:

ImageMeasurements[image,"MeanIntensity"]

Global image entropy:

ImageMeasurements[image,"Entropy"]//N

Color Transformation

Transfer colors between images using HistogramTransform:

img=;ref=;
HistogramTransform[img,ref]

Restoration via Inpainting

inpaintingmethods={"TextureSynthesis","Diffusion"};

img=; mask = ;

Inpaint[img,mask]

Inpaint[img,mask,Method->"Diffusion"]

Texture-Based Segmentation

Sine and cosine parts of GaborFilter when ϕ==0 and ϕ==(π/2):

cosGabor[image_,r_,k_]:=GaborFilter[image,r,k,0]
sinGabor[image_,r_,k_]:=GaborFilter[image,r,k,Pi/2]

Define the magnitude of the Gabor filter as the squared norm of the sine and cosine components:

gaborMagnitude[image_,r_,k_]:=With[{data=ImageData[image]},Image[cosGabor[data,r,k]^2+sinGabor[data,r,k]^2]]

Extract regions with the dominant direction parallel to -45 ° and 45 °:

i=;

λ=.5;r=10;
{t1,t2}=Table[ImageAdjust[gaborMagnitude[i,r,λ k]],{k,{{1,-1},{1,1}}}]

Highlight two texture segments on the image:

ImageCompose[i,{ImageAdd@@MapThread[ImageMultiply[MorphologicalBinarize[#1],#2]&,{{t1,t2},{Red,Green}}],.3}]

Cell Segmentation

i=;

Binarize and fill:

b=ColorNegate[Binarize[i]];
filled=FillingTransform[b]

Use watershed components with markers detected by peaks of distance transform map:

distT=DistanceTransform[filled,Padding->0];
w=WatershedComponents[g=GradientFilter[filled,3],MaxDetect[ImageAdjust[distT],0.02],Method->"Rainfall"];
Colorize[w]

Select cells based on their size and annotate each one with an approximately same-size disk:

cells=SelectComponents[w,"Count",500<#<3000&];
Colorize[cells]

circles=ComponentMeasurements[cells,{"Centroid","EquivalentDiskRadius","Label"}];
Show[i,Graphics[{Blue,Thick,Circle@@#&/@(Most/@circles[[All,2]]),MapThread[Text,{circles[[All,2,3]],circles[[All,2,1]]}]}]]

3D Volumetric Image Processing

3D Volumetric Image Processing (Image3D)

Importing & Creating

3D volumes can be created from a list of imported images, or directly from 3- or 4-dimensional arrays of data. Image3D can be of any supported data type and color space and can have an arbitrary number of channels.

Create a volume from a series of CT scans:

engine=ExampleData[{"TestImage3D","CTengine"}]

Create an RGB volume from a list of images:

Image3D[{,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,}]

Rendering & Exploring

ColorFunction:

Image3D[engine,ColorFunction->"GrayLevelOpacity"]

Image3D[engine,ColorFunction->"RainbowOpacity"]

Image3D[engine,ColorFunction->(Blend[{{0,RGBColor[0,0,1,0]},{0.5,RGBColor[0,0.5,0.25,0.1]},{0.6,RGBColor[1,1,0,0.15]},{1,RGBColor[1,0.5,0,1]}},#1]&)]

Background:

Image3D[engine,ColorFunction->"XRay",Background->Black]

ClipRange:

Image3D[engine,ClipRange->{{128,256},{0,256},{50,1100}}]

Processing

Gradient filtering in 3D:

GradientFilter[engine,2]

High-Performance Image Processing

Image processing functions in Mathematica are already very much optimized using libraries, multithreading, or fast implementations. Nevertheless, better performance may be achieved using GPU implementation of algorithms or out-of-core processing of large images.

CUDALink

CUDA implementation of intensity thresholding:

Needs["CUDALink`"]

binarizeImageCode="__global__ void binarize(unsigned char * in, unsigned char * out, mint threshold, mint width, mint height, mint channels) {
   	mint xIndex = threadIdx.x + blockIdx.x*blockDim.x;
   	mint yIndex = threadIdx.y + blockIdx.y*blockDim.y;
   	mint index = channels*(xIndex + yIndex*width);
   	if (xIndex < width && yIndex < height) {
  		 mint accum = 0;
   		for (mint ii = 0; ii < channels; ii++)
   			accum += in[index+ii];
  		 out[xIndex + yIndex*width] = accum > channels*threshold ? 255 : 0;
   	}
}";

binarizeImage=CUDAFunctionLoad[binarizeImageCode,"binarize",{{"UnsignedByte","Input"},{"UnsignedByte","Output"},_Integer,_Integer,_Integer,_Integer},{16,16}]

input=;
{height,width}=ImageDimensions[input];
channels=ImageChannels[input];
output=CUDAMemoryAllocate["UnsignedByte",{width,height}];

binarizeImage[input,output,150,width,height,channels,{width,height}]

Image[output,ImageSize->Medium]

Compare with built-in Binarize:

Binarize[ImageApply[Mean,input],150/255]

Out-of-Core Processing

When processing very large images that do not fit in the memory, or for processing a very large number of images, subimages or images can be read into the memory, processed, and put back into a file.

ExampleData[{"TestImage","Man"},"ImageSize"]

filename=ExampleData[{"TestImage","Man"},"FilePath"]

Import only a few rows of the image:

Import[filename,"FirstRow"->100,"LastRow"->150]

Out-of-core pixel scanner:

total=0;
ImageFileScan[total+=#&,filename];
total

Out-of-core pixel operator:

ImageFileApply[1-#&,filename,"output1.jpg"]

Thumbnail[Import[%]]

Out-of-core local filter:

ImageFileFilter[Max,filename,10,"output2.jpg"]

Thumbnail[Import[%]]

Thank You!

http://www.wolfram.com

http://www.wolfram.com/imageprocessing

Initialization

Color

Clear[m8red];
m8red[a_]:=RGBColor[0.811765,0.117647,0.145098,a];

C Compiler

Needs["CCodeGenerator`"]

Special Cells

Picture Cell

PicturePrint[expr_,opts:OptionsPattern[Cell]]:=CellPrint[TextCell[expr,"Picture",opts,ShowStringCharacters->False]]

Overview Mouseover

OverviewMouseover[header_,subs_]:=CellPrint@Cell[BoxData[ToBoxes@Mouseover[Framed[Style[header,"OverviewSection"],FrameMargins->{{15,15},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5],

Evaluate@If[subs=="",Framed[Style[header,"OverviewSection",White],FrameMargins->{{15,15},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5,Background->m8red[1]],Row[{Framed[Style[header,"OverviewSection",White],FrameMargins->{{15,15},{5,5}},FrameStyle->None,RoundingRadius->5,Background->m8red[1]]," » ",Framed[Style[subs,"OverviewSection"],FrameMargins->{{15,15},{5,5}},FrameStyle->m8red[1],BoxFrame->2,RoundingRadius->5,Background->White]}]]]],"OverviewSection"]

(*MapThread[OverviewMouseover,Transpose[{{"Chapter One","Chapter One topics, …"},{"Chapter Two","Chapter Two topics, …"},{"Chapter Three","Chapter Three topics, …"},{"Chapter Four","Chapter Four topics, …"}}]]*)

Large Overview Mouseover

LargeOverviewMouseover[header_,subs_]:=CellPrint@Cell[BoxData[ToBoxes@Mouseover[Framed[Style[header,"LargeOverviewSection"],FrameMargins->{{10,10},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5],

Evaluate@If[subs=="",Framed[Style[header,"LargeOverviewSection",White],FrameMargins->{{10,10},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5,Background->m8red[1]],Row[{Framed[Style[header,"LargeOverviewSection",White],FrameMargins->{{10,10},{5,5}},FrameStyle->None,RoundingRadius->5,Background->m8red[1]]," » ",Framed[Style[subs,"LargeOverviewSection"],FrameMargins->{{10,10},{5,5}},FrameStyle->m8red[1],BoxFrame->2,RoundingRadius->5,Background->White]}]]]],"LargeOverviewSection"]

(*MapThread[LargeOverviewMouseover,Transpose[{
{"Image Object","Creation, properties, ..."},
{"Interactive Image Assistance","Image toolbar, Predictive Interface, ..."},
{"Image Processing Capabilities","Filtering, geometry, morphology, ..."},
{"3D Volumetric Image Processing","Creation, rendering, analysis, ..."},
{"High-Performance Image Processing","Out-of-core, CUDALink, ..."}}]];*)

The End