Woo hoo! The Image Processing Toolbox team has just created a new product:Medical Imaging Toolbox™Shipped with the R2022b release a couple of months ago, this product provides apps, functions,... read more >>

]]>Woo hoo! The Image Processing Toolbox team has just created a new product:

Shipped with the R2022b release a couple of months ago, this product provides apps, functions, and workflows for designing and testing diagnostic imaging applications. You can perform 3D rendering and visualization, multimodal registration, and segmentation and labeling of radiology images. The toolbox also lets you train predefined deep learning networks (with Deep Learning Toolbox). I'm looking forward to writing about this product and its capabilities.

I've been talking recently with Sailesh, one of the product developers, about the Medical Image Labeler. This app is for labeling ground truth data in 2-D and 3-D medical images. With this app, you can:

- Import multiple 2-D images or 3-D image volumes.
- View images as slice planes or volumes with anatomical orientation markers and scale bars.
- Create multiple pixel label definitions to label regions of interest. Label pixels using automatic algorithms such as flood fill, semi-automatic techniques such as interpolation, and manual techniques such as painting by superpixels.
- Write, import, and use your own custom automation algorithm to automatically label ground truth data.
- Export the labeled ground truth data as a groundTruthMedical object. You can use this object to share labels with colleagues or for training semantic segmentation deep learning networks.

The Medical Image Labeler supports 2-D images and image sequences stored in the DICOM and NIfTI file formats. An image sequence is a series of images related by time, such as ultrasound data. The app supports 3-D image volume data stored in the DICOM (single or multifile volume), NIfTI, and NRRD file formats. See Get Started with Medical Image Labeler.

I asked Sailesh for his thoughts about what to tell people in this blog about Medical Image Labeler. He commented that the app is intended for people working towards any kind of AI-assisted CAD (computer-aided diagnosis). It doesn't replace a doctor or diagnostician; rather, it assists with tasks related to their research.

Sailesh also mentioned some specific capabilities that he'd like to highlight. I've prepared a 3-minute video of one of those capabilities: Using active contours to automate object labeling. This tool is one of several in Medical Image Labeler that are intended to save time in the labeling process. Check it out:

There are a few other things about the Medical Image Labeler that I hope to show you soon.

If you have suggestions for the Medical Image Labeler, or for Medical Imaging Toolbox in general, add a comment below.

]]>In the R2022b release, Image Processing Toolbox includes several new geometric transformation objects, such as rigidtform2d, affinetform3d, etc., that use the premultiply matrix convention instead of... read more >>

]]>In the R2022b release, Image Processing Toolbox includes several new geometric transformation objects, such as rigidtform2d, affinetform3d, etc., that use the premultiply matrix convention instead of the postmultiply convention. Several other related toolbox functions, such as imregtform, now preferentially use these new objects. There are functions in the Computer Vision Toolbox and Lidar Toolbox that now use the premultiply convention and new objects. These changes also improve the design consistency with Robotics System Toolbox and Navigation Toolbox.

For the key information in the documentation, see:

- Matrix Representation of Geometric Transformations
- Migrate Geometric Transformations to Premultiply Convention

In today's post, I'll explain how and why this all came about and what difference it makes to users.

Table of Contents

The matrices in question define affine and projective transformations, including translations, rotations, rigid, and similarity transformations. I'll focus on affine transformations in the following discussion, but the same concepts apply to projective transformations.

For two dimensions, an affine transformation matrix is a $ 3 \times 3 $ matrix that maps a two-dimensional point, $ (u,v) $, using matrix multiplication, like this:

$$\left[\begin{array}{c}\mathit{x}\\ \mathit{y}\\ 1\end{array}\right]=\left[\begin{array}{ccc}\mathit{a}& \mathit{b}& \mathit{c}\\ \mathit{d}& \mathit{e}& \mathit{f}\\ 0& 0& 1\end{array}\right]\times \left[\begin{array}{c}\mathit{u}\\ \mathit{v}\\ 1\end{array}\right]=\mathit{A}\times \left[\begin{array}{c}\mathit{u}\\ \mathit{v}\\ 1\end{array}\right]$$

When the affine transformation is written as above, the third row of A is always $\left[\begin{array}{ccc}0& 0& 1\end{array}\right]$. Because the matrix appears before the vector it is multiplying, I'll call this the premultiply convention.

There's another way to write this operation. You can transpose everything, like this:

$$\left[\begin{array}{ccc}\mathit{x}& \mathit{y}& 1\end{array}\right]=\left[\begin{array}{ccc}\mathit{u}& \mathit{v}& 1\end{array}\right]\times \left[\begin{array}{ccc}\mathit{a}& \mathit{d}& 0\\ \mathit{b}& \mathit{e}& 0\\ \mathit{c}& \mathit{f}& 1\end{array}\right]=\left[\begin{array}{ccc}\mathit{u}& \mathit{v}& 1\end{array}\right]\times \mathit{B}$$

In this form, the matrix appears after the vector, and so I'll call this the postmultiply convention. Note that A and B are related by matrix transposition: $ A = B^{T} $, and $ B = A^{T} $.

The first Image Processing Toolbox release that included general-purpose geometric image transformation functions was developed from about 1999 to 2001. At that time, many of the most useful publications discussing the geometric transformation of images were in the computer graphics literature. Both of the matrix conventions, premultiply and postmultiply, were in use. Which convention you learned depended on which books and papers you read, or which graphics software framework you used, such as OpenGL or DirectX.

I was influenced at the time by the book Digital Image Warping, by George Wolberg, which used the postmultiply convention. I also thought that the postmultiply convention worked well with the usual MATLAB convention of representing P two-dimensional points as a $ P \times 2 $ matrix.

Because of these influences, the initial toolbox functions, maketform and imtransform, as well as the next generation of functions, imwarp and affine2d and others, used the postmultiply convention.

Within a few years, it became apparent that this design choice was confusing people. I mentioned the confusion in my 07-Feb-2006 blog post on affine geometric transformations.

In the years since 2001, the premultiply convention has become far more widely used than the postmultiply convention. The most popular sources of information, such as Wikipedia, use the premultiply convention. As a result, our use of the postmultiply convention was confusing many more people. We could see this confusion in many posts on MATLAB Answers, as well as in tech support requests. Developers on Image Processing Toolbox and Computer Vision Toolbox teams concluded that we should try to do something about it, even though it was likely to be difficult and time-consuming. The design effort began in spring 2021. The final push, in the winter and spring of this year, was a group effort (see below) involving developers, writers, and quality engineers on multiple teams.

The R2022b release of Image Processing Toolbox introduces these new types:

- projtform2d - 2-D projective geometric transformation
- affinetform2d - 2-D affine geometric transformation
- simtform2d - 2-D similarity geometric transformation
- rigidtform2d - 2-D rigid geometric transformation
- transltform2d - 2-D translation geometric transformation
- affinetform3d - 3-D affine geometric transformation
- simtform3d - 3-D similarity geometric transformation
- rigidtform3d - 3-D rigid geometric transformation
- transltform3d - 3-D translation geometric transformation

We encourage everyone to start using these instead of the earlier set of types: projective2d, affine2d, rigid2d, affine3d, and rigid3d.

When you make one of the new transformation types from a transformation matrix, use the premultiplication form. For an affine matrix in premultiplication form, the bottom row is $\left[\begin{array}{ccc}0& 0& 1\end{array}\right]$.

A = [1.5 0 10; 0.1 2 15; 0 0 1]

tform = affinetform2d(A)

tform.A

To ease the transition, the new types are intended to be interoperable, as much as possible, with code that was written for the old types. For example, let's take a look at the old function, affine2d:

T = A'

tform_affine2d = affine2d(T)

tform_affine2d.T

For the old types, the T property is the transformation matrix in postmultiply form. For the new types, the A property is the transformation matrix in premultiply form.

Although it is hidden, the new types also have a T property, and this property contains the transformation matrix in postmultiply form.

tform

tform.A

tform.T

This hidden property is there so that, if you have code that gets or sets the T property on the old type, you will be able to use the new type without changing that code. Setting or getting the T property will automatically set or get the corresponding A property.

tform.T(3,1) = 100;

tform.T

tform.A

In addition to the generic affine transformation, the new types include the more specialized transformations translation, rigid, and similarity. You can create these using parameters that may be more intuitive than the affine transformation matrix. For example, a rigid transformation is a combination of rotation and translation, and so you can create a rigidtform2d object by specifying a rotation angle (in degrees) and a translation vector directly.

r_tform = rigidtform2d(45,[0.2 0.3])

If you ask for R (the rotation matrix) or A (the affine transformation matrix), it is computed directly from the rotation and translation parameters.

r_tform.R

r_tform.A

You can change these matrices directly, but only if the result would be a valid rigid transformation. The following assignment, which only changes the horizontal translation offset, is permitted because the result is still a rigid transformation:

r_tform.A(1,3) = 0.25

