What’s New in Image & Signal Processing

Shadi Ashnai, WRI

Join the Conversation #WolframTechConf

In this talk

Audio Processing — new in-core and out-of-core, plus a lot of processing and analysis capabilities

Signal Processing — irregularly sampled time series, canonical time warping, ...

Image Processing — photography, optical flow, mesh conversion

Integration with Neural Networks — tackle the problems in a modern way

Audio [new]

The Wolfram (Language^™) allows for immediate representation of an Audio object, whether it is in-memory or referencing a local or remote audio file.

New audio capabilities covering audio synthesis, processing, and analysis combined with comprehensive support for signal processing, statistical analysis, and machine learning enable easy prototyping of applications in various fields, including music, speech, communication, and broadcasting.

In-Core & Out-of-Core

Import["ExampleData/rule30.wav"]

Audio["ExampleData/rule30.wav"]

Audio["http://exampledata.wolfram.com/bach.mp3"]

What to do with them

Generate, visualize, add effects, process and analyze audio signals using tens of audio-specific functions as well as hundreds of filters and arithmetic operations.

a=AudioGenerator[{"Sin",TimeSeries[440. (2^((50+RandomInteger[24,50]-69)/12)),{0,10}]}]

Spectrogram[a,PlotRange->{All,{0,1000}}]

MatrixPlot[Transpose[AudioLocalMeasurements[a,"MFCC"]["Values"]],PlotTheme->"Minimal"]

MatrixPlot[DistanceMatrix[AudioLocalMeasurements[a,"MFCC"]["Values"]],PlotTheme->"Minimal"]

AudioReverb[a,ExampleData[{"Audio","IRStMarysChurch"}]]//AudioNormalize

[Workshop] Introduction to Audio in the Wolfram Language — Today 1-2pm

[TALK] Audio Effects — Tomorrow 9:30-10am

[TALK] Audio Synthesis — Tomorrow 10-10:30am

[TALK] Audio Analysis — Tomorrow 10:30-11am

Signal [dynamic and canonical time warping]

Compare Stock Prices Using DTW

Use WarpingCorrespondence to compare the first quarter of 2016 of the HPQ stock prices with historical data from 2010 to 2015.

recent=FinancialData["HPQ",{{2016,1,1},{2016,3,31}},"Value"];
{histDates,hist}=Transpose[FinancialData["HPQ",{{2010,1,1},{2015,1,31}},"DateValue"]];

Find the best matching subsequence of historical data.

{corrHist, corrRecent}=WarpingCorrespondence[hist,recent,Method->{"MatchingInterval"->"Flexible"}];

Detect the historical interval most similar to quarter one of 2016.

{m,n}=corrHist[[{1,-1}]];
histDates[[{m,n}]]

Visually compare recent data and the best historical match.

DateListPlot[{AssociationThread[Take[histDates[[m;;n]], Length[recent]], recent], AssociationThread[histDates[[m;;n]], hist[[m ;; n]]]}, PlotTheme->"Business",PlotLegends->{"2016 data", "Best historical match"}, DateTicksFormat->{"MonthNameShort", " ", "YearShort"}]

Predict the stock prices for the next 30 days based on historical data.

l=Length[recent];
colDat=ColorData["Atoms"];
offset=Last[recent]-hist[[n]];
hist30d=hist[[n;;n+30]];
ListLinePlot[{recent,{l+Range[31],hist30d+offset}\[Transpose], hist[[m ;; n]],{n-m+Range[31],hist30d}\[Transpose]},PlotStyle -> {colDat["Rh"], {colDat["Mo"],Dotted},colDat["Yb"],{colDat["Tb"],Dotted}},PlotLegends ->Placed[{"2016 data", "30 days prediction", "Best historical match","Next 30 days of historical data"},Below],Axes->{False,True},Epilog->{Dashed,{Arrow[{{n-m,hist[[n]]},{l+1,Last[recent]}}],Arrow[{{n-m+30,hist[[n+30]]},{l+31,hist[[n+30]]+offset}}]}}]

Canonical Time Warping (CTW)

Canonical time warping (CTW) iteratively performs spatial transformations and dynamic time warping on the two sequences to find the alignment with minimal distance between sequences.
The sequences can be lists of numeric scalars or vectors. In contrast to dynamic time warping, elements can also be of different dimensions.

ref=Table[{Cos[t],Sin[t],0.5t},{t,0,3Pi,3Pi/100}];
Graphics3D[Point[ref]]

query=Table[{Cos[t],Sin[t+RandomReal[.2]]},{t,Pi,2Pi,2Pi/100}];
ListPlot[query]
CanonicalWarpingDistance[ref,query]

