Contents

Here's a GUI Using the Tool

Here I show the tool embedded in a GUI (graphical user interface). There's a slider that allows the user to choose a viewing plane through the x-axis of the 3-D volume.

[x,y,z,v] = flow();
volvisGUI(x,y,z,v)

Code Architecture

The code is split into two files, the actual GUI, in volvisGUI.m and the viewer that shows planes slicing through the volume, in volumeVisualization.m. Let's first look at a skeleton of the code for volumeVisualization.

dbtype volumeVisualizationPseudo

1     function s = volumeVisualizationPseudo(x,y,z,v)
2     %volumeVisualizationPseudo   Pseudo-code for volumeVisualization.
3     
4     %   Copyright 2007 The MathWorks, Inc.
5     
6     %% Initialize
7     initDisplay(x,y,z,v)
8     
9     %% Nested Functions
10        %% Add slice plane
11        function addSlicePlane(xLoc)
12        end
13        %% Delete Slice plane
14        function deleteLastSlicePlane()
15        end
16        %% Initialize DIsplay
17        function initDisplay(x,y,z,v)
18        end
19    
20    % Export structure
21    s.addSlicePlane = @addSlicePlane;
22    s.deleteLastSlicePlane = @deleteLastSlicePlane;
23    s.xMin = min(x(:));
24    s.xMax = max(x(:));
25    end

Initialization

In this code, you can see that the figure is initialized, and a structure is returned. This structure contains the minimum and maximum values possible for a plane along the x-axis as well as function handles to two functions, one to add a slice plane, and one to delete the last slice plane. The two function handles are the functionality I may then use if I want to build on top of this visualization.

Let me try out the engine now. First I create the visualization.

close all
s = volumeVisualization(x,y,z,v)

s = 
           addSlicePlane: @volumeVisualization/addSlicePlane
    deleteLastSlicePlane: @volumeVisualization/deleteLastSlicePlane
                    xMin: 0.1000
                    xMax: 9.9000

Next, let me add and remove some slices.

s.addSlicePlane(5)

s.addSlicePlane(3.2)

Now delete the last plane added.

s.deleteLastSlicePlane()

And again.

s.deleteLastSlicePlane()

And again. Or can I? I only added two slices and I am now up to deleting the third!

s.deleteLastSlicePlane()

The Engine Code

Whew! No seg-fault or other crash. Let's look at the actual code to see what's going on. I'll discuss it in segments.

dbtype volumeVisualization

1     function s = volumeVisualization(x,y,z,v)
2     %volumeVisualization Engine for choosing y-z planes to view.
3     %   s = volumeVisualization(X,Y,Z,V) returns a structure 
4     %   containing information about the visualization of volume
5     %   data described by X,Y,Z,V.  The fields are:
6     %          addSlicePlane -- function handle to add a slice
7     %                           plane at location x
8     %   deleteLastSlicePLane -- function handle to delete the
9     %                           last slice plane added
10    %                   xMin -- minimum x location for a plane
11    %                   xMax -- maximum x location for a plane
12    %
13    %      Example:
14    %      [x,y,z,v] = flow;
15    %      s = volumeVisualization(x,y,z,v);
16    %      s.addSlicePlane(3.7)
17    %      s.addSlicePlane(7.5)
18    %      pause
19    %      s.deleteLastSlicePlane()
20    %      pause
21    %      s.deleteLastSlicePlane()
22    
23    %   Copyright 2007 The MathWorks, Inc.
24    
25    %% Store handles to the various planes
26    %initialize handle to axis
27    %initialize handle to slice plane
28    hAxis = [];         
29    hSlicePlanes = [];  
30    
31    %% Create data for generic slice through yz-plane
32    [yd,zd] = meshgrid(linspace(min(y(:)),max(y(:)),100), ...
33        linspace(min(z(:)),max(z(:)),100));
34    
35    %% Plot the volume initially
36    initDisplay()
37    
38    %% Nested Functions
39        function addSlicePlane(xLoc)
40        %addSlicePlane   Add a slice plane xLoc.
41            xd            = xLoc*ones(size(yd));
42            newSlicePlane = slice(hAxis, x, y, z, v, xd, yd, zd);
43            hSlicePlanes   = [ hSlicePlanes, newSlicePlane ];
44            set(newSlicePlane,'FaceColor'      ,'interp',...
45                              'EdgeColor'      ,'none'  ,...
46                              'DiffuseStrength',.8       );
47        end
48    
49        function deleteLastSlicePlane()
50        %deleteLastSlicePlane Delete the last slice plane added.
51            if ~isempty(hSlicePlanes)
52                delete(hSlicePlanes(end));
53                hSlicePlanes = hSlicePlanes(1:end-1);
54            end
55        end
56    
57        function initDisplay()
58        %initDisplay  Initialize Display.
59    
60            % Draw back and bottom walls
61            if isempty(hAxis) || ~ishandle(hAxis)
62                hAxis = gca;
63                hold on;
64            end
65            hx = slice(hAxis, x, y, z, v, ...
66                max(x(:)),       [],       []) ;
67            hy = slice(hAxis, x, y, z, v, ...
68                [],       max(y(:)),       []) ;
69            hz = slice(hAxis, x, y, z, v, ...
70                [],              [],min(z(:))) ;
71    
72            % Make everything look nice
73            set([hx hy hz],'FaceColor','interp',...
74                'EdgeColor','none')
75            set(hAxis,'FontSize',18,'FontWeight','Bold');
76            xlabel('X');ylabel('Y');zlabel('Z')
77            daspect([1,1,1])
78            axis tight
79            box on
80            view(-38.5,16)
81            colormap (jet(128))
82        end
83    
84    s.addSlicePlane = @addSlicePlane;
85    s.deleteLastSlicePlane = @deleteLastSlicePlane;
86    s.xMin = min(x(:));
87    s.xMax = max(x(:));
88    
89    end

