Craig DeForest's PDL demos

PDL Demos

Here are some demos of things that I've been doing lately with PDL.

Image warping, revisited

(22-Apr-2005) A couple of weeks ago I got interested in demoing the image warping package, and quickly ginned up a nice movie showing a family of coordinate transformations. The resampling code is now much faster than it was back when the galaxy demo (bottom of this page) was written. Here is an image of the shuttle Atlantis docked with Mir (before Mir re-entered), together with a family of distorted images. The distorted images are being rotated in conformal azimuthal coordinates (CAC) in CA coordinates, the horizontal direction corresponds to azimuthal angle, and the vertical direction corresponds to the logarithm of distance. Hence the rotation couples azimuth to scale, and distance from the origin couples back to azimuthal angle. Features inside of a critical radius rotate counterclockwise; features outside of that radius rotate clockwise. Likewise, stuff that is left of center is constantly shrinking, while stuff that is right of center is constantly growing. There is a scale singularity along a vertical line from the center of the images. Because the transformation is conformal, all of the image features are recognizable.

I believe that this family of distortions is what M.C. Escher used in his famous "Gallery" print of a gallery that contains a picture that, in turn, contains the town that contains the gallery. Putting a scale ambiguity in the image, across the singular line, would give you print with the same basic theme.

(click the images to see a large version of each)

The original image

Some distortion (10^o rotation)

90^o rotation:
(maximal scale weirdness)

180^o rotation:
(inside-out)

I also generated a movie in 1 degree increments of rotation; it's sort of trippy to watch. The movie was created with the set of PDL commands over there --->. Very straightforward to do with the PDL::Transform syntax for manipulating coordinate transformations. Notice that the antialiased sampling does a really good job of hiding the original pixels and also preventing 'shimmer' in the reduced areas of the image.
The transforms themselves are a little complicated as they are working directly in the pixel coordinate system (without autoscaling). The $t variable contains the transformation into CA coordinates (with an initial offset to move the origin and a rotation to make the singular line vertical). $ts just reduces the image from 1500x1200 to just 500x400 -- it's a simple scaling coordinate transform.
use PDL::Transform; $a=rim('shuttle-mir.jpg'); # rim loads RGB image as (1500x1200x3) # Define coordinate transforms $t= t_offset([-$pi,0]) x t_radial(r0=>100) x t_linear(pre=>[-250,-200],rot=>90); $ts = t_linear(scale=>0.33333); # Generate rotated movie frames for $i(0..359) { $b = $a->map( !$t x t_linear( rot=>$i ) x $t x $ts, [500,400], {pix=>1, m=>'g', b=>'mirror', big=>5 } ); wim( $b, sprintf("frame-%3.3d.png",$i )); } # make sample frame - larger image; unscale at end of transform $b = $a->map( !$ts x !$t x t_linear(rot=>180) x $t x $ts, [1500,1200], {pix=>1, m=>'g', b=>'mirror', big=>5 } ); wim($b,'shuttle-mir-180.jpg');

Image despiking

This isn't PDL-specific, just a nice technique that I've recently implemented in PDL. A lot of solar physics data is in the form of image sequences, some of which contain glitches from cosmic rays or radiation that impacted the detector during the observation. Usually spatial filters are used to remove the glitches before analysis. That technique fails when the image itself is too spiky: the filter has difficulty distinguishing solar data from cosmic ray noise.

Zspike is an algorithm that uses the slower temporal evolution of most features to eliminate spikes. It ignores completely the spatial information in the image sequence, treating each pixel as an independent time series of values. Spikes are identified by a strong second derivative component.

Here's a demonstration of a single frame (of a longer image sequence) being despiked by zspike. Note that the spikes blend in with the image data. Left: original data; center: bad-pixel mask from ZSPIKE; right: corrected image. To appreciate the algorithm better, you should view this MPEG movie (32 MB, sorry).

Cartography

When I wrote PDL::Transform, I really wanted just three things:

A way to co-register data with correct FITS pointing information (like the ZMATCH function I wrote in IDL so long ago);
An image warper for correlation tracking codes;
An image-to-map converter that would produce maps of the Sun given images from an arbitrary place (or vice versa).

The first two things are accomplished by PDL::Transform::t_fits and PDL::Transform::t_lookup, but mapping requires some care. Of course, image-to-map conversion is reminiscent of cartography in general, so over Christmas break 2002 I wrote PDL::Transform::Cartography as a sort of overkill for the third problem. The Cartography projections all convert (lon,lat) coordinates on a spherical surface to some sort of projected coordinate, and most of the standard map projections are there.

There's a vector coastline map of the world included, too. The method that returns it is PDL::Transform::Cartography::earth_coast.

Here are some sample command sequences and what they produce.

Commands Output
use PDL::Transform::Cartography; # load module $a = graticule(10,1)->glue(1,earth_coast); # Load coastline map $a->((2)) += (7-3*$a->((2)))*($a->((2))>0); # Make it more colorful $a->((2)) -= 3*($a->((2))<0); $ts = 'longitude '; $ps = 'latitude '; $as = "azimuth "; # Shorthand strings $ms = 'projection '; $rs = '(radians)'; $ds = '(degrees)';
[ None (setup) ]
A simple projection: Plate Careé (lon/lat)

@wopts = (size=>[4,2.5],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'plate_caree.gif/gif',@wopts); $w->lines($a->clean_lines,{title=>"Plate Caree $ms",xt=>$ts.$ds,yt=>$ps.$ds}); $w->close;

Mercator projection

@wopts = (size=>[4,3],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'mercator.gif/gif',@wopts); $opt = {title=>"Mercator $ms",xt=>$ts.$rs,yt=>"tan $ps"} $w->lines($a->apply(t_mercator)->clean_lines,$opt); $w->close;

Oblique Mercator
All the projections are easily rendered oblique by specifying a new origin. Here's an oblique mercator centered on the USA:

@wopts = (size=>[4,3],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'mercator-2.gif/gif',@wopts); $ps = "projected "; $opt = {title=>"Oblique Mercator $ms",xt=>"$ps theta $rs",yt=>"$ps tan phi"} $w->lines($a->apply(t_mercator(o=>[-95,40]))->clean_lines,$opt); $w->close;

Gnomonic map
Gnomonic maps are useful for showing great circle routes -- which project to straight lines. Here's one showing routes around North America and Western Europe. Notice that the meridians are straight -- that's because they're great circles themselves.
@wopts = (size=>[4,4],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'gnomonic.gif/gif',@wopts); $ps = "projected "; $opt = {title=>"Oblique Gnomonic $ms",xt=>"$ps X",yt=>"$ps Y"} $w->lines($a->apply(t_gnomonic(o=>[-30,45],c=>50)),$opt); $w->close;

View from space
This assumes an altitude of 0.1 Earth radii (380 miles), observer longitude 123 degrees west, 32 north (off the San Diego coast); the camera is pointed due north and 33 degrees down from horizontal. You can see the San Francisco Bay at the bottom of the picture; Puget Sound, the Canadian coast, and southern Alaska are visible at the top.

@wopts = (size=>[4,3],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'aerial.gif/gif',@wopts); $t = "Pacific Coast from ISS over San Diego"; $opt = {xr=>[-25,25],yr=>[-20,10],title=>$t,xt=>$ds,yt=>$ds}; $w->lines($a->apply(t_perspective(r0=>1.1,o=>[-123,32],c=>[0,-33])),$opt); $w->close;