But if you try to change the upper-left $ 2 \times 2 $ submatrix so that it is not a rotation matrix, you'll get an error:

Transforming points and images works in the same way as with the old objects.

[x,y] = transformPointsForward(r_tform,2,3)

[u,v] = transformPointsInverse(r_tform,x,y)

The following code generates 100 pairs of random points, transforms them using r_tform from above, and then plots line segments from the original points to the transformed ones.

xy = rand(100,2) - 0.5;

uv = transformPointsForward(r_tform,xy);

clf

hold on

for k = 1:size(xy,1)

plot([xy(k,1) uv(k,1)],[xy(k,2) uv(k,2)])

end

hold off

axis equal

And imwarp interprets the new transformation types using the same syntax as before.

A = imread("peppers.png");

B = imwarp(A,r_tform);

imshow(B)

These geometric transformation objects are widely used in several MathWorks toolboxes. There are currently 20 or so different documentation examples that use rigidtform3d. The examples are from these products:

- Image Processing Toolbox
- Computer Vision Toolbox
- Automated Driving Toolbox
- Lidar Toolbox

Here's a sampling:

It took a big group effort, earlier this year, to make all the changes in Image Processing Toolbox, Computer Vision Toolbox, Automated Driving Toolbox, and Lidar Toolbox. The Image Processing Toolbox design and implementation, which I worked on, was relatively straightforward, but the Computer Vision Toolbox changes were extensive and quite complicated. Thanks go to Corey, Paola, and Qu for their implementation and design work. (Corey, I'm sorry that this project swallowed up your entire internship! I'm glad that you have officially joined the development team now.) Witek leaped in to help revise the Computer Vision Toolbox designs to incorporate lessons learned over the past several years. (Witek is co-author of the upcoming 2023 edition of Robotics, Vision, and Control: Fundamental Algorithms in MATLAB, which uses the premultiply convention.) From the Lidar Toolbox team, Kritika helped with design and implementation, and Hrishikesh updated some examples.

Alex, thanks for being my Image Processing Toolbox design buddy; your experience was invaluable. Vignesh, thanks for advising me about code generation. Ashish, thanks for rescuing me at the very last minute with critical implementation help. Jessica did a wonderful job with the extensive documentation and example updates across all four products. And Abhi helped stage and qualify the final multiproduct integration under a tight deadline.

THANKS TO ALL!

]]>

About a month ago, I posted about my Code Trace for MATLAB tool, which is posted to GitHub and the File Exchange. This tool is useful for exploring and troubleshooting executing MATLAB code. I... read more >>

]]>I just published an update, version 2.0.0. There are now several new options for addCodeTrace, including:

- Evaluate code trace expressions with no output arguments

NumExpressionOutputs = 0 or 1 - Write code trace displays to file instead of showing them in the Command Window

OutputFile = filename - Suppress the stack depth indentation

ShowStackDepth = true or false - Suppress code trace displays

DisplayTrace = true or false

Note that version 2.0.0 no longer does stack depth indentation by default.

See my earlier post for more information about this tool. ]]>Back in June, my answer to a question on MATLAB Answers used a technique called double thresholding, also known as hysteresis thresholding. Although the question raised was a 1-D problem, I wanted to... read more >>

]]>Back in June, my answer to a question on MATLAB Answers used a technique called double thresholding, also known as hysteresis thresholding. Although the question raised was a 1-D problem, I wanted to share with you the general idea in two dimensions, making use of the function imreconstruct.

Double thresholding can be useful when a single threshold produces unsatisfactory results because it either includes too many pixels in the foreground (with a low threshold), or it includes too few (with a high threshold). The idea is to use two thresholds, one that is more selective and one that is less selective, and then let the more selective result "grow" by adding touching pixels from the less selective result.

Here's a simple, contrived example to illustrate what I mean.

A = mat2gray(peaks(200));

imshow(A)

Suppose I want to compute a mask of the peaks that exceed 0.9:

T_high = 0.9;

mask_high = A > T_high;

imshowpair(A,mask_high)

That's good, but I'd like for the mask to include more of the peak, so I try a lower threshold.

T_low = 0.65;

mask_low = A > T_low;

imshowpair(A,mask_low)

Now the mask includes more of the peak of interest, but it also includes some of the smaller peaks that I don't want.

The double thresholding idea is to use the output from the higher threshold to select all the touching regions in the output from the lower threshold. Morphological reconstruction can do exactly that. The Image Processing Toolbox function to use is imreconstruct.

mask_double = imreconstruct(mask_high,mask_low);

