Contouring
From ComputingForScientists
Contents 
1. References
 Overviews of algorithm: [1] [2] [3]
 Paper on Marching Cubes Algorithm: [4]
 MATLAB implementation of contouring: [5]
 Stack Overflow Discussion: [6]
 Chapter 14 of Mathematical Principles for Scientific Computing and Visualization: PDF
 Chapter 5 of Data Visualization: Principles and Practice: PDF
2. Overview
Consider a set of measurements of altitude z made on a uniform grid that extends in the x and y directions. These measurements are stored in an 6x6
array as shown in the image below. Note that the bottom and left axes have labels for the x and y position of the measurement and the top and right axes have labels indicating the location of the measurement of z in the matrix.
Countouring algorithms seek to draw lines that correspond to locations of constant altitude, even in the absence of measurements at a given location.
In the following image, a horizontal line was drawn to connect a point i,j
when its value was 1
and the value at i1,j
is 1
. This occurred at two locations: i=4,j=5
and i=4,j=2
.
To add additional lines, we would need to find locations where a given value in the matrix is 1
and the value of any of its neighbors (including diagonal neighbors) is 1
.
To finish off this example, we would need to search for the all of the other possible scenarios (the first one below is the one considered):

z(i,j) = 1
andz(i1,j) = 1
(Left) 
z(i,j) = 1
andz(i+1,j) = 1
(Right) 
z(i,j) = 1
andz(i,j+1) = 1
(Up) 
z(i,j) = 1
andz(i,j1) = 1
(Down) 
z(i,j) = 1
andz(i1,j1) = 1
(Leftdown) 
z(i,j) = 1
andz(i1,j+1) = 1
(Leftup) 
z(i,j) = 1
andz(i+1,j1) = 1
(Rightdown) 
z(i,j) = 1
andz(i+1,j+1) = 1
(Rightup)
Any time one of the scenarios was found, we would need to draw a line connecting the points associated with that scenario. For example, at i,j = 2,3
z = 1
, and at i+1,j+1
, z=1
. These two points match the last scenario, so we could connect a line from i,j=2,3
to i,j=3,4
.
The point i,j=2,5
does not match any of these scenarios, so no line will be drawn connecting this point to another point.
In problem I below, you are asked to generalize the program so that it checks for all eight cases. The sample code only looks for two cases and draws lines accordingly.
3. More general algorithm
In the previous example, we selected a contour altitude that corresponded to a value that existed in the set of measurements. More generally, we will want to select an arbitrary value, say z=1.5
and draw lines corresponding to the location where we would expect this to be approximately true.
To do this, we visit each point and check all 8 of the following scenarios:

z(i,j) <= 1.5
andz(i1,j) > 1.5
. 
z(i,j) > 1.5
andz(i1,j) <= 1.5
. 
z(i,j) <= 1.5
andz(i+1,j) > 1.5
. 
z(i,j) > 1.5
andz(i+1,j) <= 1.5
. 
z(i,j) <= 1.5
andz(i,j+1) > 1.5
. 
z(i,j) > 1.5
andz(i,j+1) <= 1.5
. 
z(i,j) <= 1.5
andz(i,j1) > 1.5
. 
z(i,j) > 1.5
andz(i,j1) <= 1.5
.
When any of these scenarios are matched, we draw a dot halfway between the two points that were considered. The following image shows the result of the calculation for this problem. Note that technically, we should not always place the dot halfway between the two points considered  for example, the dot between z = 1.0 and z = 1.8 should be closer to z = 1.8 than it is to z = 1.0. The location can be determined by using the formula
ydot = y(i) + dy*(1.01.5)/(1.81.0)
where dy is the distance between the two dots. In the case for the magneta dot drawn below the blue dot at x = 80, y = 40, the position would be
ydot = 40 + 20*(1.01.5)/(1.81.0) = 27.5
instead of at at y = 30, which is halfway between.
4. Problems
4.1. Manual Contour Line Drawing I
The following program first checks each element of the matrix z
at positions i
and j
for a match to the value 1
.
If it finds a position with this value, it checks
 if the value at the position to the left (
i1
) is also1
. If it is, a line is drawn to connect the point ati,j
toi1,j
.  if the value at the position to the right and up (
i+1,j+1
) is also1
. If it is, a line is drawn to connect the point ati,j
toi+1,j+1
.
Modify this program so that it connects all points that satisfy the 8 scenarios given at #Overview.
When executed, your code should produce lines corresponding to the shape of an octagon.
clear; x = [0:20:100]; % An array of x coordinates. y = [0:20:100]; % An array of y coordinates. % Create a matrix of altitude values. for i = 1:length(x) for j = 1:length(y) z(i,j) = (x(i)50)^2 + (y(j)50)^2; end end z = z/1000; % Open a figure window if it does not exist and then clear it. figure(1);clf; % Set axis limits wider than default axis([20 120 20 120]); % Each time a plot command is issued, MATLAB clears the figure window % unless hold on was executed. hold on; % Draw a dot at position x(i), y(i). % Display text at position x(i), y(i). % Offset text by 2 units to it does not overlap dot. for i = 1:length(x) for j = 1:length(y) plot(x(i),y(j),'b.','MarkerSize',20); text(2+x(i),2+y(j),['z=',num2str(z(i,j))],'Color','blue'); end end % Add axes labels. ylabel(gca, 'y','Rotation',0); xlabel(gca,'x'); % Make the aspect ratio of the axes square. axis square % Set tick locations manually. set(gca,'XTick',[0:20:100]) set(gca,'YTick',[0:20:100]) % Show a background grid. grid on; % Add a label indicating the function that is plotted. text(25,114,'Values of z=[(x50)^2+(y50)^2]/1000') % Draw contour lines for z = 1. Two of the possible eight cases are checked for below. for i = 2:length(x)1 for j = 2:length(y)1 if (z(i,j) == 1) if (z(i1,j) == 1) % Left % The following plot command will only execute if z(i,j) == 1 and z(i1,j) == 1 are true plot([x(i),x(i1)],[y(j),y(j)]); % Draw line that connects x(i), y(i) to x(i1), y(i) end if (z(i+1,j+1) == 1) % RightUp % The following plot command will only execute if z(i,j) == 1 and z(i+1,j+1) == 1 are true plot([x(i),x(i+1)],[y(j),y(j+1)]); % Draw line that connects x(i), y(i) to x(i+1), y(i+1) end end end end % Don't worry about understanding this part. It just makes % the labels look sensible. % Add index axes. % Based on technique from http://stackoverflow.com/questions/19569134/overlayingtwoaxesinamatlabplot % Create a new axes in the same position as the first one, overlaid on top h2 = axes('position', get(gca, 'position')); % Put the new axes' y labels on the right, set the x limits the same as the % original axes', and make the background transparent set(h2, 'YAxisLocation', 'right', 'YLim',[0,length(y)+1],'Xlim', [0,length(x)+1], 'Color', 'none'); ylabel(h2, 'j','Rotation',0); xlabel(h2, 'i'); set(h2,'XTick',[1:length(x)]) set(h2,'YTick',[1:length(y)]) set(h2, 'XAxisLocation', 'top'); axis square; hold on; print dpng contour_algorithm_basics.png print deps contour_algorithm_basics.eps
Answer 