Help

We start off with the function declaration, followed by lines 2-21 for the help including an example.

Variable Initialization

We next initialize variables to hold the handle to the axes in which the data are plotted and another array for the handles to the slice planes. They start off empty. We also make an array of values for the y and z data for a generic slice. And then we initialize the plot, by calling the nested function initDisplay. You can see that we plot several slices along different directions and alter aspects of the plot such as colormap, aspects of the axes, etc.

Finally, we set up a structure to return the minimum and maximum values for the slices along the x-axis and function handles to two functions, for adding and deleting slices. These function handles refer to nested functions in the file volumeVisualization.m and I will discuss these functions in the next section.

The Main Work

The main work for this application is done after the initialization and is embodied in the two nested functions addSlicePlane and deleteLastSlicePlane.

When I call addSlicePlane, I insert a plane at a particular x-location. The function creates data for the slice plane, instantiates a slice, makes it pretty by altering some color and shading properties, and adds the handle to the latest slice to the list of slice handles. As simple as that.

Calling deleteLastSlicePlane again manages the list of slice planes by deleting the last element on the list from the plot, and removing it from the list of slice planes. If that list is empty, the function essentially does nothing, which is why it had the good graces to not error out on me earlier.

The GUI

I am now poised to use this M-file to create the initial GUI I showed. That GUI has a slider that maps to positions for planes along the x-axis. When a user selects a value, the slider callback simply deletes the last slice plane and adds a new one at the chosen location. Here's the code.

dbtype volvisGUI

1     function volvisGUI(x,y,z,v)
2     % volvisGUI Interactive volume visualization.
3     %
4     %   Example:
5     %   [x,y,z,v] = flow;
6     %   volvisGUI(x,y,z,v)
7     
8     %   Copyright 2007 The MathWorks, Inc.
9     
10    %% Initalize visualization
11    figure;
12    s = volumeVisualization(x,y,z,v);
13    s.addSlicePlane(s.xMin);
14    
15    %% Add uicontrol
16    hSlider = uicontrol(...
17        'Units','normalized', ...
18        'Position',[.75 .05 .2 .05], ...
19        'Style','slider', ...
20        'Min',s.xMin, ...
21        'Max',s.xMax, ...
22        'Value',s.xMin, ...
23        'Callback',@updateSliderPosition);
24    
25    %%
26        function updateSliderPosition(varargin)
27            s.deleteLastSlicePlane();
28            x = get(hSlider,'Value');
29            s.addSlicePlane(x);
30        end
31    
32    end
33    