imshowpair(A,mask_double)

Now the mask includes more of the high peak, without also including portions of the lower peaks.

Let's try a slightly less contrived example that uses double thresholding on the gradient magnitude of a grayscale image.

B = imread("coins.png");

imshow(B)

Bg = mat2gray(imgradient(B));

imshow(Bg)

Can we just get the enclosing edges for each coin using a simple threshold?

T_high = 0.6;

Bg_mask_high = Bg > T_high;

imshow(Bg_mask_high)

Try a lower threshold:

T_low = 0.3;

Bg_mask_low = Bg > T_low;

imshow(Bg_mask_low)

Now use a double threshold:

Bg_mask_double = imreconstruct(Bg_mask_high,Bg_mask_low);

imshow(Bg_mask_double)

y = [1 2 0.6 2 3 4 5 4 3 2 0.4 0.3 0.4 0.4 0.6 0.4 0.5 0.4 0.6 0.5 ...

0.3 0.5 0.6 0.2 1 2 3 0.6 1 2 3 2 3 2 5 4];

plot(y)

grid on

The problem was to extract low portions below 0.5, but using some kind of loose threshold that allowed portions of the clump to be a bit higher than 0.6. The double threshold does the trick nicely.

T_low = 0.5;

mask_low = y < 0.5;

T_high = 0.7;

mask_high = y < T_high;

plot(y)

hold on

plot(T_low*mask_low, LineWidth = 3, Color = "red")

plot(T_high*mask_high, LineWidth = 3, Color = "green")

hold off

Compute the double threshold using imreconstruct:

mask_double = imreconstruct(mask_low,mask_high);

plot(y)

hold on

plot(T_low*mask_double, LineWidth = 3, Color = "red")

hold off

Today's post departs from image processing, my usual subject. Instead, I want to tell you about something that I just put up on GitHub and the MATLAB File Exchange.Earlier this year, I was debugging... read more >>

]]>Today's post departs from image processing, my usual subject. Instead, I want to tell you about something that I just put up on GitHub and the MATLAB File Exchange.

Earlier this year, I was debugging some code, and I wanted to know what was happening inside a certain function while the function was running. Typically, the way to do this is to modify the code by adding something like a call to disp:

disp("Made it to this line, and the value of n is " + n)

But modifying the code is often undesirable. You have to remember to restore the code to its original state when you are done, and you might be working with code that is stored in a read-only location.

So, I started tinkering with ways to achieve a similar effect while leaving the code untouched, and I hit on the notion of using conditional breakpoints. When you set a conditional breakpoint on a code line in a file, then code execution stops at that line if an expression that you specify evaluates to true (1). Here is a screenshot of setting a conditional breakpoint using the MATLAB Editor. Here, I'm creating a breakpoint that will stop execution at line 2 of fib.m if the value of n is 4.

Well, we can use this to replace our use of disp if we can create an expression that:

- Always returns false (0)
- Has a side effect of printing information to the command window

That's the basis for my Code Trace for MATLAB, which is available on github.com/mathworks as well as on MATLAB File Exchange. Let me illustrate it using this function fib, which is a recursive implementation of a Fibonacci number generator.

function out = fib(n)

if n == 0

out = 0;

elseif n == 1

out = 1;

else

out = fib(n-1) + fib(n-2);

end

For the simplest example, I'll add a code trace by specifying the function name and the line number.

addCodeTrace("fib",2)

Now, a "trace" message will get displayed in the Command Window every time line 2 is about to be executed. To demonstrate, I'll call fib:

result = fib(3)

In each trace message, you see the function that is called, fib, and the line number that is about to be executed, 2. The level of indentation varies with the depth of the call stack, which helps us to visualize the recursive nature of this function.

For my next example, I'll display the value of n each time line 2 is about to executed.

clearCodeTraces

addCodeTrace("fib",2,Expression = "n")

result = fib(3)

I can also display a label. For my last example, I'll create several code traces that collectively show several details about the function fib as it executes.

clearCodeTraces

addCodeTrace("fib",2,Label = "function entry", Expression = "n")

addCodeTrace("fib",7,Label = "recursive function call")

addCodeTrace("fib",3,Label = "function exit, returning 0")

addCodeTrace("fib",5,Label = "function exit, returning 1")

result = fib(4)

Clean up.

clearCodeTraces

