So you already know that I ♥ my pancreas. I hope that soon I’ll be able to ♥ both my real and “artificial” pancreases. Kerri at Six Until Me has a great write-up of an announcement between JDRF, Johnson & Johnson, and DexCom which should speed my new love’s arrival. (Update: The Diabetes Mine has even more details and an interview.)

And I ♥ my job. I won’t try to take credit for any of the hard work that JDRF-funded scientists and industry R&D folks are doing to make an artificial pancreas a reality. They deserve all of the credit and so much more. But I know for certain that many people at those institutions and others are using our tools as part of the day-to-day toil of making this a reality. (MATLAB’s graphics have a very *ahem* distinctive appearance, which I’ve seen in presentations about the artificial pancreas and diabetes self-management findings.)

The MathWorks‘ mission statement is that we “accelerate the pace of engineering and science worldwide.” Although it might just sound like nice words on a web site or T-shirt, at times like this it makes me very proud to know that we really do make great things happen faster. And knowing that — even as a peripheral actor — I help with progress, well, that inspires me to do my job that much better, by putting a little more care and polish into my features (especially those related to medical imaging), by rooting out all of the hidden software defects that I might otherwise overlook, and by working hard to get as many of those customer-requested features into the product that I can. That may sound a little corny, but I know that everything I do to accelerate the development of our products helps other people accelerate the creation of their products. And I’m not the only one in the office who feels this way.

Now, back to (indirectly) helping other people fight the good fight. . . .

read more from this topic.....

A couple weeks ago, I presented a teaser of my diabetes application. I’ve been working on it a lot lately, creating GUIs to enter most of the data on my insulin pump and blood glucose meter. But mostly I’ve been creating charts and reports. (Putting data in is much less interesting than analyzing it, in an effort to improve my “control.”)

Here’s one of the reports. It has a bit of everything, including the kitchen sink. The one I’m taking with me to my endocrinologist tomorrow is much more focused. Visualizing all of the data associated with diabetes is kind of tricky, so I’m working on that. More features to come. . . .

(Please, don’t judge me based on my readings. :^) Right now I’m using food to cover my insulin, which is the opposite of how it should work — it seems I never really got switched from the NPH mindset. Getting past that is one of the main reasons why I’m doing in this project.)

I’m committed to releasing everything under some kind of open source license. Everything, that is, except the data and anything that would involve a diagnosis or recommend a particular kind of treatment. I don’t want to be on the hook for FDA approval and all that involves.

read more from this topic.....

Just to get this disclaimer out of the way, let’s remind everyone that this is my personal website. Even though I may write some of the articles at work about things I’ve done at work, they aren’t in anyway supervised or endorsed by my employer. I’m only including this disclaimer because this article is meant to be self-deprecating in that “I had the best of intentions, but now I feel a little silly” way.

Several years ago I read Greenspun’s Tenth Rule of Programming: “Any sufficiently complicated C or Fortran program contains an ad-hoc, informally-specified, bug-ridden, slow implementation of half of Common Lisp.” Basically it’s a joke about feature creep in software. I was reminded of this quotation today when I started reading Joshua Kerievsky’s Refactoring to Patterns. He describes the perils of over-engineering, which is basically the same as wasting money.

You see, about seven or eight years ago, I added a Lisp parser to a feature in MATLAB’s Image Processing Toolbox. If you have the toolbox, just take a look in toolbox/images/iptformats/private/dicom_create_IOD.m. There you’ll see a part of the code that looks like this:

function tf = check_condition(expr, metadata)
%CHECK_CONDITION  Determine whether a particular condition is true.
%
%   Conditions are LISP-style cell arrays.  The first element is a
%   conditional operator, remaining cells are arguments to the operator.
%
%   Conditions can be nested.  Each cell array in expr indicates a new
%   conditional expression.

[snip]

