Author Archives: Andrew Collier

#MonthOfJulia Day 38: Imaging

Julia has a few packages aimed at image processing. We’ll start by looking at the TestImages package, which hosts a selection of sample images, then briefly visit the ImageView package before moving onto the Images package, which implements a range of functions for image manipulation.

Test Images

The TestImages package currently provides 25 sample images, which form a convenient basis for experimentation.

julia> using TestImages
julia> readdir(joinpath(homedir(), ".julia/v0.4/TestImages/images/"))
25-element Array{ByteString,1}:
 "autumn_leaves.png" 
 "cameraman.tif"     
 "earth_apollo17.jpg"
 "fabio.png"         
 "house.tif"         
 "jetplane.tif"      
 "lake.tif"          
 "lena_color_256.tif"
 "lena_color_512.tif"
 "lena_gray_256.tif" 
 "lena_gray_512.tif" 
 "lighthouse.png"    
 "livingroom.tif"    
 "mandril_color.tif" 
 "mandril_gray.tif"  
 "mandrill.tiff"     
 "moonsurface.tiff"  
 "mountainstream.png"
 "peppers_color.tif" 
 "peppers_gray.tif"  
 "pirate.tif"        
 "toucan.png"        
 "walkbridge.tif"    
 "woman_blonde.tif"  
 "woman_darkhair.tif"

We’ll load the archetypal test image (the November 1972 Playboy centerfold of Lena Söderberg).

julia> lena = testimage("lena_color_256.tif");

Of course, now that we’ve loaded that image, we’ll want to take a look at it. To do that we’ll need the ImageView package.

julia> using ImageView
julia> view(lena)
(ImageCanvas,ImageSlice2d: zoom = Graphics.BoundingBox(0.0,256.0,0.0,256.0))

You can optionally specify the pixel spacing as a parameter to view(), which then ensures that the aspect ratio of the image is conserved on resizing. There are various other bobs and whistles associated with view(): you can click-and-drag within the image to zoom in on a particular region; various simple transformations (flipping and rotation) are possible; images can be annotated and multiple images can be arranged on a canvas for simultaneous viewing.

Image Representation

Outside of the test images, an arbitrary image file can be loaded using imread() from the Images package. Naturally, there are also functions for writing images, imwrite() and writemime().

julia> using Images
julia> earth = imread(joinpath(homedir(), ".julia/v0.4/TestImages/images/earth_apollo17.jpg"))
RGB Images.Image with:
  data: 3000x3002 Array{ColorTypes.RGB{FixedPointNumbers.UfixedBase{UInt8,8}},2}
  properties:
    IMcs: sRGB
    spatialorder:  x y
    pixelspacing:  1 1

The default representation for the Image object tells us its dimensions, storage type and colour space. The spatial order indicates that the image data are stored using row major ordering. It’s also possible to specify physical units for the pixel spacing, which is particularly important if you are analysing images where absolute scale matters (for example, medical imaging). There are convenience methods for a few image properties.

julia> colorspace(earth)
"RGB"
julia> height(earth)
3002
julia> width(earth)
3000

We can examine individual pixels within the image using the indexing operator.

julia> earth[1,1]
RGB{U8}(0.047,0.008,0.0)

Each pixel is of type RGB (defined in the Colors package), which encapsulates a tuple giving the proportion of red, green and blue for that pixel. The underlying image data can also be accessed via the data() method.

The image can be split into its component colour channels using separate().

julia> earth_rgb = separate(earth)
RGB Images.Image with:
  data: 3002x3000x3 Array{FixedPointNumbers.UfixedBase{UInt8,8},3}
  properties:
    IMcs: sRGB
    colorspace: RGB
    colordim: 3
    spatialorder:  y x
    pixelspacing:  1 1

Note that the result is a three-dimensional Array. The spatial order has also changed, which means that the data are now represented using column major ordering. The data are thus effectively transposed.

Simple Image Processing

Kernel-based filtering can be applied using imfilter() or imfilter_fft(), where the latter is better suited to larger kernels. There’s a variety of helper functions for constructing kernels, like imaverage() and gaussian2d().