I'm interested to see whether people find this useful. Post a comment here if you have thoughts. If you encounter a problem or have a suggestion, feel free to create an issue on the GitHub repository.

]]>color segmentation, opening by reconstruction, boundary tracing, polyshapesEarlier this summer, Cleve sent me this picture of a puzzle.He asked for tips on turning the individual "cheese slices" into... read more >>

]]>color segmentation, opening by reconstruction, boundary tracing, polyshapes

Earlier this summer, Cleve sent me this picture of a puzzle.

He asked for tips on turning the individual "cheese slices" into patch objects for plotting and manipulation, perhaps something like this:

In my post today, I'll elaborate on this step.

The first thing that occurred to me was to use color to segment the puzzle pieces from the rest of the image. I used the Color Thresholder app to determine some threshold values in CIELAB space.

rgb = imread("Cheese_puzzle.png");

lab = rgb2lab(rgb);

% Threshold values chosen with the help of colorThresholder.

[L,a,b] = imsplit(lab);

mask = ((30 <= L) & (L <= 98)) & ...

((-20 <= a) & (a <= 16)) & ...

((26 <= b) & (b <= 83));

imshow(mask)

Next, we need to get rid of the extra unwanted foreground pixels. I'll use a method that the mathematical morphology folks call opening by reconstruction. The first step is to erode the image in a way that eliminates all the unwanted pixels, while maintaining at least a portion of all the objects we want to keep. An erosion by a vertical line will get rid of the thin horizontal lines, and an erosion by a horizontal line will get rid of the thin vertical lines.

mask2 = imerode(mask,ones(21,1));

imshow(mask2)

mask3 = imerode(mask2,ones(1,21));

imshow(mask3)

In the image above, the extraneous foreground pixels are gone, but the puzzle pieces have been shrunk. I'll use morphological reconstruction to recover the entire puzzle pieces.

mask4 = imreconstruct(mask3,mask);

imshow(mask4)

Great. Now, how do we turn these puzzle pieces into some kind of easily plottable representations? Cleve had suggested patch objects, but it is complicated to create patch objects containing holes. I suggested using polyshape instead. A polyshape is an object that can represent shapes that are composed of polygons. You can create a polyshape from a collection of polygons that bound individual regions and holes, and the Image Processing Toolbox function bwboundaries can produce just such a list of polygons.

b = bwboundaries(mask4)

The polyshape function can take a collection of bounding polygons and automatically figure out which polygons bound regions and which bound holes, but the form of the polyshape input arguments is a bit different from what bwboundaries produces. Here is some code to convert the bwboundaries output into something that polyshape can handle.

for k = 1:length(b)

X{k} = b{k}(:,2);

Y{k} = b{k}(:,1);

end

ps = polyshape(X,Y)

Well, that's a bit messy. In my image processing work, I often see this warning message from polyshape. Usually it is because I'm passing a bunch of colinear vertices to polyshape. They are colinear because the vertices are on the image pixel grid. I generally ignore the warning.

There's another issue, though. Why does the output of polyshape say it has 6 regions instead of 5? I had to search a bit for the reason. We can see it by zooming into one of the corners of the lower right puzzle piece and superimposing the boundary traced by bwboundaries.

figure

imshow(mask)

hold on

plot(b{5}(:,2),b{5}(:,1),'b')

hold off

xlim([1120 1140])

ylim([990 1010])

There's a spot where the traced boundary is self-intersecting, and that causes polyshape to treat that tiny triangle at the bottom as a separate region. We can do a little area-based processing of the polyshape to get rid of that triangle.

Split the polyshape into separate regions:

ps_regions = regions(ps)

Find the area of each region:

region_areas = area(ps_regions)

And get rid of the tiny region:

ps_regions(region_areas < 1) = []

Now we can plot our puzzle pieces. If you pass an array of polyshapes to plot, it will color each one separately.

plot(ps_regions)

axis ij

axis equal

I love polyshapes. They come with a rich collection of useful functions (see "Object Functions" on the polyshape reference page.) You can modify them, take them apart, join them using Boolean operations, measure them, query point locations, etc.

For example, the polybuffer function can expand (or shrink) a polyshape based on the idea of creating a "buffer zone" around the shape. To illustrate, let's take one of the puzzle pieces and give it a 50-pixel buffer zone.

ps5 = ps_regions(5);

ps5_50 = polybuffer(ps5,50);

figure

plot(ps5_50)

axis ij

axis equal

hold on

plot(ps5)

hold off

Are polyshapes useful in your work? I'd like to hear about it. Please leave a comment.

]]>Today I want to show you a handy little function, imsplit, for splitting color images into their components.Also, for no reason other than WOW, I want to show one of the images just released by the... read more >>

]]>Today I want to show you a handy little function, imsplit, for splitting color images into their components.