%
% Process conditional expressions recursively.
%
switch (lower(operator))
case 'and'

    % This AND short circuits.
    for p = 1:numel(operands)
        tf = check_condition(operands{p}, metadata);

        if (~tf)
            return
        end
    end

[snip]

So why did I do this? The DICOM file format has conditional elements, which only show up in the file if other elements are present, are absent, or have a particular value. I wanted to create a general-purpose, reusable, extensible mechanism for describing these conditions so that I didn’t actually have to write code for all of those conditionals. Moreover, I wanted our customers to be able to write the conditions in a nice tabular form rather than writing code, because — you know — one day they would want to write more of the dozens of kinds of DICOM information objects from scratch than we supported right out of the box.

You probably see where this is going. We ended up extending the product in a different direction after we learned more about how our customers actually wanted to use the DICOM export functionality. So, if you need a Lisp-like algebraic Boolean logic parser in MATLAB, just go look in the previously mentioned file. (Also take a look at the dicom_modules.m file in the same directory to see how to express those conditions.) But for g-d’s sake, ask yourself if you really, really need that, or if you aren’t just fooling yourself into over-engineering your code.

read more from this topic.....

I am attending Electronic Imaging 2009 in San Jose. Well, today was my last day here. Tomorrow I return to snowy Massachusetts. (It’s currently about 60º here.) But I’ve had my fish taco fix for another year, and I’m excited to get back to Lisa and the kitty.

As usual, I learned something new and talked to a lot of people about imaging and MATLAB. Last year, the new-to-me things concerned image quality for mobile devices and image forensics. This year’s big surprises were 3-D imaging, which apparently is going mainstream. (According to Andrew Woods, the chair of the Stereoscopic Displays and Applications conference, there are already two million televisions in people’s homes which are capable of displaying 3D video.) Every year that I attend, I seem to get a different view of the conference; so I’m cautious to read too much into any given year. But somehow, 3D felt different this time.

The exhibitors behind us were demonstrating their more advanced 3D displays and video capture systems, and the person across from us was selling 3D notecards and prints. So my coworker and I thought we’d get in on the action during the waning hours of the exhibit. We decided to try our hand at making an anaglyph, one of the easier kinds of 3D images to make (along with stereograms, which are a bit lo-fi for MATLAB, after all).

The recipe calls for combining two images of the same scene from slightly different perspectives using a simple function that combines the red channel from the left image and the green and blue channels from the right image. And of course, you need the funny glasses. First, the images that we made with my camera phone:

Left image

Right image

Using some MATLAB code that we found on MATLAB Central, you get the following image:

The anaglyph

(Unfortunately, I can’t help you with the funny glasses. Go get yer own.)

I am told the results are “okay for a cha-cha picture.” I don’t have binocular vision, so what do I know?

Have fun. (And don’t be like the guy who said that the 1.57% of the population with monocular vision aren’t really worth worrying about.)

read more from this topic.....

Peter Webb, one of my coworkers, published an article entitled “Rendering High Dynamic Range Images on the Web” in the July 2008 issue of the MATLAB Digest. Enjoy!

read more from this topic.....

Comrade programmers!

About five years ago I wrote a little program to help me find numeric patterns in arrays. Basically I needed to find byte patterns like 0xFFFE 0xE000 in the data from DICOM files. Something like MATLAB’s strfind function, except more general than looking for substrings in strings.

So I wrote a rather naïve algorithm to do just that. It wasn’t very fast, but it worked for exploratory use. But then today I needed to do something like that in production code. So I asked one of my coworkers for tips with a faster implementation.

Turns out the strfind works just as well on numeric arrays as it does on strings. And it’s really fast, too.

% "data" is a 1,200,000 element UINT8 array.
% It contains 29 instances of the pattern.
tic; offsetTable = strfind(data, [254 255 0 224]); toc

Elapsed time is 0.027390 seconds.

Update: Loren Shure posted an even better solution on her blog.

read more from this topic.....

