Image Processing

These notes follow the steps of Mathematica's tutorial included in the image processing package, although additional or different processes are used for some tasks as edge finding and color segmentation.

This loads the image processing package.

<<ImageProcessing`

Read an image (Kevin Raskoff's "Orchistoma") and display it.

OrchistomaColor = ImageRead["/home/salvador/images/Raskoff_orchistoma.jpg"] ;

Show[Graphics[OrchistomaColor]] ;

[Graphics:HTMLFiles/imageprocessing_4.gif]

Many image processing tasks are performed on gray level images, so we obtain a gray level version of the original.

Orchistoma = ToGrayLevel[OrchistomaColor] ;

Show[Graphics[Orchistoma]]

[Graphics:HTMLFiles/imageprocessing_7.gif]

-Graphics -

Edges are readily found with the Sobel filter. The first matrix is a derivative operator on the horizontal direction and the second on the vertical direction. Both are convolved with the image when obtaining the gradient magnitude and the orientation map.

MatrixForm/@SobelFilter[]

{( {{-0.25, 0., 0.25}, {-0.5, 0., 0.5}, {-0.25, 0., 0.25}} ), ( {{-0.25, -0.5, -0.25}, {0., 0., 0.}, {0.25, 0.5, 0.25}} )}

{OrchistomaMag, OrchistomaDir} = {EdgeMagnitude[Orchistoma, SobelFilter[]], EdgeOrientation[Orchistoma, SobelFilter[]]}

{-ImageData-, -ImageData-}

Show[GraphicsArray[Graphics/@{OrchistomaMag, OrchistomaDir}]] ;

[Graphics:HTMLFiles/imageprocessing_14.gif]

The edges can be obtained by thresholding the gradient magnitude. To set a reasonable threshold we sometimes need to inspect the histogram.

OrchistomaHis = CumulativeHistogram[ImageHistogram[Round[OrchistomaMag]]] ;

ListPlot[OrchistomaHis, PlotJoined→True, PlotRange→All] ;

[Graphics:HTMLFiles/imageprocessing_17.gif]

Length[Select[OrchistomaHis, (#1[[2]] <0.87&)]]

21

Show[GraphicsArray[(Graphics[Threshold[OrchistomaMag, #1]] &)/@{% + 7, %, % - 7}]] ;

[Graphics:HTMLFiles/imageprocessing_21.gif]

MatrixForm /@ {LaplacianFilter[]}

{( {{0.666667, -0.333333, 0.666667}, {-0.333333, -1.33333, -0.333333}, {0.666667, -0.333333, 0.666667}} )}

OrchistomaLap = EdgeMagnitude[Orchistoma, LaplacianFilter[]] ;

Show[GraphicsArray[Graphics /@ {Orchistoma, OrchistomaLap}]] ;

[Graphics:HTMLFiles/imageprocessing_26.gif]

Show[GraphicsArray[Graphics/@{OrchistomaLap, ZeroCrossing[OrchistomaLap]}]] ;

[Graphics:HTMLFiles/imageprocessing_28.gif]

(g = Exp[-(x^2 + y^2)/(2σ^2)] ; Δg = ∂_ {x, 2} g + ∂_ {y, 2} g // Simplify)

(^(-(x^2 + y^2)/(2 σ^2)) (x^2 + y^2 - 2 σ^2))/σ^4

[Graphics:HTMLFiles/imageprocessing_32.gif]

ListDensityPlot[LoGFilter[40, 40]] ;

[Graphics:HTMLFiles/imageprocessing_34.gif]

EdgeMagnitude[Orchistoma, LoGFilter[40, 40]] ;

Show[GraphicsArray[Graphics/@{%, ZeroCrossing[%]}]] ;

[Graphics:HTMLFiles/imageprocessing_37.gif]

Show[GraphicsArray[Graphics/@{Orchistoma, DiscreteConvolve[Orchistoma, UnsharpFilter[5, 5]]}]] ;

[Graphics:HTMLFiles/imageprocessing_39.gif]

EdgeMagnitude[Orchistoma, UnsharpFilter[5, 5]] ;

Show[GraphicsArray[Graphics/@{%, ZeroCrossing[%]}]] ;

[Graphics:HTMLFiles/imageprocessing_42.gif]

EdgeMagnitude[Orchistoma, SmoothingFilter[7, 5, 7, 5]] ;

Show[GraphicsArray[Graphics/@{%, ZeroCrossing[%]}]] ;

[Graphics:HTMLFiles/imageprocessing_45.gif]

[Graphics:HTMLFiles/imageprocessing_47.gif]

ListDensityPlot[Reverse[Table[Evaluate[g/.{σ→1}], {y, -3, 3, 0.2}, {x, -3, 3, 0.2}]]] ;

[Graphics:HTMLFiles/imageprocessing_49.gif]

gf = Reverse[Table[Evaluate[g/.{σ→1}], {y, -3, 3, 0.2}, {x, -3, 3, 0.2}]] ;

Show[GraphicsArray[Graphics/@{Orchistoma, DiscreteConvolve[Orchistoma, gf, Centered→True]}]] ;

[Graphics:HTMLFiles/imageprocessing_52.gif]

gf1 = Reverse[Table[Evaluate[g/.{σ→0.2}], {y, -3, 3, 0.2}, {x, -3, 3, 0.2}]] ;

ListDensityPlot[gf1] ;

[Graphics:HTMLFiles/imageprocessing_55.gif]

Show[GraphicsArray[Graphics/@{Orchistoma, DiscreteConvolve[Orchistoma, gf1, Centered→True]}]] ;

[Graphics:HTMLFiles/imageprocessing_57.gif]

derGx = ∂_xg//Simplify

-(^(-(x^2 + y^2)/(2 σ^2)) x)/σ^2

[Graphics:HTMLFiles/imageprocessing_61.gif]

fDerG1x = Reverse[Table[Evaluate[derGx/.{σ→0.2}], {y, -3, 3, 0.2}, {x, -3, 3, 0.2}]] ;

ListDensityPlot[fDerG1x] ;

[Graphics:HTMLFiles/imageprocessing_64.gif]

derGy = ∂_yg//Simplify

General :: spell1 : Possible spelling error: new symbol name \"derGy\" is similar to existing symbol \"derGx\".  More…

-(^(-(x^2 + y^2)/(2 σ^2)) y)/σ^2

fDerG1y = Reverse[Table[Evaluate[derGy/.{σ→0.2}], {y, -3, 3, 0.2}, {x, -3, 3, 0.2}]] ;

General :: spell1 : Possible spelling error: new symbol name \"fDerG1y\" is similar to existing symbol \"fDerG1x\".  More…

ListDensityPlot[fDerG1y] ;

[Graphics:HTMLFiles/imageprocessing_71.gif]

Show[GraphicsArray[Graphics/@{Orchistoma, DiscreteConvolve[Orchistoma, fDerG1x, Centered→True]}]] ;

[Graphics:HTMLFiles/imageprocessing_73.gif]

Show[GraphicsArray[Graphics/@{Orchistoma, DiscreteConvolve[Orchistoma, fDerG1y, Centered→True]}]] ;

[Graphics:HTMLFiles/imageprocessing_75.gif]

gradientOrchistoma = Sqrt[DiscreteConvolve[Orchistoma, fDerG1x, Centered→True]^2 + DiscreteConvolve[Orchistoma, fDerG1y, Centered→True]^2] ;

Show[GraphicsArray[Graphics/@{Orchistoma, gradientOrchistoma}]] ;

[Graphics:HTMLFiles/imageprocessing_78.gif]

Max[gradientOrchistoma]

2300.46

Min[gradientOrchistoma]

0.

The edges can also be found by applying a threshold on the gradient magnitude image.

Show[Graphics[Threshold[gradientOrchistoma, 210]]] ;

[Graphics:HTMLFiles/imageprocessing_84.gif]

Another method would be to rescale the gradient magnitude image and compute the zero crossings.

Show[Graphics[ZeroCrossing[ScaleLinear[gradientOrchistoma, {-1, 10}]]]] ;

[Graphics:HTMLFiles/imageprocessing_86.gif]

RegionProcessing[0&, OrchistomaColor, Where[Threshold[gradientOrchistoma, 210], _ ? (# == 1&)]] ;

These are the homogeneous regions where the magnitude of the gradient is low.

Show[GraphicsArray[Graphics/@{OrchistomaColor, %}]] ;

[Graphics:HTMLFiles/imageprocessing_89.gif]

These areas are where most variation in intensity is found. We see that all feature finding based on derivatives will almost always lead to very sparse sets.

RegionProcessing[0&, OrchistomaColor, Where[Threshold[gradientOrchistoma, 210], _ ? (#≠1&)]] ;

Show[GraphicsArray[Graphics/@{OrchistomaColor, %}]] ;

[Graphics:HTMLFiles/imageprocessing_92.gif]

We normalize the filters to range from -1 to 1.

{NfDerG1x, NfDerG1y} = (2 * Rescale[#, {Min[#], Max[#]}] - 1) &/@{fDerG1x, fDerG1y} ;

General :: spell1 : Possible spelling error: new symbol name \"NfDerG1x\" is similar to existing symbol \"fDerG1x\".  More…

General :: spell : Possible spelling error: new symbol name \"NfDerG1y\" is similar to existing symbols  {fDerG1y, NfDerG1x} .  More…

({Min[#], Max[#]}) &/@{NfDerG1x, NfDerG1y}

{{-1., 1.}, {-1., 1.}}

The gradient magnitude and the orientation of the edges are found with Mathematica's primitives contained in the Image Processing package.

{OrchistomaMagG, OrchistomaDirG} = {EdgeMagnitude[Orchistoma, {NfDerG1x, NfDerG1y}], EdgeOrientation[Orchistoma, {NfDerG1x, NfDerG1y}]} ;

General :: spell1 : Possible spelling error: new symbol name \"OrchistomaMagG\" is similar to existing symbol \"OrchistomaMag\".  More…

General :: spell1 : Possible spelling error: new symbol name \"OrchistomaDirG\" is similar to existing symbol \"OrchistomaDir\".  More…

({Min[#], Max[#]}) &/@{OrchistomaMagG, OrchistomaDirG}

{{0., 758.563}, {-180., 180.}}

Show[GraphicsArray[Graphics/@{OrchistomaMagG, OrchistomaDirG}]] ;

[Graphics:HTMLFiles/imageprocessing_104.gif]

Let's compare the gradient magnitude maps obtained from the Sobel filters and the normalized Gaussian derivative filters.

Show[GraphicsArray[Graphics/@{OrchistomaMag, OrchistomaMagG}]] ;

[Graphics:HTMLFiles/imageprocessing_106.gif]

Let's compare the orientation maps obtained from the Sobel filters and the normalized Gaussian derivative filters.

Show[GraphicsArray[Graphics/@{OrchistomaDir, OrchistomaDirG}]] ;

[Graphics:HTMLFiles/imageprocessing_108.gif]

Now let's take a square subimage from Orchistoma.

ImageDimensions[Orchistoma]

{350, 255}

sqrOrchistoma = ImageTake[Orchistoma, {96, 350}, {1, 255}] ;

ImageDimensions[sqrOrchistoma]

{255, 255}

Show[Graphics[sqrOrchistoma]] ;

[Graphics:HTMLFiles/imageprocessing_115.gif]

Take its Fourier transform and display it.

dftOrchistoma = DFT[sqrOrchistoma] ;

[Graphics:HTMLFiles/imageprocessing_118.gif]

Reconstructing the image with the inverse transform we notice almost no difference.

Show[GraphicsArray[Graphics/@{sqrOrchistoma, ListDensityPlot[IDFT[dftOrchistoma]//Chop, Mesh→False, FrameTicks→None, DisplayFunction→Identity]}]] ;

[Graphics:HTMLFiles/imageprocessing_120.gif]

Reconstructing solely from the magnitude does not give a very recognizable pattern,

Show[GraphicsArray[Graphics/@{sqrOrchistoma, ListDensityPlot[IDFT[Abs[dftOrchistoma]]//Re, Mesh→False, FrameTicks→None, DisplayFunction→Identity]}]] ;

[Graphics:HTMLFiles/imageprocessing_122.gif]

whereas reconstructing solely from the phase information does give a recognizable result, we can see two symmetrical versions of the orchistoma, showing the well known fact that our visual system is more sensitive to phase cues.

Show[GraphicsArray[Graphics/@{sqrOrchistoma, ListDensityPlot[IDFT[Arg[dftOrchistoma]]//Im, Mesh→False, FrameTicks→None, DisplayFunction→Identity]}]] ;

[Graphics:HTMLFiles/imageprocessing_124.gif]

Let's try some color processing.

Show[Graphics[OrchistomaColor]] ;

[Graphics:HTMLFiles/imageprocessing_126.gif]

Show[Graphics[PlanarImageData[OrchistomaColor]]] ;

[Graphics:HTMLFiles/imageprocessing_128.gif]

We'll try to isolate the yellow regions. First we need to threshold each color component

Clear[orchHis] ;

orchHis = ImageHistogram[OrchistomaColor] ;

ShowImageHistogram[orchHis, Ticks→ {{{20, 60, 100, 200, 255}, Automatic}, {{20, 60, 100, 200, 255}, Automatic}, {{20, 60, 100, 200, 255}, Automatic}}] ;

[Graphics:HTMLFiles/imageprocessing_132.gif]

After a while we realize that the histograms were completely useless in this case. The trick is make your thresholded images look like the original's color planes, at least in the region of interest.

thr = Threshold[OrchistomaColor, {90, 150, 190}] ;

Show[Graphics[PlanarImageData[thr]]] ;

[Graphics:HTMLFiles/imageprocessing_135.gif]

To isolate the yellow region we put zeros on all pixels which are not red and green (yellow).

RegionProcessing[0&, OrchistomaColor, Where[thr, _ ? (#≠ {1, 1, 0} &)]] ;

Show[GraphicsArray[Graphics/@{OrchistomaColor, %}]] ;

[Graphics:HTMLFiles/imageprocessing_138.gif]

It seems that if we're segmenting by color then the HSV color space would be the one to use.

OrchistomaHSV = ToHSVColor[OrchistomaColor] ;

Show[Graphics[PlanarImageData[OrchistomaHSV]]] ;

[Graphics:HTMLFiles/imageprocessing_141.gif]

Yellow's hue numerical value should be exactly between red (0) and green (1/3), that is, it should be 0.16666666… Here we select "yellowish" hues in the vicinity of 1/6. We then get rid of the bluish streaks of the previous segmentation.

RegionProcessing[0&, OrchistomaColor, Where[ToChannels[OrchistomaHSV][[1]], _ ? (#>0.18 || #<0.04 &)]] ;

Show[GraphicsArray[Graphics/@{OrchistomaColor, %}]] ;

[Graphics:HTMLFiles/imageprocessing_146.gif]

We would like to take the white tentacles away, so we further eliminate low saturations. The isolated region will be kept in the variable "mask".

mask = RegionProcessing[0&, %%, Where[ToChannels[OrchistomaHSV][[2]], _ ? (#<0.22&)]] ;

Show[GraphicsArray[Graphics/@{OrchistomaColor, %}]] ;

[Graphics:HTMLFiles/imageprocessing_149.gif]

We would like to change the color of the yellow part of the original picture to cyan (whose hue is 0.5), but before we must explore how to do it because the documentation is not all that clear.

First, let's find a little portion lf the original,

Show[Graphics[ImageTake[OrchistomaHSV, {150, 175}, {75, 100}]]]

[Graphics:HTMLFiles/imageprocessing_151.gif]

-Graphics -

then we display its components (HSV).

Show[Graphics[PlanarImageData[ImageTake[OrchistomaHSV, {150, 175}, {75, 100}]]]]

[Graphics:HTMLFiles/imageprocessing_154.gif]

-GraphicsArray -

It seems that the results of the previous region processing stored in "mask" are given in RGB regardless of the format of the image being processed, so we need to convert it to HSV.

maskHSV = ToHSVColor[mask] ;

Now let's see how can we access the hue plane.

RawImageData[ImageTake[maskHSV, {150, 175}, {75, 100}]][[All, All, 1]]//MatrixForm

This confirms that ToHSVColor can be an 'inverse' of RawImageData.

Show[Graphics[ToHSVColor[RawImageData[maskHSV]]]]

[Graphics:HTMLFiles/imageprocessing_160.gif]

-Graphics -

This creates an array of raw data that would be identical to the hue plane save that all hues in the isolated region have been turned to a reddish hue.

res = RegionProcessing[0&, RawImageData[OrchistomaHSV][[All, All, 1]], Where[RawImageData[maskHSV][[All, All, 3]], _ ? (#≠0&)]] ;

Plotting a small part of the result shows that we got what we expected.

ListDensityPlot[Take[res, {150, 175}, {75, 100}]]

[Graphics:HTMLFiles/imageprocessing_164.gif]

-DensityGraphics -

If we want to see the whole thing we need to turn the raw data into a gray level image.

Show[Graphics[ToGrayLevel[res]]] ;

[Graphics:HTMLFiles/imageprocessing_167.gif]

We will create a new image in the HSV format having the hue plane substituted by the results of our previous region processing, and the saturation and value planes come from the original image. To do this, we create the corresponding gray level images.

res1 = ToGrayLevel[res] ;

res2 = ToGrayLevel[ToChannels[OrchistomaHSV][[2, 1]]] ;

res3 = ToGrayLevel[ToChannels[OrchistomaHSV][[3, 1]]] ;

Show[GraphicsArray[Graphics/@{res1, res2, res3}]] ;

[Graphics:HTMLFiles/imageprocessing_172.gif]

Now we combine the previous images into a single HSV image and see the result.

result = ToHSVColor[RawImageData[FromChannels[res1, res2, res3]]] ;

Show[Graphics[result]] ;

[Graphics:HTMLFiles/imageprocessing_175.gif]

Before and after.

Show[GraphicsArray[Graphics/@{OrchistomaHSV, result}]] ;

[Graphics:HTMLFiles/imageprocessing_177.gif]


Created by Mathematica  (June 15, 2006) Valid XHTML 1.1!