Also, for no reason other than WOW, I want to show one of the images just released by the James Webb Space Telescope program (home page, flickr galleries).

url = "https://live.staticflickr.com/65535/52211883534_7fe30b9955_o_d.png";

rgb = imread(url);

imshow(rgb)

text(size(rgb,2),size(rgb,1),"Image credits: NASA, ESA, CSA, and STScl",...

VerticalAlignment = "top",...

HorizontalAlignment = "right")

Anyway, back to imsplit. Long-time readers have seen me use code like this to split up color images (sometimes RGB, sometimes $ L^* a^* b^* $) into component images:

L = lab(:,:,1);

a = lab(:,:,2);

b = lab(:,:,3);

See, for example, the posts about two-dimensional histograms (23-Dec-2010) and segmenting images in $ (a^*,b^*) $ space (04-Feb-2011).

That code is not hard to write, but now I've become fond of using imsplit instead, just because it is a bit more compact, without sacrificing understandability. The function imsplit, introduced in release R2018b, is used this way:

url = "https://blogs.mathworks.com/images/steve/2010/mms.jpg";

candy = imread(url);

lab = rgb2lab(candy);

[L,a,b] = imsplit(lab);

tiledlayout("flow")

nexttile

imshow(candy)

nexttile

imshow(L,[0 100])

title("L")

nexttile

imshow(a,[-90 90])

title("a")

nexttile

imshow(b,[-90 90])

title("b")

That's it, really. There's nothing more to it. Give it a try the next time you need to do this; it'll save you 2.7183 seconds.

Before I go, let's look at that Webb telescope image one more time, just because.

[R,G,B] = imsplit(rgb);

tiledlayout("flow")

nexttile

imshow(rgb)

nexttile

imshow(R)

title("R")

nexttile

imshow(G)

title("G")

nexttile

imshow(B)

title("B")

Today's post shows how to use ismember to conveniently locate all of the highest pixel values in an image. For example, in the following Hough transform image, what are the 20 highest values, and... read more >>

]]>Today's post shows how to use ismember to conveniently locate all of the highest pixel values in an image. For example, in the following Hough transform image, what are the 20 highest values, and where are all the pixels that have those values?

I got the idea from some example code written by Image Processing Toolbox writer Megan Mancuso. Megan says she got the idea to use ismember this way from Ahmed's answer on MATLAB Answers.

Here's where the Hough transform image above comes from.

A = imread("gantrycrane.png");

imshow(A)

bw = edge(rgb2gray(A),"canny");

imshow(bw)

[H,T,R] = hough(bw);

imshow(H,[],XData=T,YData=R)

daspect auto

axis on

Now let's sort to figure out the highest 20 values of H.

sorted_H_values = sort(H(:),"descend");

highest_H_values = sorted_H_values(1:20)

Finally, I'll use ismember to compute a binary image showing all the elements of H that have one of those 20 highest values. (And I'll use xlim and ylim to zoom in so that we can see a few of those elements clearly.)

mask_20_highest_values = ismember(H,highest_H_values);

imshow(mask_20_highest_values,XData=T,YData=R)

daspect auto

axis on

xlim([-10 60])

ylim([141 452])

Any element-wise logical function, such as ismember, can be used in this way to produce a binary image. So add this to your bag of MATLAB tricks!

I'm still playing around with RGB gamut calculations in $ L^* a^* b^* $ space. (See my last post on this topic, "Visualizing out-of-gamut colors in a Lab curve.") Today's post features some new... read more >>

]]>I'm still playing around with RGB gamut calculations in $ L^* a^* b^* $ space. (See my last post on this topic, "Visualizing out-of-gamut colors in a Lab curve.") Today's post features some new gamut-related visualizations, plus some computational tricks involving gamut boundaries and rays in $ L^* a^* b^* $ space.

First, here is another way to communicate the idea that the in-gamut area in the $ (a^*,b^*) $ plane varies with $ L^* $ (perceptual lightness). For 9 values of $ L^* $, (10, 20, ..., 90), I'll compute a 2-D $ (a^*,b^*) $ gamut mask by brute-forcing it. Then, I'll use overlaid contour plots to show the variation in gamut boundaries.

a = -110:.1:110;

b = -110:.1:110;

L = 10:20:90;

[aa,bb,LL] = meshgrid(a,b,L);

hold on

for k = 1:length(L)

rgb = lab2rgb(cat(3,LL(:,:,k),aa(:,:,k),bb(:,:,k)));