clear; x = [0:20:100]; % An array of x coordinates. y = [0:20:100]; % An array of y coordinates. % Create a matrix of altitude values. for i = 1:length(x) for j = 1:length(y) z(i,j) = (x(i)50)^2 + (y(j)50)^2; end end z = z/1000; % Open a figure window if it does not exist and then clear it. figure(1);clf; % Set axis limits wider than default axis([20 120 20 120]); % Each time a plot command is issued, MATLAB clears the figure window % unless hold on was executed. hold on; % Draw a dot at position x(i), y(i). % Display text at position x(i), y(i). % Offset text by 2 units to it does not overlap dot. for i = 1:length(x) for j = 1:length(y) plot(x(i),y(j),'b.','MarkerSize',20); text(2+x(i),2+y(j),['z=',num2str(z(i,j))],'Color','blue'); end end % Add axes labels. ylabel(gca, 'y','Rotation',0); xlabel(gca,'x'); % Make the aspect ratio of the axes square. axis square % Set tick locations manually. set(gca,'XTick',[0:20:100]) set(gca,'YTick',[0:20:100]) % Show a background grid. grid on; % Add a label indicating the function that is plotted. text(25,114,'Values of z=[(x50)^2+(y50)^2]/1000') % Draw contour lines for z = 1. Two of the possible eight cases are checked for below. for i = 2:length(x)1 for j = 2:length(y)1 if (z(i,j) == 1) if (0) % Replace 0 with 1 to cause this block to execute. % This approach uses a for loop to do a similar operation as % the version below. Note the total number of iterations is % 9 (3 in ii and 3 in jj). The test % && ~(ii == 0 && jj == 0) prevents a "line" from being drawn from % x(i),y(i) to x(i),y(i), which is not one of the cases in the % original program. for ii = 1:1 for jj = 1:1 if (z(i+ii,j+jj) == 1) && ~(ii == 0 && jj == 0) % Draw line that connects x(i), y(i) to x(ii), y(jj) plot([x(i),x(i+ii)],[y(j),y(j+jj)]); end end end end if (1) % Replace 1 with 0 to prevent this block from executing. if (z(i1,j1) == 1) % LeftDown plot([x(i),x(i1)],[y(j),y(j1)]); end if (z(i1,j) == 1) % Left plot([x(i),x(i1)],[y(j),y(j)]); end if (z(i1,j+1) == 1) % LeftUp plot([x(i),x(i1)],[y(j),y(j+1)]); end if (z(i,j1) == 1) % Down plot([x(i),x(i)],[y(j),y(j1)]); end if (z(i,j+1) == 1) % Up plot([x(i),x(i)],[y(j),y(j+1)]); end if (z(i+1,j1) == 1) % RightDown plot([x(i),x(i+1)],[y(j),y(j1)]); end if (z(i+1,j) == 1) % Right plot([x(i),x(i+1)],[y(j),y(j)]); end if (z(i+1,j+1) == 1) % RightUp plot([x(i),x(i+1)],[y(j),y(j+1)]); end end end end end % Don't worry about understanding this part. It just makes % the labels look sensible. % Add index axes. % Based on technique from http://stackoverflow.com/questions/19569134/overlayingtwoaxesinamatlabplot % Create a new axes in the same position as the first one, overlaid on top h2 = axes('position', get(gca, 'position')); % Put the new axes' y labels on the right, set the x limits the same as the % original axes', and make the background transparent set(h2, 'YAxisLocation', 'right', 'YLim',[0,length(y)+1],'Xlim', [0,length(x)+1], 'Color', 'none'); ylabel(h2, 'j','Rotation',0); xlabel(h2, 'i'); set(h2,'XTick',[1:length(x)]) set(h2,'YTick',[1:length(y)]) set(h2, 'XAxisLocation', 'top'); axis square; hold on; print dpng contour_algorithm_basics.png print deps contour_algorithm_basics.eps 
4.2. Manual Contour Line Drawing II
In the following image, magenta dots were placed halfway between two points that satisfied one of the eight conditions decribed above.
Calculate (by hand) the locations of each of the dots by accounting for how close countour value of z = 1.5 is to the two points on either side of the magneta dot. (See the notes above for an example.)
4.3. Manual Contour Line Drawing III
Suppose that you wanted to give instructions to another student for drawing a contour line by hand. You give them a printout of images that look similar to that shown below, except the z
values may be different. Write out instructions for drawing the contour lines (assuming that the dots to be connected are always placed halfway between neighboring dots).
Give a printout of the image below and your instructions to a friend and ask them to draw the contour lines for z=1.5
and z=0.6
. Turn in a copy of your instructions and your friend's solution.
4.4. Manual Contour Line Drawing IV
Draw dots for z=1.5
for the matrix in Part I above using MATLAB; place the dots halfway between grid point locations (that is, create the image below).
4.5. Manual Contour Line Drawing V
Draw dots for z=1.5
for the matrix in Part I above using MATLAB, but place the dots at the appropriate location between the grid points instead of halfway between.
4.6. Manual Contour Line Drawing VI
Draw lines that connect the dots computed in the previous problem.