Launch video
This isn't really a physical launch trajectory -- it just demonstrates aerial views from a changing perspective. The loop below made 300 frames which I mpegified. Click the image to view the MPEG.
`rm -rf launch; mkdir launch`; for $i(0..300) { my $w = pgwin(dev=>sprintf("launch/%3.3d.gif/gif",$i),@wopts); $zz = [-90+$i*$i/500.0, 60-$i*60/200.0]; $t = t_perspective(r0=>1+$i*0.01,fov=>80.0,o=>$zz); $rn =[-30,30]; $tt = sprintf("lon=%5.1f, lat=%5.1f, r=%5.1f",$zz->[0],$zz->[1],(1+$i*0.01)); $w->lines($a->apply($t),{xr=>$rn,yr=>$rn,xt=>$ds,yt=>$ds,title=>$tt}); $w->close; print "$i "; }

Many other transformations are present -- say "??Cartography" in perldl to find out which ones exist. They all can handle both vector data (as shown here) and image data (via the ->map method; see below).

Arbitrary transforms & image resampling

The PDL::Transform module (added to CVS in Nov 2002) implements arbitrary coordinate transforms. A large number of simple transformations, such as linear scaling, offsets, and FITS pixel-to-scientific conversions, are predefined. You can string together a bunch of transformations using a compose() operator.

The main purpose of PDL::Transform is to remap images. The impetus to develop it was for some image regularization and distortion work that was needed for my stereoscopic magnetograph, but (in keeping with the perl ideal of ambitious laziness), I implemented a fully general module that can do general coordinate transformations in N dimensions. It was at least partially inspired by Dominic Zarro's plot_map image-morphing package for IDL.

The centerpiece of PDL::Transform is the arbitrary-grid resampling code in the map method. The resampling code can do rapid sampling, interpolation, or fully anti-aliased variable-grid resampling of input data. The anti-aliasing is especially important for transformations that shrink part of the image.

For example, some years ago I used a custom piece of code to do conformal 2-D radial expansion on SOHO/LASCO C-3 images, and image very faint structures some 29 solar radii away from the Sun (DeForest, Plunkett, & Andrews 2001). That work was possible because we were able to average together huge numbers of pixels where the radial transformation shrunk the image (far from the Sun).

PDL::Transform::map uses tracking of the Jacobian of the transformation to average over as much of the input data as is appropriate, at each part of the output image: it implements true variable-sized spatial filtering, based only on the measured Jacobian of the arbitrary transformation you feed in.

Here is a demonstration using the M51 galaxy and a radial coordinate transform similar to the one we used in 2001.

Orginal image
This is a 256x256 FITS image of the M51 galaxy that is distributed with PDL itself.
The images below were produced with the following snippet of PDL code, using periodic boundary conditions on the original image (to demonstrate the effects of reduction as well as enlargement).
use PDL::Transform; $a = rfits('m51.fits'); $t = t_compose(t_scale(256/3.14159/2),t_radial(o=>[128,128],r0=>3)); $w = pgwin(device=>'m51.gif/gif',size=>[4,4],HardLW=>1); @opt = (XTitle=>'angle',YTitle=>'log(radius)'); $w->imag( $t->map($a,{m=>'sample',b=>'p'}) , {title=>'Sampled',@opt}); $w->imag( $t->map($a,{m=>'linear',b=>'p'}) , {title=>'Interpolated',@opt}); $w->imag( $t->map($a,{m=>'jacob', b=>'p'}) , {title=>'Anti-aliased',@opt}); $w->close;

Sampling is the fastest method -- this distortion was accomplished in less than 300ms by a 1GHz Athlon workstation. The original pixels are visible in the lower part of the image, and aliasing is present in the upper part. Interpolation gets rid of the original pixel images in the lower part of the image, but retains the aliasing in the upper part. Most of the original pixels are ignored because they fall between the centers of the inverse-transformed pixels of the output plane. This distortion required ~600ms. Anti-aliasing uses variable-sized spatial filtering to eliminate artifacts from the transformed image. It is slower than the other two methods, because each input pixel is transformed several times and the Jacobian is inverted at each output pixel. This distortion required about 10 seconds.

Movie: this MPEG (5.4 MB) shows the effects of aliasing a little more dramatically. It results from looping the above code over position: the origin moves vertically across the image and through the center of the galaxy, changing the pixel map enough to generate twinkling as the aliasing function changes.