mask = all((0 <= rgb) & (rgb <= 1),3) * 2 - 1 + L(k);

contour(a,b,mask,[L(k) L(k)],LineColor=[.8 .8 .8],LineWidth=1.5,ShowText=true,...

LabelSpacing=288)

end

hold off

axis equal

grid on

box on

xlabel("a*")

ylabel("b*")

title("Gamut boundary in the (a,b) plane for several values of L*")

Here's another visualization concept. People often show colors on the $ (a^*,b^*) $ plane, to give an idea of the meaning of $ a^* $ and $ b^* $, but that doesn't communicate very well the idea that there are usually multiple colors, corresponding to various $ L^* $ values, at any one $ (a^*,b^*) $ location. Below, I show both the brightest in-gamut color and the darkest in-gamut color at each $ (a^*,b^*) $ location.

a = -110:110;

b = -110:110;

[aa,bb] = meshgrid(a,b);

L_max = zeros(size(aa));

L_min = zeros(size(aa));

for p = 1:size(aa,1)

for q = 1:size(bb,1)

[L_min(p,q),L_max(p,q)] = Lrange(aa(p,q),bb(p,q));

end

end

rgb = lab2rgb(cat(3,L_max,aa,bb));

figure

tiledlayout(1,2)

nexttile

imshow(rgb,XData=a([1 end]),YData=b([1 end]))

axis xy

axis on

xlabel a

ylabel b

title("Brightest in-gamut color")

rgb_min = lab2rgb(cat(3,L_min,aa,bb));

nexttile

imshow(rgb_min,XData=a([1 end]),YData=b([1 end]))

axis on

axis xy

xlabel a

ylabel b

title("Darkest in-gamut color")

Next, I find myself sometimes wanting to draw a ray in $ L^* a^* b^* $ space and find the gamut boundary location along that ray. To that end, I wrote a simple utility function (findNonzeroBoundary below) that performs a binary search to find where a function goes from positive to 0. Then, I wrote some anonymous functions to find the desired gamut boundary point. Specifically, I was interested in this question: For a given $ L^* $ value and a given $ (a^*,b^*) $ plane angle, h, what is the in-gamut color with the maximum chroma, c, or distance from $ (0,0) $ in the $ (a^*,b^*) $ plane?

Fair warning: the code below gets tricky with anonymous functions. You might hate it. If so, I totally understand, and I hope you'll forgive me. :-)

I'll start by creating an anonymous function that converts from $ L^* c h $ to $ L^* a^* b^* $:

lch2lab = @(lch) [lch(1) lch(2)*cosd(lch(3)) lch(2)*sind(lch(3))];

Next, here is an anonymous function that returns whether or not a particular $ L^* a^* b^* $ point is in gamut.

inGamut = @(lab) all(0 <= lab2rgb(lab),2) & all(lab2rgb(lab) <= 1,2);

Finally, a third anonymous function uses findNonzeroBoundary to find the gamut boundary point that I'm interested in.

maxChromaAtLh = @(L,h) findNonzeroBoundary(@(c) inGamut(lch2lab([L c h])), 0, 200);

Let's exercise this function to find a high chroma dark color at $ h=0^{\circ} $.

L = 35;

h = 0;

c = maxChromaAtLh(L,h)

And here's what that color looks like.

rgb_out = lab2rgb(lch2lab([L c h]));

figure

colorSwatches(rgb_out)

axis equal

axis off

(The function colorSwatches is from Digital Image Processing Using MATLAB, and it is included in the MATLAB code files for the book.)

What happens when we try to find a high chroma color, at the same hue angle, that is bright instead of dark?

L = 90;

h = 0;

c = maxChromaAtLh(L,h)

You can can see that the maximum c value is much lower for the higher value of $ L^* $. What does it look like?

rgb_out = lab2rgb(lch2lab([L c h]));

figure

colorSwatches(rgb_out)

axis equal

axis off

When I was doing these experiments to prepare for this blog post, my original intent was to just show examples for a couple of different values of h and $ L^* $. But I couldn't stop! It was too much fun, and I kept trying different values.

After about 15 minutes or so, I decided it would be best to write some simple loops to generate a relatively large number of $ (L^*,h) $ combinations to look at. Here's the code to generate the high c colors for a variety of $ (L^*,h) $ combinations.

dh = 30;

h = -180:dh:150;

L = 35:15:95;

dL = 15;

rgb = zeros(length(h),length(L),3);

for q = 1:length(L)

for k = 1:length(h)

