# Just Launched: You can now import projects and releases from Google Code onto SourceForge

We are excited to release new functionality to enable a 1-click import from Google Code onto the Allura platform on SourceForge. You can import tickets, wikis, source, releases, and more with a few simple steps.

## [84a6ed]: PDL / Book / graphics_3d.pod Maximize Restore History

Parent: [3f1bb4] (diff)

Child: [b94816] (diff)

Download this file

### graphics_3d.pod    586 lines (446 with data), 21.2 kB

=pod

=head1 3D Graphics with OpenGL

=head2 Introduction

=for html <img WIDTH=300 src="graphics_3d/gnuplot-eg.png">

=for man
+-----------------------------+
|,,'........................',|
|                             |
|                             |
|           ...  ...          |
|   .  .',,'.  .'......       |
|   . ..  .';,'... .....      |
|   .      .......'           |
|             ...             |
|                             |
|                             |
|                             |
|. . . . . .                  |
+-----------------------------+

Figure 3.1: A 3D surface graph plotted using gnuplot,
using the commands:

set isosamples 30; splot [0:7] [0:7] sin(x)*sin(y).

There are lots of programs that let you plot so-called 3D surface
graphs, such as the one shown in Fig. 3.1. However, from
the beginning, PDL's 3D graphics have had something different
that we feel is really useful: motion, or as we call it "twiddling".
Dragging the 3D image with the mouse rotates the image, at the
speed allowed by your display hardware. This turned
out to be quite useful for displaying functions: the human eye is
able to grasp the presented 3D surface much better when it moves,
especially in response to the mouse.

=for html <img WIDTH=200 src="graphics_3d/pdltrid-eg1.png">

=for html <img WIDTH=200 src="graphics_3d/pdltrid-eg2.png">

=for man
+-----------------------------++-----------------------------+
|                             ||                             |
|                             ||                             |
|                             ||        . :c.                |
|     .                       ||          .cd,               |
|          ,.     .           ||          .'ldl.'            |
|     .   .od   ..c:          ||           .;l.'od.          |
|        .'dkl....,;.         ||       c.   .. .;ox,         |
|     .    ..,;..'. dc        ||     .dxc.      .:oxdd,      |
|              .,   :.        ||       ,o:.     .'lxxx:      |
|     .   ..   .              ||         ,ol.    .;dkd;      |
|    ...                      ||      .    ;c'.    .. .      |
|           .                 ||             ..              |
|                             ||                             |
|                             ||                             |
+-----------------------------++-----------------------------+

Figure 3.2: The same 3D surface, plotted using the
PDL::TriD module. The two different images were obtained
literally by grabbing the image in the window opened by
PDL and dragging it with the mouse to rotate.

Let's start with plotting the surface we showed using gnuplot in
the beginning:

pdl> use PDL::Graphics::TriD;
pdl> $x = xlinvals(zeroes(30), 0, 7); pdl> imag3d [ sin($x) * sin($x->dummy(0)) ]; This should produce a new window with the image seen in Figure 3.2. Notice that your console window is now frozen: it is waiting for you to twiddle in the graphics window using the mouse and to press C<q> in that window once you're done. If the above commands produce an error instead of a new window, it might be that your PDL wasn't compiled with the option to include the 3D graphics library. (See the Perl Data Language web site at L<http://pdl.perl.org> for information on installing and using PDL.) That above expression is a bit more difficult than the gnuplot version, and there's a simple reason for that: gnuplot is primarily meant for plotting functions; PDL is meant for handling and plotting numerical data. So to plot a function, we have to create the data for the function first which is a bit more difficult. Now let's go through that part by part. pdl> use PDL::Graphics::TriD; The first line simply tells Perl to load the PDL::Graphics::TriD module. The name comes from the fact that you can't have parts of module names starting with numbers, unfortunately. The second line pdl>$x = xlinvals(zeroes(30), 0, 7);

creates a one-dimensional piddle with 30 elements that has linear
values from 0 to 7:

pdl> p $x [0 0.24137931 0.48275862 0.72413793 0.96551724..... The C<xlinvals> and the corresponding C<ylinvals> and C<zlinvals> are useful for exactly this purpose: creating piddles of equally spaced values. The final line, pdl> imag3d [ sin($x) * sin($x->dummy(0)) ]; is what draws the actual image. The expression inside, uses the variable C<$x> for both the X and Y coordinates, via a clever use of
the dummy operation. (See chapter [chap_slice] for some
explanation).  This results in a 2-dimensional piddle with the
values for the Z coordinate. So far you've already seen all this.
And the final part, C<imag3d [vals]> is the call that creates the 3D
plot and opens the new window for it. The brackets around the
parameter may be slightly surprising: the 2-D commands work well
without those but there is a good reason for this, as you'll
learn later on: otherwise there would be a bad ambiguity.

=cut

=head2 Parametric Graphics

We alluded in the introduction that allowing

pdl> imag3d $piddle; could be ambiguous and should be written pdl> imag3d [$piddle];

if C<$piddle> is intended to be the Z axis values of a rectangular 2D plot. Now is the time to find out why. The simple truth is that pdl> imag3d$piddle;

is in fact legal code---if and only if the first dimension of
C<$piddle> has exactly three elements. As you probably have already guessed, these three elements are X, Y and Z. So what you can do is pass C<$piddle> with shape C<[3,t,u]> which is the same as a
2-dimensional C<[t,u]> lattice with a 3-vector at each point. This
piddle will then be interpreted parametrically: the mesh will be
drawn as a function of C<$t> and C<$u>.

Let's have an example: a curve that is not possible to plot with
just Z axis values, say the surface of a torus, with colors
coming from somewhere. First, set up the piddles and the
parameter variables:

use PDL;
use PDL::Graphics::TriD;
use PDL::NiceSlice;

$torus = zeroes(3, 60, 20);$x = $torus((0));$y = $torus((1));$z = $torus((2));$t = xlinvals $x, 0, 6.28;$u = ylinvals $x, 0, 6.28; Note that the coordinate separation can be done in just one line: ($x, $y,$z) = map { $torus(($_)) } 0..2;

Next, we color the torus. Let's put stripes on it:

$r = (1 + sin(2*$t +   $u))/2;$g = (1 + cos(2*$t + 2*$u))/2;
$b = (1 + sin(2*$t + 3*$u))/2; Then, we choose the outer and inner radii and put the coordinates into the slices. We'll let the torus lie in the XY plane so the parametric coordinates can be easily derived.$r_o = 3;
$r_i = 1;$x .= ($r_o +$r_i * sin($u)) * sin($t) ;
$y .= ($r_o + $r_i * sin($u)) * cos($t);$z .=         $r_i * cos($u);
imag3d_ns $torus, [$r, $g,$b];

And here's our colorful torus!

=for html <img WIDTH=300 src="graphics_3d/color-torus.png">

It looks a bit more like a barrel
because TriD automatically scales the axes but there it is. Note
how we use C<imag3d_ns> to get the colors instead of the shaded
version.

Now, there is more than one way to do it. If your data is not by
default in the three-vector format (as ours wasn't above), it is
probably easier to do

imag3d_ns [$x,$y, $z], [$r, $g,$b];

which will produce the same results. Also, we could concatenate
the RGB piddles to form a single C<[3,60,20]> piddle that could be
used without square brackets:

$rgb = cat($r,$g,$b)->mv(-1,0);  # $rgb is [3,60,20] imag3d_ns$torus, $rgb; Now, since PDL does its best to make dimensions usable anywhere, we can easily plot several parametrics of the same parameters at once, if we pack all the surfaces into a piddle of shape C<[3,n_t,n_u,...]> where the three periods in the end indicate the beginning of the extra parameters. For example, we can plot a family of shrinking toruses by adding an extra dimension into C<$torus>:

$cone =$torus->dummy(3, 4)->copy();
$fac = axisvals($cone, 3);
$cone *=$fac + 2;
$cone(2) += 4 *$fac;
imag3d $cone; =for html <img WIDTH=300 src="graphics_3d/torus-stack.png"> And further, if we want to distort them, it's perfectly possible:$x = $cone(0);$cone(2) += 0.1 * $x ** 2; imag3d$cone;

=for html <img WIDTH=300 src="graphics_3d/torus-stack-warp.png">

Any other kind of mutilation is also possible but we leave you to
discovering the interesting things that are possible by yourself,
because we have to move to something else that's important to
cover: coordinate systems. So far, all the examples you've seen
have happened in the Euclidean coordinate system where the
coordinates are specified as measures X, Y and Z on three
orthogonal axes.

Or actually this is not true: in fact, we have used two kinds of
coordinates, the explicit X, Y and Z given in this section but in
the preceding sections, only Z has been given and X and Y have
been assumed by the system from the context.

Of course, since PDL tries to follow "simple things simple,
complicated things possible", it is possible to override the
default context.

=head2 Types of 3D Graphical Objects

So far, we've only been toying with surfaces. However, PDL can do
much more. We can plot points; here's a picture of two samples
from different (overlapping) probability distributions, plotted
with different colors:

use PDL::Graphics::TriD;

$i = zeroes(8000);$which = random($i) < 0.5;$x = grandom($i) * (1 +$which);
$y = grandom($i) * (0.5 + $which);$z = grandom($i) * (2 -$which);
$x +=$which * $y;$y += $which *$z; # Make it oblique
points3d [$x,$y, $z], [$which, 0.5*(1-$which), 1-$which];

And the result:

=for html <img WIDTH=300 src="graphics_3d/two-prob-distributions.png">

A lot of fun things can be done with points but we'll go into
that later.

Then, there are---of course---lines. As a fun demo of lines,
let's plot a number of flowlines moving in the Lorenz attractor.
As you may know, the Lorenz attractor is described by

=for asciitex
\begin{array}{ccl}
\frac{dx}{dt} & = & \sigma (y-x)\\
\frac{dy}{dt} & = & (r-z) x - y\\
\frac{dz}{dt} & = & (y-b) z\\
\end{array}

dx
--    =    sigma (y - x)
dt

dy
--    =    (r - z) x  -  y
dt

dz
--    =    (y - b) z
dt

where sigma=10, r=28 and b=8/3. Because we're just doing this
as a simple demo, we'll use the extremely unstable d=Delta
method integration. We'll plot six trajectories that start close
to each other.

use PDL::Graphics::TriD;
$n = 500;$nstart = 0;
$nc = 6;$delta = 0.015;
# $x = pdl(1, 1, 1, 1, 1); #$y = pdl(1, 1, 1, 1, 1);
# $z = pdl(1, 1.01, 1.02, 1.03, 1.04);$xs = zeroes($n,$nc);
$ys = zeroes($n, $nc);$zs = zeroes($n,$nc);
$x = -23 * ones($nc);
$y = -2 * ones($nc);
$z = 20 * ones($nc) + 0.02 * xvals($nc);$sigma = 10; $r = 28;$b = 8.0/3.0;
for (-$nstart..$n-1) {
if($_ >= 0) {$xs(($_)) .=$x;
$ys(($_)) .= $y;$zs(($_)) .=$z;
}
$dx =$sigma * ($y -$x);
$dy = ($r - $z)*$x - $y;$dz = $x*$y - $b *$z;

$x +=$delta * $dx;$y += $delta *$dy;
$z +=$delta * $dz; }$col = yvals(1, $nc) / ($nc-1);
$tim = xvals($n) / ($n-1); line3d [$xs, $ys,$zs], [$col,$tim , 1-$col]; =cut =pod =for html <img WIDTH=300 src="graphics_3d/busy-lorenz-attr.png"> =cut =pod Figure: Busy Lorenz Attractor Unfortunately, this plot has too much stuff going on so it's difficult to see where the functions diverge even though they have different colors at different times. This is an excellent time to change variables: let's get rid of X and plot the time step instead: line3d [$tim, $ys,$zs], [$col,$tim , 1-$col]; This yields a much clearer plot of the chaotic behaviour when the lines diverge with time. =for html <img WIDTH=300 src="graphics_3d/timeseries-lorenz.png"> In the latest versions of PDL it is possible to adjust the line width as well: line3d [$tim, $ys,$zs], [$col,$tim , 1-$col], {LineWidth => 10} gives the same plot but with much thicker lines. =for html <img WIDTH=300 src="graphics_3d/timeseries-lorenz-thick.png"> The basic rectangular surface you already saw in the preceding sections. It also has an option to turn off the lines. There is also a command C<mesh3d> similar to the C<imag3d> surface which just draws the surface as a wire mesh instead of a solid surface. On slow machines this can be of great help. Finally, there are two commands for quickly painting strictly rectangular truecolor images: C<imagrgb> and C<imagrgb3d>. This can be demonstrated by Tuomas J. Lukka's 4-liner: use PDL; use PDL::Graphics::TriD;$a=zeroes 300,300;$r=$a->xlinvals(-1.5,
0.5);$i=$a->ylinvals(-1,1);$t=$r;$u=$i;for(1..30){$q=$r**2-$i**2+$t;$h=2 *$r*$i+$u;$d=$r**2+$i**2;$a=lclip($a,$_*($d>2.0)*($a==0));($r,$i)=map{$_ ->clip(-5,5)}($q,$h);}imagrgb[$a/30];

This, as odd as it may sound, plots a grayscale Mandelbrot. If
you work your way through the code, you'll see that it simply
iterates the standard Mandelbrot iteration formula

z <-- z**2 + C

where C<C> is the original point. Then it uses C<lclip> to keep the
numbers in a reasonable range and colors the points according to
the iteration when the point crossed the distance C<sqrt(2)> from
the origin. The piddle C<$a> is two-dimensional so just like for coordinates, it is enclosed in an array ref. It is also possible to use imag3gb [$r, $g,$b];
imag3gb $colors; where the RGB piddles are two-dimensional and C<$colors> has three
dimensions, the first of which is of length three.

The command C<imagrgb3d> does the same but allows the user to place
the rectangle anywhere in 3-space. This is useful e.g. for
putting an image underneath a plotted surface of the same
function, as we shall see in the next section.

=head2 More than one Image

If you have used the PDL PGPLOT interface for plotting multiple
graphs then TriD is not going to surprise you: the commands
C<hold3d> and C<release3d> work just like their PGPLOT counterparts.
Before going further, however, let me remind you that for many
plots, it is not necessary to explicitly plot several points,
lines, surfaces or whatever: it can be easier just to use extra
dimensions, like we used for the torus cone in the first section.

However, if you want to put objects of more than one type, or
objects of more than one resolution on the same graph, then you
do need to do so explicitly. As an example we'll use some fractal
mountain code by Tuomas J. Lukka from the 3D Gallery. Unlike with
the mandelbrot that has a well-known algorithm, this code we'd
just better format clearly from the start (the parameters have
also been slightly modified and the code has been modified to
plot all the iterations on top of each other).

use PDL; # XXX FIX - LOOKS BAD.
use PDL::Image2D;
use PDL::Graphics::TriD;
$k = ones(3,3) / 9;$a = 20;
$b =$a*(random(2,2)-0.5);
hold3d(); # Set the coordinate system: XXX hack!!! FIX TriD
line3d pdl([[0, 0, 0,], [0, 0, 10]]);
for (0..4) {
if ($_ != 0) {$c = $b->dummy(0,2)->clump(2)->xchg(0,1)-> dummy(0,2)->clump(2)->xchg(0,1)->copy;$c += $a*($c->random-0.5);
$a /= 1.5;$b = conv2d($c,$k);
}
imag3d [xlinvals($b,0,1), ylinvals($b,0,1), $b + 2.0*$_], {Lines => 0};
}
release3d();

Even laid out bare, this code is a mouthful with that big double
dummy-clump-xchg thing in the middle. But in fact the function is
really simple: the dummy-clump-xchg thing simply doubles the
length of each dimension, copying each value to two consecutive
locations. After doubling the resolution, we add some noise from
the random function (the magnitude of the noise is diminished
each time). Finally, we pull in PDL::Image2D for the C<conv2d>
routine that does 2-dimensional convolutions (optimized for small
kernels like ours). We use a 5x5 kernel to smooth our data at
each step by convolution. That's the numerical part, now here is
the sequence of images created:

=for html <img WIDTH=400 src="graphics_3d/fractal-mountain-sequence.png">

=head2 Putting it all together---cool hacks

Here's one where the original idea is by Robin Williams, done for
the 3D Gallery. This gallery is available in the PDL distribution
in the file Demos/TriDGallery.pm. The idea is to put interesting
scripts that do a lot using just 4 lines of 72 characters. The
crux of the idea is to use OpenGL points to perform volume-like
rendering.  This is just a quick hack. However, the principles
are interesting enough that we thought you might
enjoy them. Let's start with a function of three variables, whose
zeroes are a sphere and an ellipsoid inside the sphere, with the
Y axis slightly distorted to form a parabola with the Z axis:

sub f {
my($x,$y, $z) = @_;$y = $y + 0.04 *$z**2;
return (($x**2 +$y**2 + $z**2) - 100) * ((2*$x**2 + 4*$y**2 + 4*$z**2) - 100);
}

Note here that we can't use the += operator for C<$y> since below we use the same piddle for the three coordinates (with a simple dummy transformation). Now, we want to picture approximately where the function crosses zero, but since there are two separate zero surfaces we can't just use an algorithm that finds a zero and creates an isosurface. Besides, an isosurface renderer wouldn't be able to show both the sphere and the ellipsoid simultaneously. So rather, let's first calculate the sign of the function in a 50x50x50 lattice. The radius of the sphere is C<sqrt(100)=10> so we make the coordinate system slightly larger. use PDL::Graphics::TriD;$x = xlinvals(zeroes(float,50), -11, 11);
$f = f($x, $x->dummy(0),$x->dummy(0)->dummy(0));
$sign = byte($f>0);

Now that we have the sign, why don't we simply find the set of
points where the sign has changed. It is simplest to do this over
just one dimension:

$df = ($sign(0:-2) != $sign(1:-1)); points3d whichND($df);  # for PDL-2.4.10

And indeed, we get a rotatable set of points in 3-space that are
in the shape of a sphere with an ellipsoid inside, slightly
distorted, just as ordered.

=for html <img WIDTH=400 src="graphics_3d/nest-sphere-01.png">

This is not yet a good picture: there
is a hole in the point set where the surface is parallel to the X
axis, naturally, since there there is no difference between the
sigh between the points next to each other on X axis.

B<NOTE>:  For PDL-2.4.9 and earlier, you'll need to use
C<points3d [ whichND($df) ];> since previous to PDL-2.4.10 C<whichND> returned a list of piddles in list context. That behavior is now deprecated. To do a more complete job, we need to compare the signs not only along X but other dimensions as well. This is possible due to the wonderful invention by Robin Williams:$a = $sign; foreach (1,2,4) {$t=($a(0:-2)<<$_);
$t+=$a(1:-1);
$a=$t->mv(0,2);
}
points3d [whichND(($a != 0) & ($a != 255))];

It's a bit cryptic but truly beautiful so bear with us while we
go through it. The loop is executed thrice, once for each
dimension. In the beginning, we know that all the values in C<$a> are either 0 or 1. The first line of the loop takes a slice from C<$a>, leaving the last element of dimension one out and shifts if
by the loop index C<$_>. The second line takes another slice, this time leaving out the first element and adds it to the first. Finally, the dimensions are rotated for the next invocation. Choosing the shifts to be 1,2,4 is the key: this way after the first round, the piddle contains values 0,1,2,3 after the second it contains 0...15 and after the third, 0..255. None of the shifts shift anything on top of each other so the plus operation could be replaced with a bitwise or. So after the loop, we have a three-dimensional piddle with one index less in each dimension, and each value in that piddle contains in its 8 bits the 8 corners of a small cube. Finally, to find whether the function crosses zero at that cube, we simply check whether all the bits are equal, i.e. whether the number is 255 or 0 and if it isn't we know the function changes sign. =for html <img WIDTH=400 src="graphics_3d/nest-sphere-02.png"> The image quality can be slightly improved by removing the Moire effect through randomization: points3d (map {$_+$_->float->random} whichND(($a != 0) & (\$a != 255)))

=for html <img WIDTH=400 src="graphics_3d/nest-sphere-moire.png">

Now, to further improve image quality we could add
different-color pixels but that would require alpha blending to
the OpenGL parameters and this would get into complications we
don't necessarily want here. So now we're going to KISS*
this topic away and move to the next one.

* Keep It Simple, Stupid