Of all the MATLAB or MathWorks toolbox functions that have been mentioned in this blog since 2006, imshow and imread have appeared the most. As in last week's post, they often appear together, right at the beginning: ... read more >>

]]>When I say they have "appeared the most," I am counting each function once for each blog post in which it has been used or mentioned.

I'll send a MATLAB t-shirt to the first two people who can name another function in the top five. Make your guesses in the comments. One guess per person, please.

]]>Note: this blog post describes an image display feature that is new in R2019b. You can download this release from https://mathworks.com/downloads.... read more >>

]]>*Note: this blog post describes an image display feature that is new in R2019b. You can download this release from* https://mathworks.com/downloads.

This is my new favorite image:

But that's just a 512x512 version of it. Click here for the full resolution (2048x2048) version. *Image credit: DeFacto [CC BY-SA 4.0 (https://creativecommons.org/licenses/by-sa/4.0)]*

This image caught my eye because I'm always looking for a good sample image to demonstrate *aliasing*. Also, I wanted to write a blog post about a new MATLAB image display feature in R2019b: bilinear interpolation and antialiasing.

Before I get into image display details, though, let's zoom into the center of the full-resolution image to see why it interests me.

```
url = 'https://blogs.mathworks.com/steve/files/Decimal_Clock_face_by_Pierre_Daniel_Destigny_1798-1805.jpg';
A = imread(url);
A_cropped = A((-250:250) + 1024, (-250:250) + 1024, :);
imshow(A_cropped)
```

Next, I want to compute a pixel profile along a path that shows the high-frequency spatial variation. The code below shows the path in blue.

hold on x = [120 80]; y = [106 188]; plot(x,y,'b','LineWidth',4) hold off

And here's the pixel profile (converted to gray-scale) along that path.

close all c = improfile(A_cropped,x,y); c = rgb2gray(squeeze(c)/255); h = plot(c(:,1)); xlabel('Distance along the profile path') ylabel('Gray-scale pixel value') datatip(h,13,0.9614); datatip(h,25,0.9271); datatip(h,60,0.9282); datatip(h,72,0.9133);

I have placed a few datatips to show that the spatial periodicity along the path is about 4 pixel units. This is *very close* to the highest spatial frequency that can be represented with this pixel spacing. In other words, it is very close to what I would call *critically sampled*.

Why does this matter? Well, it means that if you shrink this image down just by using nearest-neighbor interpolation, which is the way MATLAB has always displayed images previously, you'll see some strange artifacts. Below, I'll use `imshow` to display the image at 25% magnification.

```
imshow(A_cropped,'InitialMagnification',25)
```

With R2019b, though, you can specify that you want an image to be displayed using bilinear interpolation. And, as explained in the documentation, when bilinear interpolation is specified, MATLAB also automatically applies an antialiasing technique.

imshow(A_cropped,'InitialMagnification',25,'Interpolation','bilinear')

You can still see some spatial artifacts in this result, but they are significantly reduced.

Now let's examine the effect of specifying bilinear interpolation when you magnify an image instead of shrinking it. First, I'll crop out an even smaller piece of the original image.

B_cropped = A_cropped(50:150,115:215,:); imshow(B_cropped)

Next, let's display it with 500% magnification, with and without the bilinear interpolation option.

imshow(B_cropped,'InitialMagnification',500) title('500% magnification using nearest-neighbor interpolation')

imshow(B_cropped,'InitialMagnification',500,'Interpolation','bilinear') title('500% magnification using bilinear interpolation')

You can see the improved quality in the version using bilinear interpolation.

Note that the default is still nearest-neighbor interpolation. This was a choice that we discussed at some length. Ultimately, we decided to leave the default behavior unchanged, for two reasons:

- to avoid incompatibilities
- for some image arrays, such as those containining segmentation labels or other kinds of categories, bilinear interpolation would not be appropriate

If you already have MATLAB R2019b, give this image display option a try and let me know what you think in the comments.