c = maxChromaAtLh(L(q),h(k));

rgb(k,q,:) = lab2rgb(lch2lab([L(q) c h(k)]));

end

end

And here's the code to view all those colors on a grid, with labeled h and $ L^* $ axes.

rgb2 = reshape(fliplr(rgb),[],3);

p = colorSwatches(rgb2,[length(L) length(h)]);

p.XData = (p.XData - 0.5) * (dh/1.5) + h(1);

p.YData = (p.YData - 0.5) * (dL/1.5) + L(1);

xticks(h);

xlabel("h")

yticks(L)

ylabel("L*")

grid on

title("Highest chroma (most saturated) colors for different L* and h values")

function [L_min,L_max] = Lrange(a,b)

arguments

a (1,1) {mustBeFloat}

b (1,1) {mustBeFloat}

end

L = 0:0.01:100;

a = a*ones(size(L));

b = b*ones(size(L));

lab = [L ; a ; b]';

rgb = lab2rgb(lab);

gamut_mask = all((0 <= rgb) & (rgb <= 1),2);

j = find(gamut_mask,1,'first');

k = find(gamut_mask,1,'last');

if isempty(j)

L_min = NaN;

L_max = NaN;

else

L_min = L(j);

L_max = L(k);

end

end

function x = findNonzeroBoundary(f,x1,x2,abstol)

arguments

f (1,1) function_handle

x1 (1,1) {mustBeFloat}

x2 (1,1) {mustBeFloat}

abstol (1,1) {mustBeFloat} = 1e-4

end

if (f(x1) == 0) || (f(x2) ~= 0)

error("Function must be nonzero at initial starting point and zero at initial ending point.")

end

xm = mean([x1 x2]);

if abs(xm - x1) / max(abs(xm),abs(x1)) <= abstol

x = x1;

elseif (f(xm) == 0)

x = findNonzeroBoundary(f,x1,xm);

else

x = findNonzeroBoundary(f,xm,x2);

end

end

[TLDR: To create a fast GIF animation using imwrite, use a DelayTime of 0.02 instead of 0.]In my 28-Feb-2022 post, I showed some code that created an animated GIF file using a DelayTime of 0.im =... read more >>

]]>[TLDR: To create a fast GIF animation using imwrite, use a DelayTime of 0.02 instead of 0.]

In my 28-Feb-2022 post, I showed some code that created an animated GIF file using a DelayTime of 0.

im = rand(750,750)>0.8;

for k=1:500

t = conv2(im(:,:,k),[2,2,2;2,1,2;2,2,2],'same');

im(:,:,1,k+1) = (t > 4) & (t < 8);

end

imwrite(~im,'life.gif','DelayTime',0,'LoopCount',1)

In response, reader Mikhail commented that one shouldn't use a DelayTime of 0 for the purpose of creating fast GIF animations. He pointed me to an blog post by Ben Phelps called "The Fastest GIF Does Not Exist." (Note: If you visit this site, you may get a security warning in your browser. Apparently, the site's certificate expired a few days ago. Visit at your own risk. I have attempted to contact the author to let them know.)

Many years ago (about 25!), I used to work on image file format support in MATLAB, so I usually start investigating such questions by looking at the specification for the format in question. The relevant spec is GIF89a, and here's the key paragraph:

Notice that this language doesn't specify exactly what is supposed to happen if the Delay Time is 0.

One thing I've learned by experience is that the actual behavior of commonly used software applications sometimes varies from the spec, and implementers must pay attention to that. Ben has looked at the publicly available source code for several browsers and other applications and discovered a Delay Time of 0 or 1 (1/100 s) is generally not honored. Instead, for various practical and historical reasons, a very low Delay Time like 0 or 1 is often bumped up arbitrarily to a value of 10 (1/10 s).

Therefore, setting Delay Time to 0 is not going to give you the fastest possible animation! I tested this in Safari, the browser I use on my Mac, and this does appear to be true.

This GIF was generated using DelayTime of 0:

This one was generated using a DelayTime of 0.02 (imwrite uses seconds instead of hundredths of seconds for this parameter):

Indeed, the second version does animate much more quickly.

One implication is that we may need to change the imwrite documentation, which specifically says to use DelayTime of 0 to get the fastest animation:

I actually think it may be time to change the imwrite interface. Perhaps we could give this parameter a name, such as FrameRate, (in frames per second) and document a maximum value of 50 (which corresponds to a DelayTime of 0.02 s (or 2 hundredths of a second).

I will make both suggestions to the development team responsible for imwrite.

]]>