The first MathWorks general product release of the year, R2013a, shipped a couple of months ago. I've already mentioned it once here in my 12-Mar-2013 post about the new MATLAB unit test framework.

With each new release, I peruse the release notes for MATLAB to see what things I find particularly interesting. (This helps me remember which product features have actually been released, as opposed to still being in development. My memory needs all the help it can get.)

The first thing to note is the reappearance of the table of contents for navigating in the Help Browser and in the online Documentation Center. This is a direct result of helpful feedback we received from many of you about the R2012b release.

My favorite "make-it-go-faster-without-sacrificing-accuracy" people (the MATLAB Math Team, that is) have been busy again. People with computers based on Intel or AMD chips using the AVX instruction set should see their calls to fft speed up. Anybody running permute on 3-D or higher-dimensional arrays should also get a nice boost. I've done a lot of development work related to image and scientific format support, so I know that a fast permute can be pretty useful when reading image and scientific data. That's because most of these formats store array elements in the file in a different order than MATLAB uses in memory.

In the small-but-nice category, the MATLAB Math Team also simplified a common programming pattern in my own neck of the woods (image processing). Specifically, it's a bit easier to initial an array of 0s or 1s whose type is based on existing array. Here's an example to illustrate:

clear  % Let's start with a fresh workspace.
rgb = imread('peppers.png');
imshow(rgb)
title('Obligatory image screenshot')

Now I want a 100-by-100 matrix of 0s with the same data type as rgb.

A = zeros(100,100,'like',rgb); % Make a 100-by-100 matrix that's "like" rgb.
whos

  Name        Size                Bytes  Class    Attributes

  A         100x100               10000  uint8              
  rgb       384x512x3            589824  uint8              

My developer friend Tom Bryan really "likes" this (ahem) because it enables much easier solutions to some common programming tasks for users of Fixed-Point Designer.

I have occasionally done a little web scripting in MATLAB, so it's nice to see urlread and urlwrite get a little love. These functions can now handle basic authentication via the 'Authentication', 'Username', and 'Password' parameters.

Do you use a Mac? You can now write MPEG-4 H.264 files using VideoWriter (requires Mac OS 10.7 or later).

A couple of handy new string functions have appeared, strsplit and strjoin. Based on how often users have submitted their own versions to the MATLAB Central File Exchange, I'm sure these will be popular.