query=Table[{t,Sin[5t]},{t,Pi,2Pi,2Pi/100}];
ListPlot[query]
CanonicalWarpingDistance[ref,query]

Canonical time warping is translation, scale and rotation invariant.

s1=Table[{Cos[3t]*Cos[t],Cos[3t]*Sin[t]},{t,0,2Pi,2Pi/100}];

s2=TranslationTransform[{1,1}]/@s1;
GraphicsRow[{With[{s2=TranslationTransform[{1,1}]/@s1},ListLinePlot[{s1,s2},PlotLabel->"(translation) ctw: "<>ToString[Chop[CanonicalWarpingDistance[s1,s2]]],ImageSize->Medium]],
With[{s2=2 s1},ListLinePlot[{s1,s2},PlotLabel->"(scaling) ctw: "<>ToString[Chop[CanonicalWarpingDistance[s1,s2]]]]],
 With[{s2=RotationTransform[Pi]/@s1},ListLinePlot[{s1,s2},PlotLabel->"(rotation) ctw: "<>ToString[Chop[CanonicalWarpingDistance[s1,s2]]]]]}]

Compare the Shape of Different Countries

Use CTW to find similarity between shapes of different countries.

names={"USA","UK","Italy","Poland","Portugal","Vietnam","Brazil","Finland"};
borders=<|#->ArrayResample[CountryData[#,"Polygon"][[1,1,1]],200]&/@names|>;