julia> lena_smooth = imfilter(lena, imaverage([3, 3]));
julia> lena_very_smooth = imfilter_fft(lena, ones(10, 10) / 100);
julia> lena_gauss_smooth = imfilter_gaussian(lena, [1, 2]);

The effects of the above smoothing operations can be seen below, with the original image on the left, followed by the 3-by-3 and 10-by-10 boxcar filtered versions and finally the Gaussian filtered image.

The imgradients() function calculates gradients across the image. You can choose from a set of methods for calculating the gradient. The morphological dilation and erosion operations are available via dilate() and erode().

julia> (lena_sobel_x, lena_sobel_y) = imgradients(lena, "sobel");
julia> lena_dilate = dilate(lena);
julia> lena_erode = erode(lena);

Below are the two components of the image gradient calculated using the Sobel operator followed by the results of dilate() and erode().

Other Packages

The ImageMagick package implements further imaging functionality. If, in the future, it provides an interface to the full functionality on the ImageMagick suite then it will be a truly phenomenal resource. Also worth looking at is the PiecewiseAffineTransforms package which implements a technique for warping portions of an image.

If you’d like to exercise your image processing and machine learning skills in Julia, take a look at the First Steps With Julia competition on kaggle.

The post #MonthOfJulia Day 38: Imaging appeared first on Exegetic Analytics.

#MonthOfJulia Day 37: Fourier Techniques

The Fourier Transform is often applied to signal processing and other analyses. It allows a signal to be transformed between the time domain and the frequency domain. The efficient Fast Fourier Transform (FFT) algorithm is implemented in Julia using the FFTW library.

1D Fourier Transform

Let’s start by looking at the Fourier Transform in one dimension. We’ll create test data in the time domain using a wide rectangle function.

julia> f = [abs(x) <= 1 ? 1 : 0 for x in -5:0.1:5];
julia> length(f)
101

This is what the data look like:

We’ll transform the data into the frequency domain using fft().

julia> F = fft(f);
julia> typeof(F)
Array{Complex{Float64},1}
julia> length(F)
101
julia> F = fftshift(F);

The frequency domain data are an array of Complex type with the same length as the time domain data. Since each Complex number consists of two parts (real and imaginary) it seems that we have somehow doubled the information content of our signal. This is not true because half of the frequency domain data are redundant. The fftshift() function conveniently rearranges the data in the frequency domain so that the negative frequencies are on the left.

This is what the resulting amplitude and power spectra look like:

The analytical Fourier Transform of the rectangle function is the sinc function, which agrees well with numerical data in the plots above.

2D Fourier Transform

Let’s make things a bit more interesting: we’ll look at the analogous two-dimensional problem. But this time we’ll go in the opposite direction, starting with a two-dimensional sinc function and taking its Fourier Transform.

Building the array of sinc data is easy using a list comprehension.

julia> f = [(r = sqrt(x^2 + y^2); sinc(r)) for x in -6:0.125:6, y in -6:0.125:6];
julia> typeof(f)
Array{Float64,2}
julia> size(f)
(97,97)

It doesn’t make sense to think about a two-dimensional function in the time domain. But the Fourier Transform is quite egalitarian: it’s happy to work with a temporal signal or a spatial signal (or a signal in pretty much any other domain). So let’s suppose that our two-dimensional data are in the spatial domain. This is what it looks like:

Generating the Fourier Transform is again a simple matter of applying fft(). No change in syntax: very nice indeed!

julia> F = fft(f);
julia> typeof(F)
Array{Complex{Float64},2}
julia> F = fftshift(F);

The power spectrum demonstrates that the result is the 2D analogue of the rectangle function.

Higher Dimensions and Beyond

It’s just as easy to apply the FFT to higher dimensional data, although in my experience this is rarely required.

Most of the FFTW library’s functionality has been implemented in the Julia interface. For example:

it’s possible to generate plans for optimised FFTs using plan_fft();
dct() yields the Discrete Cosine Transform;
you can exploit conjugate symmetry in real transforms using rfft(); and
it’s possible to run over multiple threads using FFTW.set_num_threads().

Watch the video below in which Steve Johnson demonstrates many of the features of FFTs in Julia.

The post #MonthOfJulia Day 37: Fourier Techniques appeared first on Exegetic Analytics.

#MonthOfJulia Day 36: Markdown

