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.

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.....