Range operator

Indexing

range is a rank-reducing index operator, a generalization of index and indexND. Like those operators, it has two-way dataflow: you can get values out of an array by using the result of range in an expression, or you can put them back in to the original array via elementwise assignment to the result.

You start with a source array S (with N dimensions of order S₀...S_N-1) and an index array I with M dimensions of order I₀...I_M-1. If I₀ == N, then you can make the expression:

  $S->range($I)

and it will have dimensionality (I₁..I_M-1). range treats $I as a collection of column-vectors that index elements of $S.

So far, range sounds just like the indexND operator. But unlike indexND, range is fully threadable. If I₀ < N, then whole columns of $S are extracted, so $S->range($I) has dimensionality (I₁..I_M-1,S_M..S_N-1).

Ranging

Furthermore, range's name comes from the fact that it can extract N-D rectangles out of $S. This aspect is handy for (e.g.) vectorizing operations on disparate regions of a single data set. You feed in an additional size vector, like so:

  $S->range($I,$D);

$D must be a scalar or an M-vector; scalar values are treated as if they are repeated over all (M-1) non-collapsed dimensions of the index. The output is then a list of M-cubes extracted from $S. In that case the output has dimensionalty ( I₀...I_M-1,D₀...D_M-1,S_M...S_N-1). Those middle dimensions range over elements within the little pieces that you extracted. The beginning dimensions range over the different pieces specified by $I, and the ending dimensions range over the dimensions of $S that you didn't index.

Notice that you're free to group the dimensions however you want -- you can treat the output as a doubly-threaded collection of M-1-dimensional cubes, or as a single-threaded collection of N-dimensional cubes, because the extra $S dimensions are tacked on the end of the indexed dims.

Mix-and-match

Sometimes you want to index single hyperplanes out of a single dimension of a piddle, but get ranges in others. You can do that with range by setting elements of $D to 0. That normally wouldn't mean anything useful, because the number of elements in the output is the product of all the orders of all of its dimensions. So if you set a particular order to 0, it's instead treated as a "silent" dimension of order 1 -- you get a single sample at the corresponding index coordinate, but the dimension is omitted from the final dimension list.

Silent dimensions may seem strange -- but the limiting case where all the range dimensions are silent is just the indexing case at the top of this little screed. So all that they mean is that you can mix-and-match which dimensions you use to index and which you use to get ranges.

Boundary conditions

Indexing is a pain to get right in the general case because you have to keep track of the boundaries. Affine transformations like slice are very fast but don't let you run off the end of the source array. range takes a performance hit anyway by not being affine, so it allows you to specify what to do if you index off the end of a boundary. You access boundary conditions like so:

  $S->range($I,$D,$B)

where $B is a numeric or string code for whatever boundary condition you want to use. (See the documentation for all the abbreviations). Allowed conditions are:

`forbid': Behaves like an affine transformation: if you attempt to index out-of-bounds, you throw an exception. This is the default behavior.
`extend': Just duplicates the boundary elements of the array out to infinity.
`truncate': Inserts a virtual pointer to a missing-value element into the array, effectively duplicating BAD or 0 out to infinity.
`periodic': Calculates coordinates modulo the dimensions of $S, so that $S is effectively duplicated to fill all of Z^N

Examples

range is very powerful for general N-D operations, and is extensively used in PDL::Transform::map (see above). Several trivial examples of its use exist in the PDL documentation; I'll write some here soon.