Get
the MATLAB code

Published with MATLAB® R2019b

In my second year of grad school (1987-88), my thesis advisor asked me to develop some PC-based digital signal processing (DSP) tools to be used for computer-based lab assignments. (At the time, PC-based engineering labs were very new, laptops weren't a thing yet, and Student Edition of MATLAB was... read more >>

]]>I think about this experience whenever I interview a candidate for a software development job. The MathWorks has many products in specialized areas of engineering, science, and mathematics, so we often hire domain experts in these areas to develop software. Our biggest challenge in evaluating a candidate is to answer this critical question:

Do they have the tool builder's gene?

Here's what I mean: Many candidates can describe something that they have created that was useful to someone else. But only a few find the tool-building activity itself to be intensely rewarding. Only a few enjoy extending the tool so that it's useful not only to themselves and maybe a few others in the lab, but to a wide audience. Only a few are motivated to refine and extend a tool after it has solved the problem for which it was originally created. Those with the tool builder's gene get their energy from knowing they've helped many people do their own jobs better.

From my sophomore year in college through getting my Ph.D., I focused single-mindedly on becoming an engineering professor, and I did become one for a few years. But I was fortunate to realize quickly that academia ultimately wasn't the best place for me, and that I would be much happier at a place like The MathWorks. The reason is that I have the gene.

How about you? Are you interested in working at MathWorks? Then please visit our careers page. And here's a hint: the Engineering Development Group is a great way to get started on a technical career at MathWorks.

PS. Hey, the blog comment section has a new look and features, including better formatting, threaded replies, and notifications. Try it out: Share with us a story about your favorite tool you built for others to use.

]]>Back in the summer I had another chance to use the Color Thresholder, a very nice app that's in the Image Processing Toolbox. I happened to come across a question on MATLAB Answers - someone was looking for a way to segment the yellow region in this image:... read more >>

]]>Back in the summer I had another chance to use the Color Thresholder, a very nice app that's in the Image Processing Toolbox. I happened to come across a question on MATLAB Answers - someone was looking for a way to segment the yellow region in this image:

I proposed using the Color Thresholder and the L*a*b* color space, and I showed this screen shot to demonstrate how to apply it to this problem:

I first wrote about the Color Thresholder back in 2014, shortly after the app first appeared ("Color Thresholder App in R2014a"). In that post, I played around with picture of candy lying on my desk:

Besides figuring out how to segment all the different candy colors, I also discovered how to characterize the color of my desk:

The Color Thresholder is a very useful tool. Give it a look!

Get
the MATLAB code

Published with MATLAB® R2019b

I recently noticed a change in the way we write some of our product release notes, and I wanted to mention it to you.... read more >>

]]>I recently noticed a change in the way we write some of our product release notes, and I wanted to mention it to you.

In my quarter century at MathWorks doing toolbox and MATLAB development, there have been a few areas of focus that have been remarkably consistent over that entire time. One of those areas is performance. Specifically, computation speed.

If you have used MATLAB more than five years, it is likely that something you use in MATLAB a lot has been completely reimplemented to make it go faster in our ever-evolving computational environments.

Maybe it was new algorithms, like image resizing or Gaussian filtering. Maybe the memory access patterns were modified to exploit changing memory cache architectures, like image resizing (again), transposition (and `permute`), `conv2`, and even the seemingly straightforward `sum` function.

Possibly the functions you rely upon were modified to adapt to new core libraries, such as LAPACK or FFTW.

Many, many, many functions and operators were completely overhauled when multicore computers became common. Then they were modified again to exploit extended processor instruction sets for instruction word parallelism.

Finally, the very foundations for MATLAB language execution were completely overhauled in 2015 to make everything go faster. Since then, the MATLAB execution engine continues to be refined with almost every release to add new types of optimizations.

The curious thing about all this effort, over so many years, is how ... well ... **vague** we typically have been in describing performance improvements in our release notes.

For example, here is a snippet from the R2018b Release Notes for the Image Processing Toolbox:

Like I said: it's vague.

It was never our intent to be obscure. It's just that performance measurements are almost always challenging to report with accuracy and precision, and the experiences of individual users will almost always vary, sometimes considerably. Part of our company culture here is that we are allergic to making statements that could be perceived as inaccurate. I think that's what has been behind the history of vague statements about performance improvements in release notes. (OK, I should state this explicitly: this is my personal opinion, and not a statement of what company policy is or has been.)

Well, things are starting to change. Our documentation writers now have a new standard to follow when writing release notes about performance. Here is a sample from R2019b, which was released last month:

The release note describes what operation has been improved, how it was timed, what the times were for specific releases, and details about the computer used to measure the performance.

Look for more performance changes to be reported with this level of detail in the future. I think this is a great improvement!

Get
the MATLAB code

Published with MATLAB® R2019b

I thought this was worth sharing. ... read more >>

]]>From NASA's hubblesite.org:

This new Hubble Space Telescope view of Jupiter, taken on June 27, 2019, reveals the giant planet's trademark Great Red Spot, and a more intense color palette in the clouds swirling in Jupiter's turbulent atmosphere than seen in previous years. The colors, and their changes, provide important clues to ongoing processes in Jupiter's atmosphere.

You should definitely click the image to see the full resolution.

Image credits: NASA, ESA, A. Simon (Goddard Space Flight Center), and M.H. Wong (University of California, Berkeley)

]]>I recently saw some code that transformed the RGB pixel values of an image into a Px3 matrix, such that each row contained the red, green, and blue color components of a single pixel. Today I want to show that code fragment, explain it, and then demonstrate what I think... read more >>

]]>I recently saw some code that transformed the RGB pixel values of an image into a Px3 matrix, such that each row contained the red, green, and blue color components of a single pixel. Today I want to show that code fragment, explain it, and then demonstrate what I think is the fastest possible to perform that transformation in MATLAB.

I tend to use the peppers sample image a lot, so I'll switch things up and use NGC6543 instead:

```
RGB = imread('ngc6543a.jpg');
imshow(RGB)
```

Historical note: This was the first "truecolor" sample image to ship with MATLAB.

`RGB` is a three-dimensional array:

size(RGB)

ans = 650 600 3

The third dimension, which has length three, contains the red, green, and blue component images. Compare, for example, the red and green component images:

imshow(RGB(:,:,1))

imshow(RGB(:,:,2))

Here is the code that I saw for transforming this array into a Px3 matrix, where each row contains the three color component values of one pixel:

red = RGB(:,:,1); green = RGB(:,:,2); blue = RGB(:,:,3); A = [red(:) green(:) blue(:)];

The first five pixels are along the uppermost left edge of the image, and they are all black:

A(1:5,:)

ans = 5×3 uint8 matrix 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

The first three lines of code did not give me pause:

red = RGB(:,:,1); green = RGB(:,:,2); blue = RGB(:,:,3);

I have written code just like that approximately 2.7183 bazillion times. It was the next line that made me stop and scratch my head for a moment:

A = [red(:) green(:) blue(:)];

Hmm. The array is being split apart into three components, which are then immediately put back together into a matrix. That suggests ... `reshape`, maybe?

Yes, as it turns out, this can all be done with a call to `reshape`. And, there is a big advantage to doing so, as I'll explain in a moment. First, the new code:

B = reshape(RGB,[],3);

That line of code says, "MATLAB, reshape the array `RGB` so that it is Px3, and will you please just go ahead and figure out what P should be to make it work?"

You see in a call to `reshape`, the number of elements must not change. Here are the number of elements of the array `RGB`:

numel(RGB)

ans = 1170000

To reshape that successfully into a matrix with 3 columns, we can compute how many rows there must be:

P = numel(RGB) / 3

P = 390000

The empty matrix as the second argument to `reshape` signals to MATLAB to compute that size automatically. That works out just fine:

size(B)