shapes=Graphics[CountryData[#, "Shape"][[1]], ImageSize -> {50, 50}] & /@ Keys[borders]

Calculate and show shape distances.

dm=DistanceMatrix[Values[borders],DistanceFunction->(Chop[CanonicalWarpingDistance[##,MaxIterations->5]]&)];

MatrixPlot[dm,FrameTicks->{{MapIndexed[{#2[[1]],#1}&,shapes],MapIndexed[{#2[[1]],Rotate[#1,π/4]}&,names]},{MapIndexed[{#2[[1]],Rotate[#1,π/4]}&,names],MapIndexed[{#2[[1]],#1}&,shapes]}},Mesh->True,PlotLegends->Automatic,ImageSize->Large]

Time Series Processing

direct arithmetic and statistic operations on time series

extended support for irregularly sampled time series

a=Import["ExampleData/rule30.wav"]

t1=AudioLocalMeasurements[a,"RMSAmplitude",PartitionGranularity->{.05,.02}]

t2=AudioLocalMeasurements[a,"SpectralFlatness",PartitionGranularity->{.2,.2}]

ListPlot[{t1,t2,t1/10+2t2}]

[TALK] Time Series Processing — Tomorrow 1-1:30pm

Image [more of]

Computational Photography — photo manipulations, stacking, metadata support

Video Processing (the start) — more advanced analysis of image sequences including optical flow

Large Images (the start) — viewing large out-of-core images

Integration with Computational Geometry — convert an image to a mesh region

Fully Upgraded HighlightImage

The new HighlightImage extends the powerful syntax of graphics primitives and directives to give you maximal freedom.

i=;
HighlightImage[i,{"Boundary",EdgeForm[{Dashed,Thick,Green}],Rectangle@@@FindFaces[i]}]

HighlightImage[i,{EdgeForm[{Dashed,Thick,Green}],FaceForm[Purple],Rectangle@@@FindFaces[i]}]

HighlightImage[i,{EdgeForm[{Dashed,Thick,Green}],FaceForm[Purple],Rectangle@@@FindFaces[i],Yellow,EdgeDetect[i,10],Red,ImageCorners[i,MaxFeatures->30]}]

HighlightImage[i,{EdgeForm[{Dashed,Thick,Green}],Rectangle@@@FindFaces[i]},"Desaturate"]

HighlightImage[i,{ImageMarker[Rectangle@@@FindFaces[i],SetAlphaChannel[LinguisticAssistant,Graphics[{White,Disk[]},Background->Black]]]}]

Computational Photography: Photo Manipulation

Correct for Underwater Color Cast

Underwater photography is usually affected by a strong blue or green color cast. 
The process of correcting these unrealistic colors is called color balance (or white balance). 
Correct a strong blue color cast, assuming there should be some white in the picture.

TabView[Map[ImageCollage[{#,ColorBalance[#]},Background->None,ImagePadding->5]&,{,,,}]]

Simulate Different Time of Day Using ColorBalance

Sunlight changes color throughout the day as sun rays travel through the planet's atmosphere for varying durations. You can use ColorBalance to simulate this change in the scene illumination.

Manipulate[ColorBalance[,->(Quantity[temperature, "K", kelvins, "Kelvins"])],{temperature,3800,7000}]

Perform Color Tone Mapping

Tone mapping can be used to represent a high dynamic range (HDR) image on a low dynamic range screen by compressing the range of luminance values while preserving local spatial details.

hdr=;

TabView[{hdr,ColorToneMapping[hdr]}]

Combine Flash and No-Flash Photography

Flash is great to capture details in low-light conditions, but it can alter rather dramatically the colors in a scene.
Enhance a faithful—albeit noisy—picture using GuidedFilter to smooth the noise while keeping the details.

noflash=;flash=;
TabView[{"No-Flash"->noflash,"Flash"->flash,"Combined"->GuidedFilter[noflash,flash,10,.001]}]

Photography: Photo Stacking

In macro photography or microscopy images the depth of field (DOF) is generally very small. 
This makes it impossible to capture a single image where objects at different distances are all sharp.
With ImageFocusCombine, you can combine several images at different focal distances to generate an all-sharp image.

Source: Copper crystal 4 mm in edge © Honorio Cocera (link)

list=Import/@FileNames["*",FileNameJoin[{NotebookDirectory[],"images"}]];
ListAnimate[list]

ImageFocusCombine[list]

Photography: Metadata Support

Using pictures and their metadata, you can perform geo-computations to visualize and analyze a walking tour of Lisbon, Portugal.

Import and analyze a collection of photos taken during a visit.

files=ExampleData[#,"FilePath"]&/@ExampleData[{"TestImageSet","Lisbon2016"}]

Gather all the information required in an association.

labels={"Thumbnail","DateTime","GeoPosition"};

dataset=Dataset@Map[file|->AssociationThread[labels->Import[file,{labels},IncludeMetaInformation->False]],files]

Visualize locations where images where taken.

sorted=Normal@dataset[SortBy["DateTime"],Take[#GeoPosition,All,2]&];
TravelDistance[sorted,TravelMethod->"Walking"]

GeoGraphics[GeoMarker@sorted,GeoRangePadding-> (Quantity[100, "m", meters, "Meters"])]

GeoGraphics[{GeoMarker@sorted,Thick,Blue,Line[gps]},GeoRangePadding-> (Quantity[100, "m", meters, "Meters"])]

[TALK] Image Collections — Today 4:30-5pm

Large Images (view, compare, measure)

DynamicImage

img=ExampleData[{"TestImage","Volubilis"}];

DynamicImage[img, ImageSize-> {640,480}, ZoomCenter->{Scaled[Dynamic[xx]],Scaled[Dynamic[yy]]}]

Slider[Dynamic[xx]]

Interactive Image Measurements

DynamicImage facilitates the construction of custom interfaces to interactively measure objects in images.

Consider the measurement of a cell radius r.

img0=;

measurementCircle[Dynamic[center_],Dynamic[r_]]:={Yellow,Thick,Dynamic@Circle[center,r],Dynamic@Text[Style[StringJoin[ToString@Round@r,"px"],FontSize-> 18],Scaled@{0.1,0.9},Background-> RGBColor[0,0,0,0.5]]};

DynamicModule[
{center=ImageDimensions[img0]/2,r},
Manipulate[
EventHandler[
DynamicImage[
img0,
Epilog-> measurementCircle[Dynamic[center],Dynamic[r]]
],
{"MouseDown":>If[CurrentValue["OptionKey"],center=MousePosition["Graphics"]]},
PassEventsDown-> Dynamic[Not[CurrentValue["OptionKey"]]]
],
{{r,32},8,800},
FrameMargins->0
]
]

Dynamic Comparison of Two Images

Compare two satellite images of Venice, Italy, using an interactive graphical user interface based on two coupled DynamicImage objects. The left display shows a radar image taken by the ESA satellite Sentinel-1 on May 10, 2016; the right is an optical photograph taken by Sentinel-2 on April 29, 2016.

{img1,img2}={
Import[FileNameJoin[{NotebookDirectory[],"Venice-Sentinel-1.jpg"}]],
Import[FileNameJoin[{NotebookDirectory[],"Venice-Sentinel-2.jpg"}]]
};
TabView[{img1,img2}]

DynamicModule[
{x,y,s},
Row[{
DynamicImage[
img1,
ZoomCenter->Dynamic[Scaled[{x,y}]], 
ZoomFactor-> Dynamic[s],
AppearanceElements->{"Pan","Scrollbars","ZoomButtons","Zoom"},
ImageSize-> {360,540}
],
DynamicImage[
img2,
ZoomCenter->Dynamic[Scaled[{x,y}]],
ZoomFactor-> Dynamic[s],
AppearanceElements->{"Pan","Scrollbars","ZoomButtons","Zoom"},
ImageSize->{360,540}
]
}]
]

Optical Flow Detection & Analysis

ImageDisplacements can capture dynamic features in video sequences. From one video frame to the next, the displacement is a vector field (OverVector[M]) that describes the optical flow that is essentially the 2D motion field of objects in the video.

Vector field analysis provides a variety of tools to analyze these motion fields.

Rotation Detection

seq = Image3DSlices[];

A video sequence of rotating wind turbines.

ListAnimate[seq]

Obtain the average optical flow of the entire video sequence.

flow=Mean@ImageDisplacements[Append[seq,First[seq]]];

Show[
First[seq],
ListVectorPlot[Transpose@Reverse@flow,VectorStyle->Yellow,VectorPoints-> Fine]
]

The curl ∇×(OverVector[M]) of a vector field (OverVector[M]) determines the amount of inherent rotation.

listCurl=
flow  |-> ListConvolve[{{0.5`,0,-0.5`}},flow[[All,All,2]],{{1,2},{1,2}},"Fixed"]-ListConvolve[{{-0.5`},{0.`},{0.5`}},flow[[All,All,1]],{{1,2},{1,2}},"Fixed"];

curl=ImageAdjust@Image@(listCurl@flow)

A blurred curl distribution provides a more robust result for large-scale rotations. Black stands for negative rotation; white denotes rotation with positive orientation.

blurredCurl=ImageAdjust[GaussianFilter[curl,8{3,1}]]

Colorize the rotation map and superimpose it on top of the original image.

ImageCompose[seq[[1]],SetAlphaChannel[Colorize[blurredCurl,ColorFunction->ColorData[{"TemperatureMap","Reverse"}]],ColorNegate[blurredCurl]]]

Live Video Analysis

preProcessFrame=
img|->ImageReflect[ColorConvert[img,"Grayscale"],Left->Right];

The entire code of this augmented reality example is less than one page.

DynamicModule[
{frames,flow,curl,mask,time=0,W= 320,H=240,w=240,h=160,r=64},
(* specify a sensitive image region *)
mask=ImageCompose[
ConstantImage[0,{W,H}],
ImageAdjust[Image[GaussianMatrix[r/2]]],
{w,h}
];
(* initialize the optical flow acquisition *)
frames=Table[preProcessFrame@CurrentImage[],{2}];
flow=First@ImageDisplacements[frames];
Graphics[
{Inset[
(* read and display camera image, determine the curl *)
Dynamic[
frames[[1]]=frames[[2]];
frames[[2]]=preProcessFrame@CurrentImage[];
flow=First@ImageDisplacements[frames,0.7flow];
curl=Image@listCurl@flow;
time-= 2 ImageMeasurements[curl,"Mean",Masking->mask];
frames[[2]]
],
{0,0},{0,0},{W,H}
],
Inset[
(* display clock *)
Dynamic[
ClockGauge[
DatePlus[Date[],{Round[5 time,5],"Minute"}],
GaugeStyle-> ,
LabelStyle-> White,
GaugeFaceStyle-> Opacity[0.1]
],
UpdateInterval->1
],
{w,h},{0,0},{2r,2r}
]
},
PlotRange-> {{0,W},{0,H}},
ImageSize-> {W,H}
]
]

Image to Mesh Conversion

Waves in the Great Salt Lake

Using ImageMesh, you can convert image segments into BoundaryMeshRegion objects. These mesh regions allow you to utilize functions across other areas, such as finite element methods (FEM).

To illustrate this FEM connection, determine the major modes of surface waves on the Great Salt Lake in Utah.

img=EntityValue[(Entity["Great Salt Lake", Entity["Lake","GreatSaltLake::yw8cf"], "Entity[\"Lake\", \"GreatSaltLake::yw8cf\"]", "lake"]),"Image"]

Regularize the image with a mean-shift filter.

img2=MeanShiftFilter[img,3,0.1]

Obtain a segmentation via region growing.

mask=RegionBinarize[img2,,1/5]

Convert to a mesh.

ℛ=ImageMesh[mask]

Triangulate.

Ω=TriangulateMesh[ℛ,MaxCellMeasure-> 8]

Solve the wave equation of the lake surface by first determining the
eigenfunctions of the Laplacian inside the lake region.

ℒ = -(Subscript[∇, {x,y}]^2)φ[x,y];

ℬ = DirichletCondition[φ[x,y]==0,True];

{Λ,Φ}=NDEigensystem[{ℒ,ℬ},φ[x,y],{x,y}∈Ω,64];

Display the first six oscillation modes.

GraphicsGrid[
Partition[
Table[ContourPlot[Φ[[k]],{x,y}∈Ω,PlotRange->All,PlotLabel->Λ[[k]],PlotTheme->"Minimal"],{k,6}],
3
],
ImageSize-> 512
]

Multi-View Reconstruction (Image + Geometry + Optical Flow)

A 3D object can be reconstructed from multiple 2D views. ImageDisplacements is used to determine the parallax from one view to the next. The larger the parallax shift, the closer the object behind the corresponding pixel. With this depth information, you can extrude the vertices of an object mesh obtained through ImageMesh and TriangulateMesh. The result is displayed via texture mapping on the 3D mesh object.

imgs={,,};

Obtain an object mask of the central image.

mask=Erosion[Binarize[imgs[[2]],0],1];

Determine the parallax in the left and right images with respect to the center image.

parallaxL=First@ImageDisplacements[imgs[[{2,1}]]];

parallaxR=First@ImageDisplacements[imgs[[{2,3}]]];

Combine these displacements in a total parallax.

parallax=parallaxL-parallaxR;

In the given setting, the x component of the parallax is roughly proportional to the depth of the respective pixel source.

depth = Blur@Opening[ImageMultiply[ImageAdjust@Image[parallax[[All,All,1]]],mask],DiskMatrix[4]]

Construct a depth function.

depthFunction = ListInterpolation[Transpose@Reverse@ImageData[depth]];

Obtain an object mesh with a refinement that depends on the object's saliency.

resolution=ImageAdjust@ImageSaliencyFilter[imgs[[2]]];

resolutionFunction = ListInterpolation[Transpose@Reverse@ImageData@resolution];

Ω=TriangulateMesh[
ImageMesh[Erosion[mask,DiskMatrix[2]]],MeshRefinementFunction->Function[{vertices,area},area>32+512(1-resolutionFunction@@Mean[vertices])^6]
]

Extrude the object mesh with the depth function.

Graphics3D[
GraphicsComplex[
Apply[{##,depthFunction[##]}&,MeshCoordinates[Ω],{1}],
{EdgeForm[],MeshCells[Ω,2]}
],
PlotRange-> Append[Thread[{0,ImageDimensions[mask]}],{0,1}],
BoxRatios-> {1,1,2/3},
ViewPoint->Top
]

Extract the object texture from the central image.

texture = SetAlphaChannel[imgs[[2]],mask]

Map the texture onto the 3D object.

Graphics3D[
{Texture[texture],
GraphicsComplex[
Apply[{##,depthFunction[##]}&,MeshCoordinates[Ω],{1}],
{EdgeForm[],MeshCells[Ω,2]},
VertexTextureCoordinates-> Map[#/ImageDimensions[texture]&,MeshCoordinates[Ω]]
]
},
PlotRange-> Append[Thread[{0,ImageDimensions[mask]}],{0,1}],
BoxRatios-> {1,1,2/3},
Lighting-> {{"Ambient",White}},
ViewPoint-> Top,
Boxed->False
]

Images + Deep Neural Networks

Tackle problems including large collections of images in a modern way.

Digit Classification

Use the MNIST database of handwritten digits to train a convolutional network to predict the digit given an image.

First obtain the training and validation data.

resource = ResourceObject["MNIST"];
trainingData=ResourceData[resource,"TrainingData"];
testData=ResourceData[resource,"TestData"];

RandomSample[trainingData,5]

Define a convolutional neural network that takes in 28×28 grayscale images as input.

lenet=NetChain[
{ConvolutionLayer[20,5],Ramp,PoolingLayer[2,2],ConvolutionLayer[50,5],Ramp,PoolingLayer[2,2],FlattenLayer[],500,Ramp,10,SoftmaxLayer[]},
"Output"->NetDecoder[{"Class",Range[0,9]}],
"Input"->NetEncoder[{"Image",{28,28},"Grayscale"}]
]

Train the network for four training rounds.

lenet=NetTrain[lenet,trainingData,ValidationSet->testData,MaxTrainingRounds-> 3];

Evaluate the trained network directly on images randomly sampled from the validation set.

imgs=Keys @ RandomSample[testData,8];
Thread[imgs->lenet[imgs]]

[TALK] Neural Networks — Friday 10-10:30pm

[TALK] Deep Learning Applied to Computer Vision — Friday 3:30-4pm

New in 11 – Marketing Pages

Advances in Image & Signal Processing

Computational Audio

Computational Photography

Neural Networks

Time Series Processing

Thanks!