Here’s a little something for the MATLAB lovers out there who also program in C++. The code sample at the end of this post shows how to add a C++ iterator to an mxArray, which is the basic datatype for MATLAB arrays. This makes it harder to walk off the end of an mxArray during linear traversal of all of the elements in an array and (hopefully) improves the readability of code that accesses the data stored inside mxArrays. The code gives a couple simple-minded examples of element-wise traversal (provided the MATLAB array you give it is a UINT8 array).

If you’re familiar with the C++ Standard Library, you can start to use your wrapped mxArray like an STL container. I use the term “like” advisedly, since it’s not actually a container. To reduce memory overhead, I don’t make a copy of the data that gets put into ArrayWrapper. Consequently, it doesn’t model the container concept that element lifetimes are the same as the “container” lifetime. You’re shouldn’t count on being able to use ArrayWrapper with any ole STL function. But then again, you might be able to use it with some; I haven’t tried.

#include "mex.h"
#include <cstddef>

// Treat an mxArray like a container.
// (See Josuttis. "The C++ Standard Library." 1999. p. 219-220)
template<class T>
class ArrayWrapper {
  private:
    T      *v;
    mwSize thesize; 

  public:
    // Type definitions
    typedef T         value_type;
    typedef T*        iterator;
    typedef const T*  const_iterator;
    typedef T&        reference;
    typedef const T&  const_reference;
    typedef mwSize    size_type;
    typedef ptrdiff_t difference_type;

    // Constructor
    ArrayWrapper(const mxArray *theArray)
    {
        v = static_cast<T*>(mxGetData(theArray));
        thesize = mxGetNumberOfElements(theArray);
    }

    // Iterator support
    iterator begin() { return v; }
    const_iterator begin() const { return v; }
    iterator end() { return v + thesize; }
    const_iterator end() const { return v + thesize; }

    // Direct element access
    reference operator[](size_type i) { return v[i]; }
    const_reference operator[](size_type i) const { return v[i]; }

    // Size is constant
    size_type size() const { return thesize; }
    size_type max_size() const { return thesize; }
};

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    const mxArray *inputArray = prhs[0];

    ArrayWrapper<uint8_T> inData(inputArray);

    // 1- Count the number of samples in the image.
    mwSize count = 0;
    ArrayWrapper<uint8_T>::const_iterator i = inData.begin();
    while (i++ != inData.end())
    {
        count++;
    }

    mexPrintf("The image has %ld samples.\n", count);

    // 2 - Halve the image values.
    mxArray *outputArray = mxDuplicateArray(prhs[0]);
    ArrayWrapper<uint8_T> outData(outputArray);

    i = inData.begin();
    ArrayWrapper<uint8_T>::iterator j = outData.begin();
    while (i != inData.end())
    {
        *(j++) = *(i++) * 0.5;
    }

    plhs[0] = outputArray;
}

read more from this topic.....

Avoid str2double and str2num. Use sscanf instead.

For scalars, you’ll see a modest improvement.

>> str = '0009';
>> tic; for p=1:1000, str2num(str); end; toc
Elapsed time is 0.126388 seconds.
>> tic; for p=1:1000, sscanf(str, '%d'); end; toc
Elapsed time is 0.022299 seconds.

>> str = '3.14159265';
>> tic; for p=1:1000, str2double(str); end; toc
Elapsed time is 0.056466 seconds.
>> tic; for p=1:1000, sscanf(str, '%f'); end; toc
Elapsed time is 0.017805 seconds.

For vectors, you’ll see a more hefty speed up.

>> str = '0009 3.14159265';
>> tic; for p=1:1000, str2double(str); end; toc
Elapsed time is 0.480512 seconds.
>> tic; for p=1:1000, sscanf(str, '%f'); end; toc
Elapsed time is 0.027449 seconds.

Also favor sprintf instead of num2str.

>> tic; for p=1:1000, num2str(p); end; toc
Elapsed time is 0.101599 seconds.
>> tic; for p=1:1000, sprintf('%d', p); end; toc
Elapsed time is 0.018325 seconds.