ans = 390000 3

So, each column of `B` has 390,000 elements. The first column of `B` contains the first 390,000 elements (in memory order) of `RGB`. The second column of `B` contains the second 390,000 elements of `RGB`, and similarly for the third column. Because of the way MATLAB arrangements elements in memory, these columns are exactly the elements from the red component, the green component, and the blue component of `RGB`.

The big advantage of doing it this way is that MATLAB makes no memory copy of data when calling `reshape`. Because of this optimization, calls to `reshape` are almost instantaneous.

Let `reshape` be your friend.

Get
the MATLAB code

Published with MATLAB® R2019a

For today's post, I'd like to welcome guest author Kushagr Gupta, a developer on the Image Processing Toolbox team. -Steve... read more >>

]]>*For today's post, I'd like to welcome guest author Kushagr Gupta, a developer on the Image Processing Toolbox team. -Steve*

The Image Processing Toolbox function `ordfilt2` can be used to perform various types of 2-D statistical filtering, including minimum or maximum filtering and median filtering. A straightforward implementation of this function is expected to scale with $O(r^2)$, where $r$ is the width of the window (assuming a square window). To compute each output pixel, we have to select the $k$-th largest value from the set of $r^2$ input image pixels, which can be done in $O(r^2)$ time with an algorithm such as quick-select.

The function `ordfilt2`, however, does better than that; it scales with $O\left(r\right)$complexity. That's because `ordfilt2` uses a sliding window-based histogram approach to filter the data. When the algorithm moves from one input image pixel to the next, only part of the histogram needs to be updated. (Reference: Huang, T.S., Yang, G.J., and Yang, G.Y., "A fast two-dimensional median filtering algorithm," *IEEE transactions on Acoustics, Speech and Signal Processing*, Vol ASSP 27, No. 1, February 1979.)

As shown in the figure above, r additions and r subtractions need to be carried out for updating the histogram. Computing the output pixel from the histogram doesn't depend on the window size, and so the overall complexity is $O(r)$.

We recently received an interesting tech support case related to performance of `ordfilt2` for uint16 random input data. One of our customers reported the time taken to filter an image with a smaller window size was **higher** than the time taken to process the same image with a relatively larger window size.

Was this just a noisy observation? Was this reproducible? Was this a bug in our implementation?

We investigated?

The customer was generous enough to give us reproduction code. And, the behavior was indeed reproducible! (Many, many times!). What was the explanation?

Here are the reproduction steps:

Im = im2uint16(rand(1e3)); n = [7 9 11 13 15 21 25 31 35 41 49 51 57 65 73]; t = zeros(length(n),1); for id=1:length(n) t(id) = timeit(@() ordfilt2(Im,ceil(n(id)^2*0.5),... true(n(id)))); end plot(n,t) grid on title('ordfilt2 performance curve (uint16)') xlabel('window size') ylabel('Time (s)')

A puzzling aspect of this case was that the performance behavior was observed only for random input data and not for real-world images. Here is a comparison:

Im1 = im2uint16(imresize(imread('coins.png'),[500 500])); Im2 = im2uint16(rand(500)); n = [11 15 21 27 31 35 41 49 51]; treal = zeros(length(n),1); trand = zeros(length(n),1); for id=1:length(n) treal(id) = timeit(@() ordfilt2(Im1,ceil(n(id)^2*0.5),... true(n(id)))); trand(id) = timeit(@() ordfilt2(Im2,ceil(n(id)^2*0.5),... true(n(id)))); end plot(n,treal) grid on hold on plot(n,trand) legend('Real World Image','Random Image') title('ordfilt2 performance comparison') xlabel('window size') ylabel('Time (s)') hold off

We also observed that the the performance behavior was not reproducible for uint8 random images. Instead, the algorithm's performance scaled linearly with the window size, as expected.