Commands	Output
use PDL::Transform::Cartography; # load module $a = graticule(10,1)->glue(1,earth_coast); # Load coastline map $a->((2)) += (7-3$a->((2)))($a->((2))>0); # Make it more colorful $a->((2)) -= 3*($a->((2))<0); $ts = 'longitude '; $ps = 'latitude '; $as = "azimuth "; # Shorthand strings $ms = 'projection '; $rs = '(radians)'; $ds = '(degrees)';	[ None (setup) ]
A simple projection: Plate Careé (lon/lat) @wopts = (size=>[4,2.5],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'plate_caree.gif/gif',@wopts); $w->lines($a->clean_lines,{title=>"Plate Caree $ms",xt=>$ts.$ds,yt=>$ps.$ds}); $w->close;
Mercator projection @wopts = (size=>[4,3],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'mercator.gif/gif',@wopts); $opt = {title=>"Mercator $ms",xt=>$ts.$rs,yt=>"tan $ps"} $w->lines($a->apply(t_mercator)->clean_lines,$opt); $w->close;
Oblique Mercator All the projections are easily rendered oblique by specifying a new origin. Here's an oblique mercator centered on the USA: @wopts = (size=>[4,3],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'mercator-2.gif/gif',@wopts); $ps = "projected "; $opt = {title=>"Oblique Mercator $ms",xt=>"$ps theta $rs",yt=>"$ps tan phi"} $w->lines($a->apply(t_mercator(o=>[-95,40]))->clean_lines,$opt); $w->close;
Gnomonic map Gnomonic maps are useful for showing great circle routes -- which project to straight lines. Here's one showing routes around North America and Western Europe. Notice that the meridians are straight -- that's because they're great circles themselves. @wopts = (size=>[4,4],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'gnomonic.gif/gif',@wopts); $ps = "projected "; $opt = {title=>"Oblique Gnomonic $ms",xt=>"$ps X",yt=>"$ps Y"} $w->lines($a->apply(t_gnomonic(o=>[-30,45],c=>50)),$opt); $w->close;
View from space This assumes an altitude of 0.1 Earth radii (380 miles), observer longitude 123 degrees west, 32 north (off the San Diego coast); the camera is pointed due north and 33 degrees down from horizontal. You can see the San Francisco Bay at the bottom of the picture; Puget Sound, the Canadian coast, and southern Alaska are visible at the top. @wopts = (size=>[4,3],j=>1,hardlw=>1,hardch=>2); # Set up window options $w = pgwin(dev=>'aerial.gif/gif',@wopts); $t = "Pacific Coast from ISS over San Diego"; $opt = {xr=>[-25,25],yr=>[-20,10],title=>$t,xt=>$ds,yt=>$ds}; $w->lines($a->apply(t_perspective(r0=>1.1,o=>[-123,32],c=>[0,-33])),$opt); $w->close;
Launch video This isn't really a physical launch trajectory -- it just demonstrates aerial views from a changing perspective. The loop below made 300 frames which I mpegified. Click the image to view the MPEG. `rm -rf launch; mkdir launch`; for $i(0..300) { my $w = pgwin(dev=>sprintf("launch/%3.3d.gif/gif",$i),@wopts); $zz = [-90+$i$i/500.0, 60-$i60/200.0]; $t = t_perspective(r0=>1+$i0.01,fov=>80.0,o=>$zz); $rn =[-30,30]; $tt = sprintf("lon=%5.1f, lat=%5.1f, r=%5.1f",$zz->[0],$zz->[1],(1+$i0.01)); $w->lines($a->apply($t),{xr=>$rn,yr=>$rn,xt=>$ds,yt=>$ds,title=>$tt}); $w->close; print "$i "; }


Sampling is the fastest method -- this distortion was accomplished in less than 300ms by a 1GHz Athlon workstation. The original pixels are visible in the lower part of the image, and aliasing is present in the upper part.	Interpolation gets rid of the original pixel images in the lower part of the image, but retains the aliasing in the upper part. Most of the original pixels are ignored because they fall between the centers of the inverse-transformed pixels of the output plane. This distortion required ~600ms.	Anti-aliasing uses variable-sized spatial filtering to eliminate artifacts from the transformed image. It is slower than the other two methods, because each input pixel is transformed several times and the Jacobian is inverted at each output pixel. This distortion required about 10 seconds.