read more from this topic.....

The CDF folks at Goddard Space Flight Center have identified a security vulnerability — a buffer overflow to be specific — that can enable the execution of arbitrary code on your machine if you open a particular malformed file. If you’re accessing CDF files via MATLAB, you can download a security patch from NASA GSFC.

Thank you. That is all.

UPDATE: You can also download an update directly from The MathWorks.

read more from this topic.....

Note: Here’s another dispatches about what happened at the Electronic Imaging symposium even though I’m back from San Jose. They will continue until I run out of useful things to write or something more interesting happens in my life.

Yasuyo Ichihara of Kogakuin University presented some research that she and her colleagues performed on what many of us call “color blindness” or “color deficiency.” She would prefer that we say these differently sighted folks have “color confusion,” which I’ll probably try since it’s slightly more evocative than “color deficiency.” And from a PC point of view, it may not be nice to stagmatize roughly 5% of the population as deficient, etc.

Anyway, they presented research of how their lab applied color universal design to remake the map and timetables of the Tokyo subway system to be friendlier for people with both protanopia and deuteranopia. It turns out that the color confusion can be quite simply modeled by placing two iso-chromaticy lines on the standard x-y chromaticity diagram. Along those lines, all colors appear to be the same hue; so it’s not too hard to pick colors that are maximally contrasting by avoiding the iso-lines. They picked four colors, which they claim are suitable for many print applications: black, orange-red, bluish-green, and a more pure blue.

Color universal design doesn’t stop at picking the right colors. Adding non-color information — such as letters, underlining, framing, or other symbols — proved to be very useful aids to comprehension, too.

All of this reminded me that awhile back I implemented a color deficiency simulator in MATLAB based on work by Hans Brettel. I’m making these functions available here in the hopes that they’re useful and that you actually use them when designing something colorful. Simply download the protanopia and deuteranopia simulator M-files and the associated data file. (You will need the Image Processing Toolbox, since these functions use rgb2ind.)

Let’s see what the peppers image I created some years ago looks like to people with color confusion:

rgb = imread('peppers.png');
figure; imshow(rgb); title('Normal color vision')
figure; imshow(deuteranopia(rgb)); title('Deuteranopia')
figure; imshow(protanopia(rgb)); title('Protanopia')

read more from this topic.....

Every once and a while someone will ask me why I still use film. First of all, I have a film camera, and to buy a new one would be a significant chunk of cash. But that’s not the main reason. I like the look of film; I know how the various emulsions work and know what I need to do to get the best image from them; and (above all) I like the fact that I have an artifact when I get my film back from the lab.

Digital does have a lot of advantages, though. It provides instant feedback at the point of making the photograph. It’s also hard to beat the speed and flexibility of a digital workflow. And for 35mm-sized sensors, it’s hard to tell the difference from film.

This article mixes the look of film and the enormous possibilities of digital. It also continues the trend of notes about high dynamic range images that I’ve been posting recently.

Color film is an analog tone reproduction operator. It renders a vast range of illumination values to a particular color image. That is, it maps the high dynamic range world to a lower dynamic range representation. The major factors that influence this tone mapping “function” include the exposure settings (shutter time, aperture, and ISO speed), the film response curves, and development technique. A film response curve for Fuji Velvia, the color film I use, appears below. Notice how the curve is more-or-less flat in the highlights and shadows and that it has an S-shaped curve in between.

If you take some liberties with that S-shaped part of response curve where the real detail in the image is, you might almost say that the curve has a flat part for the shadows, a flat part for the highlights, and a line that connects the two.

When I first started puttering around with high dynamic range images, I began by translating what I knew about how film renders a scene with a lot of tones:

• Film is sensitive across a particular exposure range.
• Parts of the image with less than a certain amount of exposure are rendered as black.
• Overexposed parts of image are rendered as white (technically clear).
• You can’t change the film sensitivity (except with some processing tricks that we’ll ignore).
• You can change the amount of light that is recorded during the exposure.

It wasn’t hard for me to combine this information with the generalization about the film response “curve” to produce a very simple tone-renderer which simulates color slide film in MATLAB. You can download the M-file from my page on MATLAB Central, where you can find some other good stuff. (I’ve also reprinted it below.) I ran it on the HDR image of my office, varying the exposure from 0 EV to 12 EV in half-stop increments, and was very glad to see the same kind of results that I do with film.

(Click for larger image…)

But film is so twentieth-century. So I went on to more interesting forms of tone mapping.

function rgbSimulated = simulateFilm(rgbRadiance, nStops)
%simulateFilm   Perform film-like tone mapping.
%   LDR = simulateFilm(HDR, middleEV) converts the floating-point high
%   dynamic range image HDR to a UINT8 low dynamic range image LDR using a
%   method that recreates the sensitivity of film.  The middleEV value sets
%   the middle exposure value (EV) for the rendering.  Larger EV numbers
%   will generate brighter images; smaller values of middleEV will result
%   in darker images.
%
%   The film sensitivity is set to 6 EV, which is comparable to many
%   late-model transparency (positive) films.
%
%   Example
%   -------
%
%   Render the same HDR image using different "exposure settings."
%   Essentially, simulate exposing the same scene from 0 to 12 EV in 1/2
%   stop steps.
%
%      % Read the HDR image and create a buffer for the rendered images.
%      hdr = hdrread('office.hdr');
%      s = size(hdr);
%      ldr = ones(s(1), s(2), 3, 25, 'uint8');
%
%      % Perform the tone mapping.
%      for p = 0.5:0.5:12.5
%          ldr(:,:,:,p*2) = simulateFilm(hdr, p);
%      end
%      figure; montage(ldr)

% Author: Jeff Mather
% Copyright 2006-2007 The MathWorks, Inc.

% Set the film sensitivity and the mid-tone log-radiance.
sensitivity = 6;
midPoint = 3;

minLogExposure = midPoint - sensitivity/2;
maxLogExposure = midPoint + sensitivity/2;

% Film works in stops.  Convert radiance to base 2 to compute perception.
% Values that are outside the film sensitivity are lost (min or max).
rgbRadiance = rgbRadiance * 2^(nStops);

rgbSimulated = rgbRadiance;
rgbSimulated(rgbRadiance ~= 0) = log2(rgbRadiance(rgbRadiance ~= 0));
rgbSimulated(rgbSimulated < minLogExposure) = minLogExposure;
rgbSimulated(rgbSimulated > maxLogExposure) = maxLogExposure;

% Convert to RGB.
rgbSimulated = uint8(255 * (rgbSimulated - minLogExposure) ./ ...
                     sensitivity);

read more from this topic.....

How does the tonemap function in the Image Processing Toolbox product work?

Tone mapping, or tone reproduction, compresses the enormous amount of illumination data in a high dynamic range image to something more suitable for output on a medium that has a lower dynamic range. The recent discussion of my office image shows the results of tone mapping.

But what is dynamic range? And what makes a high dynamic range image?

In the simplest terms, dynamic range is the ratio of the brightest value in an image or scene and the darkest value. Consider the following table of illumination sources from High Dynamic Range Imaging by Erik Reinhard, Greg Ward, Sumanta Pattanaik, and Paul Debevic:

The human visual system is capable of discriminating about five orders of magnitude at most of these different illumination levels.[1] That is, if you have normal vision, you can see something with a brightness value of 1 and a value of 100,000 in the same scene. Under different conditions, you might be able to make out detail in something with a value of 0.001 at the same time that you look at an object with a brightness level of 100. What you can’t do is see detail in something with a value of 0.001 and another thing with brightness 100,000 simultaneously. That would require being able to distinguish between eight orders of magnitude at once.