out = strsplit(pwd,'\')

out = 

    'B:'    'published'    '2013'

You can now do extrapolation with both scattered and gridded interpolation. For extrapolation with scattered interpolation, use the new scatteredInterpolant. Here's an example I lifted from the doc.

Query the interpolant at a single point outside the convex hull using nearest neighbor extrapolation.

Define a matrix of 200 random points.

P = -2.5 + 5*gallery('uniformdata',[200 2],0);

Sample an exponential function. These are the sample values for the interpolant.

x = P(:,1);
y = P(:,2);
v = x.*exp(-x.^2-y.^2);

Create the interpolant, specifying linear interpolation and nearest neighbor extrapolation.

F = scatteredInterpolant(P,v,'linear','nearest')

F = 

  scatteredInterpolant with properties:

                 Points: [200x2 double]
                 Values: [200x1 double]
                 Method: 'linear'
    ExtrapolationMethod: 'nearest'

Evaluate the interpolant outside the convex hull.

vq = F(3.0,-1.5)

vq =

    0.0031

Disable extrapolation and evaluate F at the same point.

F.ExtrapolationMethod = 'none';
vq = F(3.0,-1.5)

vq =

   NaN

I encourage you to wander over to the R2013a release notes for MATLAB or any other product that you use and see what's new that might be helpful to you.

There are also lots of new things in the image processing and computer vision worlds, of course. I'll look at those next time.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

I was reviewing enhancement requests recently when I came across this one: Add support for a 'Quality' parameter when using imwrite to write a PNG file, just like you can currently do when writing a JPEG file.

Well, there is actually a pretty good reason why there is no 'Quality' parameter for writing PNG files, but it's not an obvious reason unless you know more about the differences between various image compression methods.

A very important characteristic of a compression method is whether it is lossless or lossy. With a lossless compression method, the original data can be recovered exactly. When you make a "zip" file on your computer, this is what you certainly expect.

Here's an example: I can use gzip to compress the MATLAB script file I'm using for this blog post.

gzip('lossless_lossy.m');
dir('lossless_lossy.*')

lossless_lossy.m     lossless_lossy.m.gz  

Let's compare their sizes to see if the output of gzip is really compressed.

d1 = dir('lossless_lossy.m');
d1.bytes

ans =

        4098

d2 = dir('lossless_lossy.m.gz');
d2.bytes

ans =

        1785

Indeed, the compressed does actually have fewer bytes. Can we get the original file back exactly?

gunzip('lossless_lossy.m.gz','./tmp')
isequal(fileread('lossless_lossy.m'),fileread('./tmp/lossless_lossy.m'))

ans =

     1

Yes.

Let's switch to image formats. The PNG image format uses lossless compression. When you save image data to a PNG file, you can read the file back in and get back the original pixels, unchanged. For a sample image I'll use my imzoneplate function on the MATLAB Central File Exchange.

I = im2uint8(imzoneplate);
imshow(I)

Let's write I out to a PNG file, read it back in, and see if the pixels are the same.

imwrite(I,'zoneplate.png');
I2 = imread('zoneplate.png');
isequal(I,I2)

ans =

     1

OK, now let's try the same experiment using JPEG.

imwrite(I,'zoneplate.jpg')
I2j = imread('zoneplate.jpg');
isequal(I,I2j)

ans =

     0

No, the pixels are not equal! It turns out the JPEG is a lossy image compression format. (Full disclosure: there is a lossless variant of JPEG, but it is rarely used.)

Why in the world would we use a compression format that doesn't preserve the original data? Because by giving up on exact data recovery and by taking advantage of properties of human visual perception, we can make the stored file a lot smaller. Let's compare the file sizes of the PNG file with the JPEG file.

z1 = dir('zoneplate.png');
num_bytes_png = z1.bytes

num_bytes_png =

      218864

z2 = dir('zoneplate.jpg');
num_bytes_jpeg = z2.bytes

num_bytes_jpeg =

       72660

The JPEG file is only one-third the size of the PNG file! But it looks almost exactly the same.

imshow('zoneplate.jpg')

So, what about that 'Quality' parameter that I mentioned at the top of today's post? It turns out that we can make the JPEG file even smaller if we are willing to put up with some visible distortion in the image. Let's try a quality factor of 25. (The default is 75 on a scale of 0 to 100.)

imwrite(I,'zoneplate_25.jpg','Quality',25)
I2j_25 = imread('zoneplate_25.jpg');
imshow(I2j_25)

You can see some distortion, especially around the high-frequency part of the pattern. How about the file size?

z3 = dir('zoneplate_25.jpg');
z3.bytes

ans =

       39544

That's about 54% of the size of the zoneplate.jpg.

Sometimes people think they can get "lossless JPEG" by using a 'Quality' factor of 100, but unfortunately that isn't the case. Let's check it:

imwrite(I,'zoneplate_100.jpg','Quality',100)
I2j_100 = imread('zoneplate_100.jpg');
isequal(I,I2j_100)

ans =

     0

So there you go -- that's why there's no 'Quality' parameter when writing PNG files. PNG files are always perfect!

When I write blog posts, sometimes I use PNG image files and sometimes I use JPEG. My choice is based on what kind of graphics I have in that particular post. And that's a blog topic for another day.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

This is the fourth and last part of my plan (my evil plan?) to rewrite an Image Processing Toolbox example from 20 years ago using more modern MATLAB language features. I got the idea from Dave Garrison's recent article on writing MATLAB apps.

Here's the old app I'm trying to reinvent:

And here's what I had working last time:

There are two things I want to wrap up before calling it good enough:

First let's get the DCT coefficient mask display working. Recall that last time I added a function to compute the reconstructed image, given the desired number of DCT coefficients. I'll add a second output argument to that function in order to return the DCT coefficient mask. Here's the code. The only changes are on the first line (to define the additional output argument) and the last three lines (to compute the mask).

function [I2,mask] = reconstructImage(I,n)
% Reconstruct the image from n of the DCT coefficients in each 8-by-8
% block. Select the n coefficients with the largest variance across the
% image. Second output argument is the 8-by-8 DCT coefficient mask.

% Compute 8-by-8 block DCTs.
f = @(block) dct2(block.data);
A = blockproc(I,[8 8],f);

% Compute DCT coefficient variances and decide
% which to keep.
B = im2col(A,[8 8],'distinct')';
vars = var(B);
[~,idx] = sort(vars,'descend');
keep = idx(1:n);

% Zero out the DCT coefficients we are not keeping.
B2 = zeros(size(B));
B2(:,keep) = B(:,keep);

% Reconstruct image using 8-by-8 block inverse
% DCTs.
C = col2im(B2',[8 8],size(I),'distinct');
finv = @(block) idct2(block.data);
I2 = blockproc(C,[8 8],finv);

mask = false(8,8);
mask(keep) = true;
end

Next I need some code to visualize the coefficient mask. I want to display it as image with gray lines drawn between the mask pixels. So I added a local function called displayCoefficientMask:

function displayCoefficientMask(mask,ax)
imshow(mask,'Parent',ax)
for k = 0.5:1.0:8.5
    line('XData',[0.5 8.5], ...
        'YData',[k k], ...
        'Color',[0.6 0.6 0.6], ...
        'LineWidth',2, ...
        'Clipping','off', ...
        'Parent',ax);
    line('XData',[k k],...
        'YData',[0.5 8.5],...
        'Color',[0.6 0.6 0.6], ...
        'LineWidth',2,...
        'Clipping','off', ...
        'Parent',ax);
end
title(ax,'DCT Coefficient Mask')
end

The last step is to call displayCoefficientMask from the update method (which gets called whenever the slider is moved). In the code below, I have modified the call to reconstructImage to use two output arguments in order to get the mask; I have assigned the various app properties; and I have added the call to displayCoefficientMask at the end.

function update(app)
% Update the computation
[recon_image,mask] = reconstructImage(app.OriginalImage, ...
    app.NumDCTCoefficients);

diff_image = imabsdiff(app.OriginalImage, recon_image);

% Update the app properties
app.ReconstructedImage = recon_image;
app.ErrorImage = diff_image;
app.DCTCoefficientMask = mask;

% Update the display
imshow(app.OriginalImage,'Parent',app.OriginalImageAxes);
title(app.OriginalImageAxes,'Original Image');

imshow(recon_image,'Parent',app.ReconstructedImageAxes);
title(app.ReconstructedImageAxes,'Reconstructed Image');

imshow(diff_image,[],'Parent',app.ErrorImageAxes);
title(app.ErrorImageAxes,'Error Image');

displayCoefficientMask(mask,app.MaskAxes);

drawnow;
end

Here's the result with the DCT coefficient mask visualization included:

The last thing I want to do with this little app is to allow users to set the app's NumDCTCoefficients property from the command line and to have the app automatically update. To do this, I'll make a couple of changes to the NumDCTCoefficients property. First, I'll make it a dependent property. Instead of being stored independently, this property will be computed on-the-fly from slider setting whenever it is queried. That requires that I define a property get method that computes the property's value on demand. And last I'll need a property set method that defines what actions should be taken whenever the user sets the property.

Here's the modified property block that indicates that NumDCTCoefficients is a dependent property.

properties (Dependent)
    NumDCTCoefficients
end

And here are the get and set functions for NumDCTCoefficients. The get function computes the number on-the-fly based on the current slider position. The set function modifies the slider position and then calls the update method to recompute and redisplay everything.

function value = get.NumDCTCoefficients(app)
    value = round(get(app.Slider,'Value') * 64);
end

function set.NumDCTCoefficients(app,num_coefficients)
    set(app.Slider,'Value',num_coefficients/64)
    update(app);
end

Here's an example of this interaction.

app = dctCompressionExample_v5

app = 

  dctCompressionExample_v5 with properties:

         OriginalImage: [240x240 double]
    ReconstructedImage: [240x240 double]
            ErrorImage: [240x240 double]
    DCTCoefficientMask: [8x8 logical]
    NumDCTCoefficients: 3

I can see the DCT coefficient mask on the app, but I can also look at the DCTCoefficientMask property.

app.DCTCoefficientMask

ans =

     1     1     0     0     0     0     0     0
     1     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0

And I can set the NumDCTCoefficients property, which causes the app to update.

app.NumDCTCoefficients = 1;

OK, I think that's enough to get all the basic ideas. If you want to play around with the final version of the code for this blog post, you can download it from here.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

Although this is the time of year that I normally mention some of the new features from the "a" release, today I thought I would focus on a feature added several releases ago that I think deserves more attention: parallel cluster-based image display. I can't remember exactly when this was released; I think it might have been R2011c.

Anyway, it turns out image display speed can be dramatically increased if we exploit the capabilities of today's massive parallel computing clusters. You start by opening a MATLAB pool from your desktop session using the Parallel Computing Toolbox function matlabpool. The Parallel Computing Toolbox has a new cluster profile that makes it easy to set things up for fast image display:

matlabpool image-display

Once the pool is configured, then the Image Processing Toolbox function imshow automatically uses the workers in the pool.

Let's try it with my blog author picture:

(I'm writing this blog post at home, so for a cluster I'll be using my home-brew video game console, which runs on a cluster of 128 TRS-80s.)

imshow

Here's what the display looks like on worker #37:

And here is the display on worker #92:

If you look carefully, you can see that the pixel block is not exactly the same size on the two workers. That's because of the automatic worker load balancing feature. Way cool!

Let's get down to some numbers. Specifically, I was interested to see how well image display performance scales with the number of workers. I used my function timeit to measure the display speed as I varied the number of workers. The plot below is "normalized speed-up." That is, 1 is the image display speed for normal MATLAB, with no use of the Parallel Computing Toolbox.

As you can see, we get almost perfectly linear scaling with the number of workers. (I'm not sure what's going on with that blip on the left. I think it might have been caused by a bad tape unit on worker #25.)

So, Parallel Computing Toolbox users, please give this a try and let us know about your experience!

Before I wrap it up for today, I can't resist giving a hint about future development directions. (Long-time readers know how much MathWorks loves to preannounce features.) We recently got our hands on one of the experimental new superfluidic clusters. Early indications are that we'll be able to get image display speed to scale with the square of the number of workers! Woo-hoo!

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

This is the third part of my plan to rewrite an Image Processing Toolbox example from 20 years ago using more modern MATLAB language features. I got the idea from Dave Garrison's recent article on writing MATLAB apps.

Here's the old app I'm trying to reinvent:

And here's what I had working the last time:

The code isn't doing any actual DCT computations yet. The idea is to compute 8-by-8 block DCTs, zero out the least significant DCT coefficients of each block, and then reconstruct the image by computing the inverse DCT of each 8-by-8 block.

Let me show how this will work on the sample image I'm using for the app:

pout = imread('pout.tif');
pout2 = pout(1:240,:);
I = im2double(adapthisteq(pout2));
imshow(I)

First, use blockproc to compute the 8-by-8 2-D DCTs.

f = @(block) dct2(block.data);
A = blockproc(I,[8 8],f);
imshow(A,[])
title('DCT coefficients')

Second, I want to compute the variances of each of the 64 coefficients. That is, compute the variance of the (1,1) DCT coefficients from each block, the variance of the (1,2) DCT coefficients from each block, and so on. I'll use the computed variances to determine which DCT coefficients to keep and which to zero out. To do this, I'll rearrange the elements from the previous step so that the corresponding DCT coefficients are together in the columns of a 64-column matrix.

B = im2col(A,[8 8],'distinct')';
vars = var(B);
plot(vars)
title('Variances of the 64 DCT coefficients')

Next, suppose we want to keep the 3 block DCT coefficients with the highest variances. I'll use sort to figure out which ones to keep.

[~,idx] = sort(vars,'descend');
keep = idx(1:3)

keep =

     1     9     2

For an 8-by-8 block of DCT coefficients, index values 1, 9, and 2 correspond to the upper-left DCT coefficient and the two coefficients just the right and down from it. Let's keep just those coefficients.

B2 = zeros(size(B));
B2(:,keep) = B(:,keep);

Now I'll rearrange the columns back into blocks and reconstruct the image by computing the inverse DCT of each block.

C = col2im(B2',[8 8],size(I),'distinct');
finv = @(block) idct2(block.data);
D = blockproc(C,[8 8],finv);
imshow(D)
title('Image reconstructed from 3 DCT coefficients per block')

This algorithm code into a local function to go at the bottom of my app code file, like this:

function I2 = reconstructImage(I,n)
% Reconstruct the image from n of the DCT coefficients in each 8-by-8
% block. Select the n coefficients with the largest variance across the
% image.

% Compute 8-by-8 block DCTs.
f = @(block) dct2(block.data);
A = blockproc(I,[8 8],f);

% Compute DCT coefficient variances and decide
% which to keep.
B = im2col(A,[8 8],'distinct')';
vars = var(B);
[~,idx] = sort(vars,'descend');
keep = idx(1:n);

% Zero out the DCT coefficients we are not keeping.
B2 = zeros(size(B));
B2(:,keep) = B(:,keep);

% Reconstruct image using 8-by-8 block inverse
% DCTs.
C = col2im(B2',[8 8],size(I),'distinct');
finv = @(block) idct2(block.data);
I2 = blockproc(C,[8 8],finv);
end

I already have an update method in my class, so I can add code to that method that will compute the reconstructed image and error image. Here's the revised update method:

function update(app)
    recon_image = reconstructImage(app.OriginalImage, ...
        app.NumDCTCoefficients);

    diff_image = imabsdiff(app.OriginalImage, recon_image);

    imshow(app.OriginalImage,'Parent',app.OriginalImageAxes);
    title(app.OriginalImageAxes,'Original Image');

    imshow(recon_image,'Parent',app.ReconstructedImageAxes);
    title(app.ReconstructedImageAxes,'Reconstructed Image');

    imshow(diff_image,[],'Parent',app.ErrorImageAxes);
    title(app.ErrorImageAxes,'Error Image');

    imshow(zeros(size(app.OriginalImage)),'Parent',app.MaskAxes);
    title(app.MaskAxes,'DCT Coefficient Mask');
    drawnow;
end

Although I don't have all the interactive behavior wired up yet, and don't have the DCT coefficient mask visualization done yet, the app runs and looks a bit more finished than it did before. It now shows the results for reconstructing the image from 3 DCT coefficients per block.

I think it'll take just one more post to finish making the rest of this app work.

Get the MATLAB code

Published with MATLAB® R2013a

Published with MATLAB® R2013a

R2013a, the latest semi-annual MathWorks product release, just went live. One of the significant new capabilities in the MATLAB release is a new unit test framework (overview video, documentation). I am very happy to see this!

Long-time blog readers might be wondering, though, what about MATLAB xUnit? This is a unit test framework that I created and put on the File Exchange. I first wrote about it in this space in 2009. It is my second most popular File Exchange contribution (over 200 downloads in the last 30 days).

Today I want to recap how MATLAB xUnit came to be, explain what will probably happen to it, and convince its users to give the new "official" unit test framework in R2013a a try.

I learned to write unit tests when I came to MathWorks in 1993 as a new software developer. I've used several generations of test frameworks created for internal use here. Until recently, our internal testing machinery was tightly coupled to the entire system we used to build, test, and release our products, and so it was impractical to use it for small projects. This became an issue for me sometime around 2002-2003, as I was working on the first edition of Digital Image Processing Using MATLAB. As the "software guy" among the three coauthors, I was responsible for overseeing the MATLAB functions in the book, including making sure that everything was working OK. So I wrote a simple test harness that exercised the book's examples and some of the functions.

Later, the MATLAB R2008a release included the new generation of object-oriented language features. I was interested in learning more about it. Coincidentally, I had just read Kent Beck's Test-Driven Development, which included a case study on using test-driven development methods to create an xUnit-style test harness. I decided that would be an excellent project with which to learn about test-driven development as well as the R2008a MATLAB language changes. (With the wisdom of hindsight, I have to say here that I don't actually recommend learning test-driven development by using test-driven development to develop a test harness. It hurt my brain too much.)

Anyway, sometime in the spring on 2008 I had something pretty basic working. At about that time, Greg Wilson (original creator of Software Carpentry) visited MathWorks to give a talk about his edited book, Beautiful Code. I met Greg then, and we talked about the needs of MATLAB scientific users. Greg encouraged me to turn my fledgling tool into something real. He also invited me to write up something for the IEEE magazine, Computing in Science and Engineering. That became the article "Automated Software Testing for MATLAB", published in the November/December 2009 issue.

The MATLAB xUnit package raised the visibility of unit testing tools in MATLAB, and its popularity helped make the case for putting something like it in MATLAB. At the same time, our testing tools team was looking to modernize our testing infrastructure. They wanted to take advantage of the latest test framework design principles, and they wanted to decouple the testing machinery from the rest of our complex systems for building and releasing products. So an idea was born: let's make a test framework that can stand on its own, that can test our own software, that can test customer code, and that we can ship in MATLAB!

And now, here it finally is!

Here are a few examples that I've taken from the documentation to illustrate. The first example shows how to write a couple of tests for a hypothetical function called quadraticSolver.

classdef SolverTest < matlab.unittest.TestCase
    % SolverTest tests solutions to the quadratic equation
    % a*x^2 + b*x + c = 0

    methods (Test)
        function testRealSolution(testCase)
            actSolution = quadraticSolver(1,-3,2);
            expSolution = [2,1];
            testCase.verifyEqual(actSolution,expSolution);
        end
        function testImaginarySolution(testCase)
            actSolution = quadraticSolver(1,2,10);
            expSolution = [-1+3i, -1-3i];
            testCase.verifyEqual(actSolution,expSolution);
        end
    end

end

If you need to define setup and teardown methods for your test cases, here's how you do it. (I've just shown the relevant portion of the entire test file.)

properties
    TestFigure
end

methods(TestMethodSetup)
    function createFigure(testCase)
        %comment
        testCase.TestFigure = figure;
    end
end

methods(TestMethodTeardown)
    function closeFigure(testCase)
        close(testCase.TestFigure);
    end
end

Someone asked me just this week about enhancing MATLAB xUnit to provide for setup and teardown methods that get executed just once for all the methods in a test file, instead of being executed once for every individual test method. MATLAB xUnit does not have this capability, but the new unit test framework in R2013a does. A class-level setup method is identified by using the TestClassSetup method attribute, like this:

methods (TestClassSetup)
    function doSomethingBeforeAllTestMethodsAreExecuted(testCase)
        disp('Hi, Mom!')
    end
end

There are lots of new capabilities in this framework. A little birdie told me that more information about matlab.unittest may appear pretty soon in Loren's blog. For now, though, I'd like to point you to the overview video and to the documentation.

That brings me back to MATLAB xUnit and its fate. Well, this little package is entering its "retirement phase." I imagine it will remain on the File Exchange for quite a while because it can take a long time for the user community to adopt a new release, and also because there's a lot of customer-written code out there that uses it. However, I do not plan to work on it anymore. Our testing tools team is an excellent group of engineers, and they can move matlab.unittest forward better and faster than I could ever do with my little side project.

So, if you're interested in unit testing and you get your hands on MATLAB R2013a, please give it a try.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

This is the second part of my plan to rewrite an Image Processing Toolbox example from 20 years ago using more modern MATLAB language features.

Recall the old app I'm trying to reinvent:

The idea is to experiment with the basic principles of DCT-based image compression. I've decided that I don't want to copy this original layout and functionality exactly. For one thing, computers (and MATLAB!) are a lot faster than they were 20 years ago, and I think we really don't need a separate "Apply" button.

Here's my first attempt at sketching what I'd like to make this time around.

The truth, is, I'm making the code and these posts up as I go along. So far I've only had a chance to get some of the initial pieces up and working. But it's enough to start exploring some of the coding ideas.

Here's what happens when I call the app class, which I've called (for now) dctCompressionExample_v2.

app = dctCompressionExample_v2

app = 

  dctCompressionExample_v2 handle

  Properties:
         OriginalImage: [240x240 double]
    ReconstructedImage: []
            ErrorImage: []
    DCTCoefficientMask: []
    NumDCTCoefficients: 10

You can see I got some output in the Command Window, and I also got a new figure on the screen.

Here's the code so far.

classdef dctCompressionExample_v2 < handle

    properties (SetAccess = private)
        OriginalImage
        ReconstructedImage
        ErrorImage
        DCTCoefficientMask
    end

    properties
        NumDCTCoefficients = 10
    end

    properties (Access = private)
        Slider
        OriginalImageAxes
        ReconstructedImageAxes
        ErrorImageAxes
        MaskAxes
    end

    methods
        function this_app = dctCompressionExample_v2
            this_app.OriginalImage = initialImage;
            layoutApp(this_app);
            update(this_app);
        end

        function layoutApp(app)
            f = figure;
            app.OriginalImageAxes = axes('Parent',f, ...
                'Position',[0.12 0.56 0.28 0.28]);
            app.ReconstructedImageAxes = axes('Parent',f, ...
                'Position',[0.60 0.56 0.28 0.28]);
            app.MaskAxes = axes('Parent',f, ...
                'Position',[0.12 0.16 0.28 0.28]);
            app.ErrorImageAxes = axes('Parent',f, ...
                'Position',[0.60 0.16 0.28 0.28]);
            app.Slider = uicontrol('Style','slider', ...
                'Value',app.NumDCTCoefficients/64, ...
                'Units','normalized', ...
                'Position',[0.12 0.08 0.28 0.04], ...
                'Callback',@(~,~,~) app.reactToSliderChange);

        end

        function reactToSliderChange(app)
            v = get(app.Slider,'Value');
            app.NumDCTCoefficients = round(64*v);
            update(app);
        end

        function update(app)
            imshow(app.OriginalImage,'Parent',app.OriginalImageAxes);
            title(app.OriginalImageAxes,'Original Image');

            imshow(app.OriginalImage,'Parent',app.ReconstructedImageAxes);
            title(app.ReconstructedImageAxes,'Reconstructed Image');

            imshow(zeros(size(app.OriginalImage)),'Parent',app.ErrorImageAxes);
            title(app.ErrorImageAxes,'Error Image');

            imshow(zeros(size(app.OriginalImage)),'Parent',app.MaskAxes);
            title(app.MaskAxes,'DCT Coefficient Mask');
            drawnow;
        end
    end

end

function I = initialImage
pout = imread('pout.tif');
pout2 = pout(1:240,:);
I = im2double(adapthisteq(pout2));
end

At this point I'll remind you of Dave Garrison's article that I mentioned last time. This article teaches about how this kind of code works. I don't intend to repeat all of that article here. I'll just point out a few things.

Here's the function that runs when you launch the app.

        function this_app = dctCompressionExample_v2
            this_app.OriginalImage = initialImage;
            layoutApp(this_app);
            update(this_app);
        end

initialImage is a little function I stuck at the bottom of the file. It's purpose is to create our sample image, which for the purpose of this simple DCT demonstration needs to be a multiple of 8-by-8 in size. (Plus, I threw in a call to adapthisteq because the original image has low contrast.)

function I = initialImage
pout = imread('pout.tif');
pout2 = pout(1:240,:);
I = im2double(adapthisteq(pout2));
end

Then there's the properties block.

    properties (SetAccess = private)
        OriginalImage
        ReconstructedImage
        ErrorImage
        DCTCoefficientMask
    end

    properties
        NumDCTCoefficients = 10
    end

    properties (Access = private)
        Slider
        OriginalImageAxes
        ReconstructedImageAxes
        ErrorImageAxes
        MaskAxes
    end

Some of the properties, such as OriginalImage, I want to make accessible but read-only. Some of them, such as OriginalImageAxes, are private because they are only needed by the guts of the app. Notice that one of them, NumDCTCoefficients, is both readable and settable. I wanted to explore the idea of giving this app a little programmable interface so that you can write code to make it change.

I would like to make it so that the NumDCTCoefficients property changes when you drag the slider. Here's how to accomplish that:

            app.Slider = uicontrol('Style','slider', ...
                'Value',app.NumDCTCoefficients/64, ...
                'Units','normalized', ...
                'Position',[0.12 0.08 0.28 0.04], ...
                'Callback',@(~,~,~) app.reactToSliderChange);

Setting the 'Callback' property of the slider like this causes the app's function reactToSliderChange to get called whenever the slider is dragged. Here's what reactToSliderChange does:

        function reactToSliderChange(app)
            v = get(app.Slider,'Value');
            app.NumDCTCoefficients = round(64*v);
            update(app);
        end

Now we're starting to get a little interactivity going in our app.

That's about as far as I got this week. Next time we'll continue wiring things up and maybe start putting in some algorithmic DCT block processing code.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Last month, Dave Garrison wrote a nice article on the MathWorks web site about coding an app (or GUI) by defining a class.

Curiously, this article made me think of a day back in September 1993. I arrived at MathWorks on that day to interview for a job as an image processing software developer. Sometime in the afternoon I found myself in an office with the two developers of version 1.0 of the Image Processing Toolbox. (One of those developers, Loren Shure, was the first employee hired by the MathWorks founders. You can find her over at Loren on the Art of MATLAB.)

Anyway, the Image Processing Toolbox was about to ship, and the developers were trying to think of some demo ideas. I joined in a brainstorming session in front of a whiteboard. One of my ideas that day became the product demo called dctdemo. Here is a screen shot.

(That's a picture of Loren, by the way.)

Prompted by Dave's article, I thought it might be a good time to revisit this oldie-but-goodie demo. It was originally written using MATLAB 4, and the programming techniques available at that time for doing this sort of thing were relatively crude. My plan is to take a fresh look at this demo and reimplement it using the techniques described by Dave.

Here's a start:

classdef dctCompressionDemo < handle

end

That's enough code to actually execute.

dctCompressionDemo

ans = 

  dctCompressionDemo handle with no properties.

And that's about all it does (for now). Next time we'll start filling in the details.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

A couple of weeks ago, Cleve wrote a post about George Forsythe, Cleve's "thesis advisor, colleague, and friend."

In that post, Cleve showed some MATLAB code to load in and display a 1964 picture of the organizing committee of the Gatlinburg conferences on numerical algebra.

load gatlin
image(X)
colormap(map)
axis image
axis off

As Cleve mentioned, this picture is one of the very first images distributed with MATLAB. (I recall getting the first MATLAB version capable of image display sometime around 1990. It was very exciting!)

As it turns out, there is another interesting method in MATLAB to display this picture (or at least a cropped, low-resolution version of it):

default_image = pow2(get(0,'DefaultImageCData'),47);
numbits = 12 - 9 + 1;
b = bitshift(default_image,-9);
b = fix(b);
b = bitand(b,bitcmp(0,4));
b = b/max(b(:));

imshow(b,'InitialMagnification',200)

Forsythe is on the right of the cropped version.

So, what's going on here? Well, it turns out that, if you call image with no input arguments (no data), an image with "default" pixel values gets displayed. Furthermore, that image has many images "hidden" within it. For the full story, see "The Story Behind the Default MATLAB Image."

Cleve knew that I was working on this little hidden-image project in the mid-1990s, and I asked him for suggestions about images to include. He suggested this Gatlinburg Conference picture, as well as the 4-by-4 Durer magic square.

c = bitshift(default_image,-23);
c = fix(c);
c = bitand(c,bitcmp(0,5));
c = c/max(c(:));

imshow(c,'InitialMagnification',200)

Thanks for the suggestions, Cleve!

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Today's post is part of an ongoing tutorial series on digital image processing using MATLAB. I'm covering topics in roughly the order used in the book Digital Image Processing Using MATLAB.

When working with images in MATLAB, it is important to understand how different numeric data types can come into play.

The most common numeric data type in MATLAB is double, which stands for double-precision floating point. It's the representation of numbers that you get by default when you type numbers into MATLAB.

a = [0.1 0.125 1.3]

a =

    0.1000    0.1250    1.3000

class(a)

ans =

double

Double-precision floating-point numbers are intended to approximate the set of real numbers. To a reasonable degree, one can do arithmetic computations on these numbers using MATLAB (and the computational hardware on your CPU or GPU) and get the same results as "true arithmetic" (or "God's math," as I've heard Cleve say) on the real numbers.

Working with floating-point numbers is very useful for mathematical image processing algorithms (such as filtering, Fourier transforms, deblurring, color computations, and many others).

Prior to 1997, the double was the only kind of data type in MATLAB. Image processing customers complained about this because of the memory required for these kinds of numbers. A double-precision floating-point number requires 64 bits, whereas many people working with image data were used to using only 8 bits (or even just 1 bit in the case of binary images) to store each pixel value.

So with MATLAB 5 and Image Processing Toolbox 2 in 1997, we introduced support for a new data type, uint8, which is an abbreviation for unsigned 8-bit integer. This data requires just 8 bits to represent a number, but the representable set of numbers is limited to the integers from 0 to 255.

You can make one in MATLAB by calling the uint8 function.

b = uint8(5)

b =

    5

class(b)

ans =

uint8

Also, you often see uint8 numbers when you call the imread to read an image from a file. That's because image file formats often use 8 bits (prior to compression) to store each pixel value.

rgb = imread('peppers.png');
rgb(1:3,1:4,1)

ans =

   62   63   63   65
   63   61   59   64
   65   63   63   66

class(rgb)

ans =

uint8

Almost immediately after MATLAB 5 and Image Processing Toolbox 2, we started hearing from customers who had scientific data stored using 16 bits for value, so 8 bits wasn't enough and 64 bits (for double) still seemed wasteful. So Image Processing Toolbox 2.2 in 1999 added support for uint16 numbers (unsigned 16-bit integers).

But still that wasn't enough. The medical imaging community, it seemed, needed signed 16-bit numbers. And, said many, what about single-precision floating-point?

For Image Processing Toolbox 3 in 2001, we stopped adding data type support piecemeal and instead added support for all the data types in MATLAB at the time. Here is a summary of the entire set:

Support for the logical data type (the only values are 0 and 1, 1 byte per number) was added a few years later.

Two other data types have appeared since then, uint64 and int64. Relatively little effort has been made to support these data types for image processing, for two reasons:

c = uint64(184467440737095516)

c =

   184467440737095516

d = double(c)

d =

   1.8447e+17

e = uint64(d)

e =

   184467440737095520

e - c

ans =

                    4

When I pursue this topic further next time, I'll talk more about data type conversions: the basic ones in MATLAB, plus the Image Processing Toolbox ones that handle additional details of data scaling.

For more information, see Section 2.5 of Digital Image Processing Using MATLAB.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b