Modular Programming Architecture

I've demonstrated here a programming architecture, using nested functions, where you can create individual tools, and incorporate them into larger applications.

I have the visualization engine and I showed how to use it both in a stand-alone way as well as how to incorporate it into a GUI. It now can operate as a tool ready to incorporate in other applications. For example, I could have chosen it to power an animation. Or I could add other degrees of freedom in a GUI by supplying a way to choose other colormaps perhaps.

The Files

You can find the running files for this example here on the File Exchange at MATLAB Cental.

How might you use such ideas for creating tools you can reuse and share easily? Post your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Some functions have an option you can set that displays information throughout the course of a calculation, for example, fminsearch via the OutputFcn in optimset. What are your choices when you want some intermediate output from your own code? I see three choices.

Graphical Output Option Using waitbar

Here's a very short illustration of the kind of control you can get using waitbar. First, create a waitbar and get its handle.

h = waitbar(0,'Initializing waitbar...');

Here's an update when the calculation is halfway done. Make sure to use the waitbar handle to update the existing waitbar.

waitbar(0.5,h,'Halfway there...')

And further along again:

perc = 75;
waitbar(perc/100,h,sprintf('%d%% along...',perc))

Close the progess meter when the calculation is done.

close(h)

There are a plethora of alternatives to waitbar on The MATLAB Central File Exchange. I count 29 today when I search for waitbar.

Textual Output Option

fprintf can print output to the command window by specifying the file identifier as the value 1. In addition, you can take advantage of escape characters, such as \b for backspace, to control where the output goes.

Look at lines 2 and 4 in the code snippet here. On line 2, I print out a line, but do not terminate it with a newline character \n. On line 4, we use the backspace capability to overwrite our previous integer, before pausing so the command window output can be updated. Finally, when the loop is complete, I use fprintf to advance the cursor to a new line \n.

dbtype cmdwinProgress
clc

1     % Example Code for Printing Progress to the Command Window
2     fprintf(1,'here''s my integer:  ');
3     for i=1:9
4         fprintf(1,'\b%d',i); pause(.1)
5     end
6     fprintf('\n')

Here's the final result from that command window, having first made sure the screen (using clc) before running cmdwinProgress.

Progress Indicators

What do you look for in a progress indicator? Do any of the options mentioned above work for you? If not, why not? Post your thoughts here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Repeated Indices with a Single Assignment

I initialize an output matrix to zeros and create two vectors, one with locations (indices) and one for the corresponding data. Now let's look at the results of an indexed assignment.

A = zeros(3);
loc = [1 3 7 9 1];
data = 1:5;
A(loc) = A(loc) + data

A =
     5     0     3
     0     0     0
     2     0     4

You see in this case that even though index 1 was repeated, its output value is only 1, and not 2. Let me reset A so I can show you what's going on.

A = zeros(3);

Here's what "essentially" happens during an assignment in MATLAB. On the right-hand side, A(loc) is created.

A(loc)

ans =
     0     0     0     0     0

This means that there are 2 copies of A(1) currently and an element of data will be added to each of these. The important thing to realize is that now the values for the right-hand side have been calculated and set.

A(loc) + data

ans =
     1     2     3     4     5

Following this, the 5-element right-hand side is assigned into A(loc) on the left-hand side. Values are currently placed into the output vector by marching down the left-hand side locations, replacing one element at a time. So the first element gets replaced not once, but twice, with the final one being the value that you see.

A(loc) = A(loc) + data

A =
     5     0     3
     0     0     0
     2     0     4

Repeated Indices using a for Loop for the Assignment

Here's another way to do the assignment, but in this case, the

B = zeros(3);
for n = 1:length(loc)
  B(loc(n)) = B(loc(n)) + data(n);
end
B

B =
     6     0     3
     0     0     0
     2     0     4

Notice that the output here for B does not match that from the vectorized assignment above. This is because the values for the right-hand side get updated one at a time, and can therefore accumulate as additions made to a particular index.

isequal(A,B)

ans =
     0

Using accumarray to Accumulate Results with Repeated Indices

What if you want the accumulation effect, but you'd prefer to use vectorized code? That's what the function accumarray is for.