Markdown is a lightweight format specification language developed by John Gruber. Markdown can be converted to HTML, LaTeX or other document formats. You probably knew all that already. The syntax is pretty simple. Check out this useful cheatsheet.

In the latest stable version of Julia support for markdown is provided in the Base package.

julia> using Base.Markdown
julia> import Base.Markdown: MD, Paragraph, Header, Italic, Bold, LineBreak, plain, term, html,
                             Table, Code, LaTeX, writemime

Markdown is stored in objects of type Base.Markdown.MD. As you’ll see below, there are at least two ways to construct markdown objects: either directly from a string (using markdown syntax) or programmatically (using a selection of formatting functions).

julia> d1 = md"foo *italic foo* **bold foo** `code foo`";
julia> d2 = MD(Paragraph(["foo ", Italic("italic foo"), " ", Bold("bold foo"), " ",
               Code("code foo")]));
julia> typeof(d1)
Base.Markdown.MD
julia> d1 == d2
true

You’ll find that Base.Markdown.MD objects are rendered with appropriate formatting in your console.

Functions html() and latex() convert Base.Markdown.MD objects into other formats. Another way of rendering markdown elements is with writemime(), where the output is determined by specifying a MIME type.

julia> html(d1)
"<p>foo <em>italic foo</em> <strong>bold foo</strong> <code>code foo</code></p>n"
julia> latex(d1)
"foo \emph{italic foo} \textbf{bold foo} \texttt{code foo}n"

Markdown has support for section headers, both ordered and unordered lists, tables, code fragments and block quotes.

julia> d3 = md"""# Chapter Title
       ## Section Title
       ### Subsection Title""";
julia> d4 = MD(Header{2}("Section Title"));
julia> d3 |> html
"<h1>Chapter Title</h1>n<h2>Section Title</h2>n<h3>Subsection Title</h3>n"
julia> latex(d4)
"\subsection{Section Title}n"

Most Julia packages come with a README.md markdown file which provides an overview of the package. The readme() function gives you direct access to these files’ contents.

julia> readme("Quandl")
  Quandl.jl
  ============
  (Image: Build Status)
  (Image: Coverage Status)
  (Image: Quandl)

  Documentation is provided by Read the Docs.

  See the Quandl API Help Page for further details about the Quandl API. This package 
  closely follows the nomenclature used by that documentation.

We can also use parse_file() to treat the contents of a file as markdown.

julia> d6 = Markdown.parse_file(joinpath(homedir(), ".julia/v0.4/NaNMath/README.md"));

This is rendered below as LaTeX.

section{NaNMath}
Implementations of basic math functions which return texttt{NaN} instead of throwing a
texttt{DomainError}.
Example:
begin{verbatim}
import NaNMath
NaNMath.log(-100) # NaN
NaNMath.pow(-1.5,2.3) # NaN
end{verbatim}
In addition this package provides functions that aggregate one dimensional arrays and ignore
elements that are NaN. The following functions are implemented:
begin{verbatim}
sum
maximum
minimum
mean
var
std
end{verbatim}
Example:
begin{verbatim}
using NaNMath; nm=NaNMath
nm.sum([1., 2., NaN]) # result: 3.0
end{verbatim}
href{https://travis-ci.org/mlubin/NaNMath.jl}{begin{figure}
centering
includegraphics{https://travis-ci.org/mlubin/NaNMath.jl.svg?branch=master}
caption{Build Status}
end{figure}
}

And here it is as HTML.

NaNMath

Implementations of basic math functions which return NaN instead of throwing a DomainError.
Example:

import NaNMath
NaNMath.log(-100) # NaN
NaNMath.pow(-1.5,2.3) # NaN

In addition this package provides functions that aggregate one dimensional arrays and ignore elements that are NaN. The following functions are implemented:

sum
maximum
minimum
mean
var
std

Example:

using NaNMath; nm=NaNMath
nm.sum([1., 2., NaN]) # result: 3.0

What particularly appeals to me about the markdown functionality in Julia is the potential for automated generation of documentation and reports. To see more details of my dalliance with Julia and markdown, visit github.

The post #MonthOfJulia Day 36: Markdown appeared first on Exegetic Analytics.

juliabloggers.com

A Julia Language Blog Aggregator