Im = im2uint8(rand(1e3)); n = [7 9 11 13 15 21 25 31 35 41 49 51 57 65 73]; t = zeros(length(n),1); for id=1:length(n) t(id) = timeit(@() ordfilt2(Im,ceil(n(id)^2*0.5),... true(n(id)))); end plot(n,t) grid on title('ordfilt2 performance curve (uint8)') xlabel('window size') ylabel('Time (s)')

To understand what was going on, it became essential to dig deeper into the inner workings of the underlying algorithm.

As it turns out, while using the sliding window histogram approach, the algorithm also keeps track of the running sum of the histogram so that it does not need to be computed for each pixel in the row. It can be updated accordingly as the window slides by traversing the histogram in the right direction towards the element of interest. This running-sum computation makes the algorithm data dependent. Moreover, the worst case scenario for this data-dependent algorithm happens when the input image contains random data.

But why is the behavior prominently present for small window sizes? It can be argued that, as the window size increases, we draw greater number of samples from the underling distribution resulting in smaller variance. As the variance goes down, the code path becomes more predictable and the algorithm obeys linear time complexity. However, when the window size reduces, the data dependency of the algorithm causes the variance to increase (which is a lot more for random input) making the code path less predictable. This causes the behavior reported by the customer of smaller windows taking more time to process as compared to larger windows.

After some consideration, we have decided to leave the implementation of `ordfilt2` as it is. This decision was made based on the fact that the behavior is primarily visible for random input data, which is not usually what customers use `ordfilt2` for.

Get
the MATLAB code

Published with MATLAB® R2019a

My wife, Geri Eddins, has been making a lot of quilts lately. A few months ago, I printed out these 10 colors, showed them to her, and asked, "Can you make a wall hanging quilt with this color scheme? And make it look kind of mathematical?" (Whatever that means!)... read more >>

]]>My wife, Geri Eddins, has been making a lot of quilts lately. A few months ago, I printed out these 10 colors, showed them to her, and asked, "Can you make a wall hanging quilt with this color scheme? And make it look kind of mathematical?" (Whatever that means!)

These are ten of the colors from the default MATLAB colormap, called parula. I've posted about this colormap several times.

So, we looked at some pattern books and found something that seemed promising. She took the color samples and bought some matching fabrics, and then she did her sewing magic. Pretty soon, she had finished two wall hangings! I hung one up at home and one at my MathWorks office. Here is the final result:

I took this picture in my office, using my phone's camera, and lit with sunlight filtered through office windows. But the picture above isn't the first thing I saw after taking the picture. Instead, this is what I got:

That's quite a difference! The yellows are almost completely gone. I can only guess that the way I framed the picture really fooled the phone camera's automatic white balancing algorithm.

After puzzling over what to do with that picture, I remembered about a new-ish (R2017b) color balancing function in the Image Processing Toolbox called `chromadapt`.

The simplest way to use this function is with the call:

out = chromadapt(in,illuminant)

where `illuminant` is chosen as a pixel color from the input that **should** be gray or white, but isn't. Notice that, in this wall hanging, Geri has conveniently provided a white region: the background of the rectangular label at the bottom.

I used `impixelinfo` to grab a pixel color ([143 170 241]/255) from this region. Then I processed the original image as follows:

corrected = chromadapt(original,[143 170 241]/255);

The resulting image is the first one displayed above, the one with the clearly visible yellow hues.

Enjoy.

Get
the MATLAB code

Published with MATLAB® R2019a

In a complete coincidence, two different coworkers today sent me links to two different works of 3Blue1Brown on YouTube. I loved both of them, and I wanted to share them with you.... read more >>

]]>The first one demonstrates how Fourier series can be used to trace out two-dimensional shapes on a plane. I've seen animations like this before, but this one might be the most beautifully done.

*Pure Fourier series animations for 12 oddly satisfying minutes*

The second is an exploration of Hilbert curves, including a connection to image processing.

*Hilbert's Curve: Is infinite math useful?*

For some related material, see my 25-Jan-2012 post on generating Hilbert curves.

I hope you enjoy these videos.

]]>