C = accumarray(loc',data,[9 1]);
C = reshape(C,[3 3])

C =
     6     0     3
     0     0     0
     2     0     4

We see the same results here as those from the for loop example.

isequal(B,C)

ans =
     1

Repeated Index Usage

Did you know this about MATLAB already? I'm curious to hear your experiences in this area. Post them here.

Get the MATLAB code

Published with MATLAB® 7.4

GOTO Humor

For those of you who wonder why MATLAB does not have a goto statement, you might look at this cartoon in addition to reading an earlier blog article I wrote. I am amused to see this confession about that post.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Code That's Not Recommended

The following code doesn't catch dimensions higher than 2.

 if any(any(input))
     output = doAlg(input);
 end

This code is also not recommended because n will be a scalar even if input has more than 2 dimensions. Depending on what is done in the loop, the output might be the wrong shape/dimension (thanks to Jeremy so I could clarify this).

 [m,n] = size(input);
 for i=1:n
     doSomething;
 end

Code That is Safe for Higher Dimensions

This code should work fine for any input.

 sz = size(input);
 if any(input(:))
     output = doAlg(input);
 end
 output = reshape(output,sz);

Benefits

If you consider N-dimensional input when you write functions, you are less likely to find bugs or wrong answers when you go to use it later.

   * The code is easier to understand without so many any statements.
   * Code is robust to N-dimensional inputs.

Your Advice?

What techniques have you used to ensure your code works with all array shapes? List them here.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

By Default, Column-wise Operations

By default, MATLAB performs many operations by treating the columns as individual vectors and acting on them. However, if the array has only a single element per column, then MATLAB performs the operation along the first non-singleton dimension.

User Issue and Solution

So, that's the default behavior. However, for this user's application, he always wanted the sum to be along the first dimension (down the columns), even if there was only a single entry per column.

To accomplish this, the altered code used the optional second input argument, dimension. This works similarly with other functions that typically reduce the dimensionality from the input to the output, such as

Here's the relevant portion of the help for sum. I happen to know that I only need to go up through the second instance of DIM.

h = help('sum');
f = strfind(h,'DIM');
disp(h(1:f(2)+5))

 SUM Sum of elements.
    S = SUM(X) is the sum of the elements of the vector X. If
    X is a matrix, S is a row vector with the sum over each
    column. For N-D arrays, SUM(X) operates along the first
    non-singleton dimension.
    If X is floating point, that is double or single, S is
    accumulated natively, that is in the same class as X,
    and S has the same class as X. If X is not floating point,
    S is accumulated in double and S has class double.
 
    S = SUM(X,DIM) sums along the dimension DIM. 

Examples

Let's see this behavior in action.

n = 4;
A = cat(3, pascal(4), magic(4), invhilb(4), hadamard(4), hilb(4))

A(:,:,1) =
     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20
A(:,:,2) =
    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1
A(:,:,3) =
          16        -120         240        -140
        -120        1200       -2700        1680
         240       -2700        6480       -4200
        -140        1680       -4200        2800
A(:,:,4) =
     1     1     1     1
     1    -1     1    -1
     1     1    -1    -1
     1    -1    -1     1
A(:,:,5) =
    1.0000    0.5000    0.3333    0.2500
    0.5000    0.3333    0.2500    0.2000
    0.3333    0.2500    0.2000    0.1667
    0.2500    0.2000    0.1667    0.1429

Sum A down the columns.

sumA = sum(A)

sumA(:,:,1) =
     4    10    20    35
sumA(:,:,2) =
    34    34    34    34
sumA(:,:,3) =
    -4    60  -180   140
sumA(:,:,4) =
     4     0     0     0
sumA(:,:,5) =
    2.0833    1.2833    0.9500    0.7595

Sum A along the rows.

sumA2 = sum(A,2)

sumA2(:,:,1) =
     4
    10
    20
    35
sumA2(:,:,2) =
    34
    34
    34
    34
sumA2(:,:,3) =
    -4
    60
  -180
   140
sumA2(:,:,4) =
     4
     0
     0
     0
sumA2(:,:,5) =
    2.0833
    1.2833
    0.9500
    0.7595

Sum A across the third dimension.

sumA3 = sum(A,3)

sumA3 =
  1.0e+003 *
    0.0350   -0.1155    0.2453   -0.1248
   -0.1125    1.2123   -2.6858    1.6912
    0.2513   -2.6888    6.4912   -4.1788
   -0.1338    1.6972   -4.1758    2.8221

Try summing A along the 7th dimension.

sumA7 = sum(A,7)

sumA7(:,:,1) =
     1     1     1     1
     1     2     3     4
     1     3     6    10
     1     4    10    20
sumA7(:,:,2) =
    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1
sumA7(:,:,3) =
          16        -120         240        -140
        -120        1200       -2700        1680
         240       -2700        6480       -4200
        -140        1680       -4200        2800
sumA7(:,:,4) =
     1     1     1     1
     1    -1     1    -1
     1     1    -1    -1
     1    -1    -1     1
sumA7(:,:,5) =
    1.0000    0.5000    0.3333    0.2500
    0.5000    0.3333    0.2500    0.2000
    0.3333    0.2500    0.2000    0.1667
    0.2500    0.2000    0.1667    0.1429

How Can That Be?

Why did I get no error for summing along dimension 7 when I only have a 3 dimensional array? The reason is that MATLAB treats all arrays as having an infinite number of dimensions, most of them trailing singletons. So if I sum along the 7th dimension here, I am only summing single elements, or effectively in this case just returning the original array.

Ever Used Dimension to Your Advantage?

Have you ever used the dimension argument to your advantage? Or tripped over the issue of having an unexpected singleton dimension? Any thoughts, please post them here.

Get the MATLAB code

Published with MATLAB® 7.4

Thanks. I look forward to meeting several of you later this summer.

Get the MATLAB code

Published with MATLAB® 7.4

Contents

Two Conventions

When calling functions in MATLAB that create or reference arrays, there are two difference conventions for specifying the first two dimensions of an array:

These different conventions get used depending on whether the array is strictly an array or when it has a probable physical interpretation in terms of Cartesian grid.

Matrix versus Cartesian Frameworks

In MATLAB, we index into matrices using row and column indices. That is, for matrix

A = magic(3)

A =
     8     1     6
     3     5     7
     4     9     2

I can reference the last element in the second row like this:

A(2,3)

ans =
     7

or

A(2,end)

ans =
     7

and to retrieve the first element of columns 1 and 3, I type this:

A(1,[1 3])

ans =
     8     6

Check out this post to get a fairly complete picture of indexing in MATLAB.

Plotting Matrices

The options xy and ij for axis allow us to control whether to view the data just the way the matrix is displayed (the (1,1) element is in the upper left corner):

colormap(jet(9))
image(A)
colorbar

or we can view the image/matrix as a top-down view of a physical "map" with the y-values now increasing along the vertical axis.

axis xy

In images, the matrices column values plot vertically. In line plots, the columns are plotted vs. their indices, making columns plot horizontally.

h = plot(A);
legend(h,{'Column 1','Column 2','Column 3'})

Creating an N-Dimensional Grid

ndgrid uses matrix indexing since it is designed to suit working in any number of array dimensions.

[ii, jj] = ndgrid(1:3,1:2)

ii =
     1     1
     2     2
     3     3
jj =
     1     2
     1     2
     1     2

Creating an 2- or 3-Dimensional Grid

On the other hand, mesh, interprets its first inputs in a Cartesian framework since it is used primarily for working with 2- and 3-dimensional plots, often depicting something over an area or volume.

[xx, yy] = meshgrid(1:3,1:2)

xx =
     1     2     3
     1     2     3
yy =
     1     1     1
     2     2     2

In the Remarks section of the documentation for meshgrid, you will find a discussion of this, especially with respect to inverting the first two input arguments compared to ndgrid. There is related discussion in the Description section for the function mesh.

Functions designed for two or three dimension, such as contour and quiver, usually work in Cartesian coordinates. Functions that are designed to work in N-dimensions generally work with a matrix indexing framework. Be sure to check the help for functions you use to see which convention applies.

Thoughts?

Any ideas on what we can do, perhaps in examples or documentation (or eslewhere?) to help people out with this potential confusion? If so, please post here.

Get the MATLAB code

Published with MATLAB® 7.4