An image with values that range from 1 to 100,000 has a dynamic range of 100,000:1. Most of the images made for display on contemporary monitors have a dynamic range of only 256:1 per color channel, because that’s all that most monitors are built to support. We’ll call the latter images low dynamic range (LDR) images and anything with the ability to store more illumination detail a high dynamic range (HDR) image. There’s nothing stopping us from creating HDR images, but we won’t be able to display them on a LDR device like a monitor or inkjet prints.

Tone mapping is the HDR to LDR conversion. There are many tone mapping operators that we might have implemented. (See “A review of tone reproduction techniques” (2002) by Kate Devlin for an overview of these technique.) We chose a technique that renders floating-point high dynamic range RGB radiance maps to low dynamic range RGB UINT8 images using a spatially uniform tone mapping operator similar to what Ward et al. described in their 1997 paper “A visibility matching tone reproduction operator for high dynamic range images” (IEEE Transactions on Visualization and Computer Graphics, 3(4):291-306). It differs in several ways, though.

In a nutshell, our tonemap function takes the base-2 logarithm of the RGB radiance map to mimic the human visual response to brightness, transforms the image into the L*a*b* colorspace, and performs contrast limited adaptive histogram equalization on the luminance channel without changing the overall color content of the image. After converting the modified L*a*b* image back to RGB, the result for most images is a fairly good LDR image suitable for output. We provide a couple of parameters to improve the “artistic” qualities of the rendering – AdjustSaturation and AdjustLuminance – as well as the ability to tune the histogram equalization.

While our technique doesn’t use much of what is known about how the human visual system performs tone mapping, we think we picked a method that balances good results, acceptable performance, and low complexity. My hope is that people doing work on HDR imagery will take a look at the two-dozen so important lines inside tonemap and see just how easy it is to get started implementing their own tone reproduction operator in MATLAB.

[1] – One “order of magnitude” is one power of ten. Engineering types are fond of talking about orders of magnitude.

read more from this topic.....

Curious how I spend my days at work? Well, my new article Two New Functions for Converting Datatypes and Changing Byte Order pretty much sums it all up. Create small applications in MATLAB and C to access file formats. Write tools to make working with those files easier. Write the occasional article about tool-building for others toolmakers.

Oh! And meetings. Can’t forget those. . . .

My other recently published article, Color Management and Color Transformations in MATLAB, shows the work that isn’t related to formats that I do at the office. Currently I’m somewhere between that stage where I think I know what I’m talking about and actually knowing it.

read more from this topic.....

“The next-generation engineer will use the 3D control scheme he has grown up with.” That’s how Jörg Buchholz introduced his Doom 1.1 MATLAB code that lets data visualizers “fly through a 3D scene like in a first-person shooter in god mode.”

Are shooter games the future of data visualization? Maybe not, but the entertainment industry will play a role.

• Radiologists are awash in so much data than they have trouble reading “films” in the traditional way. SCAR’s Transforming the Radiological Interpretation Process (TRIP) initiative draws upon expertise from Hollywood for—among other things—novel visualization techniques and data navigation tools.
• When I visited the Chicago Mercantile Exchange on a break during a business trip, I saw a trader using what looked like a PlayStation game controller. A Baseline Magazine article explains why game controllers are useful for trading: “[The Exchange is pursuing] new markets such as trading arcades, small venues in Dublin, London and Gibraltar that host a new generation of traders who can use joysticks instead of computers to swap financial derivatives. ‘They trade as if it were a video game,’ says Steve Goldman, director of network architecture at the exchange, holding up a PlayStation joystick. ‘And we want to make sure they have something to trade here.’ . . . These traders follow set systems and algorithms and chase seemingly minuscule market moves. Under these systems, rapid hand-eye coordination is critical.”

read more from this topic.....

Earlier this week The MathWorks published my article Color Management and Color Transformations in MATLAB. Enjoy!

read more from this topic.....