LABTITLE=: 'Fractals, Visualization & J'
LABCOMMENTS=: 'Based on Cliff Reiter''s book Fractals Visualization and J'
LABFOCUS=: 0

NB. =========================================================
Lab Chapter Introduction
NB. =========================================================
Lab Section Introduction
This lab is designed to provide a hands-on introduction to the visualization of mathematical objects, algorithms and data. This includes the plotting of functions in two and three dimensions, the creation of fractals using a wide variety of methods and the manipulation of data and images in general.

The lab is based on the first 2 chapters of "Fractals, Visualization & J" by Cliff Reiter.

The lab makes essential use of the standard graph utilities. If necessary, first take the lab "Graph Utilities".
) 

NB. =========================================================
Lab Section
We require the graph, numeric and trig libraries, which can be loaded as follows:
)
load 'graph numeric trig'

NB. =========================================================
Lab Section Terminology
Some terminology:

"object graphics" are drawn with the objects supported by the Windows graphic interface (GDI), for example, lines and polygons. In J, these can be drawn directly using the graph package.

"raster graphics" are drawn by plotting suitable pixels. In J, these are drawn as bitmaps.

This chapter has examples of both types: drawing a polygon using object graphics, and a Sierpinksi triangle using raster graphics.
)
PREPARE
coinsert 'jgl2'
PREPARE

NB. =========================================================
Lab Section
The lab will show drawings in one or more windows. You can rearrange and resize these windows as required.

The windows can be left open as you step through the lab. However, you can also close them if you wish - they will be re-created when they are next used.
)

NB. =========================================================
Lab Section
You are encouraged to experiment as you go through the course.

You can modify and re-enter examples by moving the cursor up to a line in the session manager, and pressing enter.

For example, the next two sections draw a pentagon. Try modifying the example to draw a hexagon.
)

NB. =========================================================
Lab Section Experiment: Plotting Polygons
We start with an example of object graphics: a drawing of a 5-sided regular polygon, i.e. a pentagon. First generate 5 points equally spaced on a unit circle, and draw a polygon from them.

The calculations are shown below. The next section draws the polygon.
)
NB. divide 2 pi into five intervals:
2p1*(i.5)%5

NB. generate the points as cosine-sine pairs:
+. r. 2p1*(i.5)%5

NB. =========================================================
Lab Section
Draw the polygon.

The commands used open a graph window, draw the polygon, then show the window.

Note that the right argument to gdpolyon should be a list of points. A matrix would be treated as a separate polygon drawing for each row.

Now try drawing a blue hexagon (6 sides) before continuing the lab.
)
gdopen ''
RED gdpolygon ,+. r. 2p1*(i.5)%5
gdshow''

NB. =========================================================
Lab Section
The same technique can be used for regular polygons of any number of sides.

The utility int01 can be used to generate the intervals. This produces an extra value at the end, but this does not affect the drawing.

Try varying the 2p1 multiplier to produce different shapes.
)
((i.5)%5) ; int01 5

gdopen ''
BLUE gdpolygon ,+. r. 6p1 * int01 11
gdshow''

NB. =========================================================
Lab Section Simplify Code
We can simplify the code a little here.

First, the utility "polygon" can be used to generate the required values for a regular polygon. The optional left argument is a scale factor on the multiplier 2p1. The result of polygon is cleaned, i.e. values near 0 are set to 0.

For example:
)
(+. r. 6p1 * int01 11) ; 3 polygon 11

NB. =========================================================
Lab Section
Next, the sequence:

   gdopen''
   gdpolygon data
   gdshow''

can be written in a single line, as:

   gdpolygon gddraw data

Here, gddraw is an adverb that wraps its left argument between gdshow and gdopen.
)

NB. =========================================================
Lab Section
For example, the line of code run below is equivalent to:

  gdopen''
  TEAL gdpolygon ,7 polygon 15
  gdshow''
)
TEAL gdpolygon gddraw ,7 polygon 15

NB. =========================================================
Lab Section
This lab will, for the most part, employ the short form using gddraw. If you have trouble following the code, you should re-enter it using the steps gdopen, gdpolygon... and gdshow''
)

NB. =========================================================
Lab Section Graph Windows
The windows used so far are named "graph".

This lab sometimes generates windows named "lab", as in the next section. The idea is that you can experiment with the "graph" window while the "lab" window remains open.
)

NB. =========================================================
Lab Section
See if you can reproduce the "lab" figure.
)
PREPARE
gdopen 'fvjlab';'lab'
YELLOW gdpolygon ,2 polygon 7
GREEN gdpolygon 0.62 * ,3 polygon 7
gdshow''
PREPARE

NB. =========================================================
Lab Section Sierpinksi Triangle
Now we look at an example of raster graphics.

A Sierpinski triangle can be generated by repeated copying.

Given a matrix m1, try to find a J expression that will create another matrix, m2, of the form:

m1 m1
m1  0

For example, if m1 is:

 1 1
 1 0

then m2 should be:

 1 1 1 1
 1 0 1 0
 1 1 0 0
 1 0 0 0
)
NB. define m1:
[m1=: 1 1 ,: 1 0

NB. =========================================================
Lab Section
A simple verb to do this is:
)
f=: ,~,.~

f m1

f ,1        NB. generates m1

f f ,1

f ^:3 ,1

NB. =========================================================
Lab Section
Try entering:

   f ^:6 ,1

You can browse this matrix using:

   wdview f ^:6 ,1
)

NB. =========================================================
Lab Section
The viewmat utility can be used to display the matrix.
)

NB. =========================================================
Lab Section
)
PREPARE
f=: ,~,.~
PREPARE
viewmat f ^: 4 ,1

NB. =========================================================
Lab Section
Exeriment with different powers of f, for example, try:

  viewmat f ^: 6 ,1

Each iteration increases the number of points 4 fold, so it is not advisable to try much higher powers.

A left argument of a color matrix can be used to set colors. Try:

   (TEAL ,: RED) viewmat f ^: 6 ,1
)
PREPARE
f=: ,~,.~
PREPARE

NB. =========================================================
Lab Section Sierpinski Carpet
The same principle can be used to construct the Sierpinski carpet.

Can you duplicate the "lab" example?
)
PREPARE
f=: ((],,),.~),.],,~
(BLUE,:YELLOW) viewmat (f ^: 5 ,: 1);'blab';'lab'
PREPARE

NB. =========================================================
Lab Chapter Fractals and Transformations
NB. =========================================================
Lab Section Introduction
In this chapter we look at organizing and applying the mathematics required for creating transformations of the plane: rotations, rescaling and translations. This is done in a uniform way by using homogeneous coordinates.

The work on transformations prepares the way for a deeper look at fractals. In particular, iterated function systems and some recursive refinement methods are described.
)
PREPARE
load'numeric trig graph'
PREPARE

NB. =========================================================
Lab Section
Note that the standard graphics box is a square from _1 _1 to 1 1. The objects defined in the lab fit in this box, but might not do so after being transformed.

If an object cannot be seen when you display it, check its coordinates. If necessary, scale down the object to fit into the standard box.
)

NB. =========================================================
Lab Section Plotting Functions and Parametric Curves
Now lets experiment with plotting functions.

Define x to be 1000 points in the range from _1 to 1, then draw x against sin pi times x.
)
steps _1 1 10
x=: steps _1 1 1000
y=: sin o. x
gdlines gddraw ,x ,. y

NB. =========================================================
Lab Section
In the previous example, we specified x, and then defined y as a function of x.

In a parametric plot, both x and y are functions of some variable, typically denoted t.

In the example, a factor of 0.9 shrinks the plot to fit comfortably in the graph box:
)
t=: steps 0 2p1 1000
x=: cos t
y=: sin t
gdlines gddraw ,0.9 * x,.y

NB. =========================================================
Lab Section
Here is an equivalent expression that makes it easier to think in terms of the two functions being graphed:
)
PREPARE
t=: steps 0 2p1 1000
PREPARE
gdlines gddraw ,0.9 * +. r. t

NB. =========================================================
Lab Section
Modify these functions for some interesting experiments.

For example, try:

  gdlines gddraw ,0.9 * (cos ,. sin*cos) t
  gdlines gddraw ,0.9 * (cos@(3&*) ,. sin) t
  gdlines gddraw ,0.9 * (cos@(3&*) ,. sin@(2&*)) t
)
PREPARE
t=: steps 0 2p1 1000
PREPARE

NB. =========================================================
Lab Section
Can you reproduce the "lab" example?
)
PREPARE
t=: steps 0 2p1 1000
gdopen 'lab'
gdlines ,0.9 * (cos@(5&*) ,. sin@(8&*)) t
gdshow''
PREPARE

NB. =========================================================
Lab Section
Note that all plots are parametric.

For example, here is the first plot written in parametric style:
)
t=: steps _1 1 1000
gdlines gddraw ,(] ,. sin@o.) t

NB. =========================================================
Lab Section
swap the x and y functions:
)
PREPARE
t=: steps _1 1 1000
PREPARE
gdlines gddraw ,(sin@o. ,. ]) t

NB. =========================================================
Lab Section The Koch Snowflake
The Koch curve and snowflake are curves that are described by iterative refinement. When this refinement is applied to an initial triangle, the Koch snowflake results.

The basic step in the iterative process is to take an edge and replace it with four edges introducing a bulge. If the given edge is trisected, then the two end thirds are retained and the center third is replaced by the other two edges of an equilateral triangle placed on the center third.
)

NB. =========================================================
Lab Section
As a specific example of refining an edge consider an edge from 0 0 to 1 0 which can be specified by the rows of the matrix, each row specifying one vertex (endpoint):

   0 0
   1 0
)

NB. =========================================================
Lab Section
Applying the basic step to this edge would result in the following vertex sequence (polygon) where y is -%:%12 :

   0    0
   1r3  0
   1r2  y
   2r3  0
   1    0
)

NB. =========================================================
Lab Section
In order to implement this subdivision we define functions that calculate the midpoint and the points which divide the segment into thirds.
)
tri1=: 2r3&*@[ + 1r3&*@]           NB. one third point
tri2=: 1r3&*@[ + 2r3&*@]           NB. two thirds point
mid=:  -:@+                        NB. midpoint
0 0 tri1 1 0
0 0 tri2 1 0
0 0 mid  1 0
|. 5 6 7 8                         NB. reverse items
nor=: _1 1&*@|.                    NB. left normal vector
nor 1 2
bulge=: mid + (%:%12)&*@nor@:-     NB. bulge point
0 0 bulge 1 0

NB. =========================================================
Lab Section
Note that bulge uses the at composition @:, so that nor@:- finds the normal to the vector difference, not the normal to each element in the difference.

Compare:
)
nor@- 1 2

nor@:- 1 2

NB. =========================================================
Lab Section
Thus, bulge computes the midpoint with a displacement in the opposite of the left normal direction. Below, segdiv divides up the given segment into pieces.

Note that the last point is not included in the result since it will be included as the first point in the next segment to be refined.
)
segdiv=: [ , tri1 , bulge ,: tri2
0 0 segdiv  1 0

NB. =========================================================
Lab Section
The function refine applies the division to neighboring pairs of vertices.
)
i. 4 3
1|.i. 4 3              NB. rotate items by one position cyclically
refine=: ,/@(] segdiv"1 (1&|.))

NB. =========================================================
Lab Section
Generate a snowflake with a single iteration of refine.

In this lab, snow is usually red.
)
[p=: |:2 1 o./ 2r3p1*i.3   NB. make a triangle
RED gdpolygon gddraw ,refine ^:1 p

NB. =========================================================
Lab Section
Generate a snowflake with four iterations of refine.
)
PREPARE
p=: |:2 1 o./ 2r3p1*i.3   NB. make a triangle
PREPARE
RED gdpolygon gddraw ,refine ^:4 p

NB. =========================================================
Lab Section
An experiment with superimposing the snowflakes.

Note that arguments to gdpolygon can be boxed lists.
)
PREPARE
p=: |:2 1 o./ 2r3p1*i.3   NB. make a triangle
PREPARE
(BLUE;RED) gdpolygon gddraw (,refine ^:4 p);,refine ^:1 p

NB. =========================================================
Lab Section
Can you reproduce the lab figure?
)
PREPARE
gdopen 'fvjlab';'lab'
RED gdpolygon , refine ^:4 [ 3 polygon 4
gdshow''
PREPARE

NB. =========================================================
Lab Section Transformations and Homogeneous Coordinates
The dot product of two vectors is computed by taking the sum of the pairwise products of the elements of two given vectors. The standard matrix product generalizes this so that the result is the dot product of each row of the left argument with each column of the right argument.

In J, the matrix and dot products are both denoted by the inner product:
  +/ . *
)
1 _1 2 +/ . * 0 2 4     NB. dot product

[a=: i.3 2
[b=: 0 1,:1 2
a +/ . * b              NB. matrix product

NB. =========================================================
Lab Section Transformations
Matrices can be used to represent transformations of the plane.

Suppose p is a vector giving the x-y coordinates of a point in the plane. If we want to expand or contract in the x-direction by a factor of r and in the y-direction by a factor of s, then multiplying p on the right by a matrix of the following form will effect the transformation.

  r  0
  0  s

For example, the matrix a rescales the x-direction by a factor of 3.
)
[a=: 2 2$3 0 0 1
[p=: 10 10
p +/ . * a

NB. =========================================================
Lab Section Rotations
If we want to rotate a point p through an angle of t counterclockwise about the origin, we can multiply by a matrix of the following form.

   cos(t)  sin(t)
  -sin(t)  cos(t)

Let nor be defined as before, then we can construct a rotation matrix a as follows.
)
nor=: _1 1&*@|.

[a=: (,: nor) 2 1 o. 1r6p1     NB. rotate by 1r6p1

1 0 +/ . * a

p=: 10 10

p +/ . * a

NB. =========================================================
Lab Section Translations
Lastly, if we want to translate, we can just add.

Thus, to move a point, p, 1 unit to the right and 2 units upward, add:

  p + 1 2
)

NB. =========================================================
Lab Section
One often wants to think about transformations in the plane as a sequence of rotations, rescalings and translations. It is somewhat awkward to do this using the above ideas since the compositions will be a matrix product plus some translation vector.

In order to help organize the desired compositions, we use a slightly different coordinate system that allows all three types of transformations to be implemented with matrix multiplication. These are called homogeneous coordinates.
)

NB. =========================================================
Lab Section Homogeneous Coordinates
Basically, we simply attach a trailing 1 to each x-y pair. More precisely, we consider two triples to represent the homogeneous coordinates of the same point in the plane if they are different by a nonzero scalar multiple, so 3 2 1 and 6 4 2 denote the same point in homogeneous coordinates.

We will always scale homogeneous coordinates so there is a trailing 1.

Now to rescale by r and s, rotate by t, or translate by u,v we multiply p,1 by a matrix of one of the following forms.

 rescaling     rotation                translation
 r  0  0       cos(t)  sin(t)  0       1  0  0
 0  s  0      -sin(t)  cos(t)  0       0  1  0
 0  0  1         0       0     1       u  v  1
)

NB. =========================================================
Lab Section
For example, if transformation matrices a1, a2, and a3 are defined as below, then the result of applying them to the point 1 1 after it is converted to homogeneous coordinates is shown.
)
NB. matrix a1 rescales x by 3
[a1=: 3 0 0, 0 1 0,: 0 0 1
10 10 1 +/ . * a1

NB. matrix a2 rotates by 45
[a2=: ((-:%:2) * 1 1 0,:_1 1 0),0 0 1
10 10 1 +/ . * a2

NB. matrix a3 translates by 1 2
[a3=: 1 0 0, 0 1 0,: 1 2 1
10 10 1 +/ . * a3

NB. =========================================================
Lab Section Experiment: Some Transformations
Begin by defining a matrix p whose rows give the homogeneous coordinates of an L-shaped polygon.
)
[p=: (0.25 * 6 2$0 0 2 0 2 1 1 1 1 3 0 3) ,.1

NB. =========================================================
Lab Section
Display the L-shape
)
PREPARE
[p=: (0.25 * 6 2$0 0 2 0 2 1 1 1 1 3 0 3) ,.1
PREPARE
RED gdpolygon gddraw ,}:"1 p

NB. =========================================================
Lab Section
Now let a1 be a matrix and t1 the transformation associated with matrix multiplication by a1 on the right.

The next draws t1 p.

How are the polygons p and t1 p related?
)
[a1=: 0.5 0 0,0 0.5 0,:0 0 1
t1=: +/ .* & a1

NB. =========================================================
Lab Section
Draw p and t1 p.
)
PREPARE
[p=: (0.25 * 6 2$0 0 2 0 2 1 1 1 1 3 0 3) ,.1
[a1=: 0.5 0 0,0 0.5 0,:0 0 1
t1=: +/ .* & a1
PREPARE
(RED;BLUE) gdpolygon gddraw (,}:"1 p);(,}:"1 t1 p)

NB. =========================================================
Lab Section
Experiment with the following transformation matrices.
)
[a2=: _1 0 0, 0 _1 0,: 0 0 1
[a3=: 2 0 0, 0 1 0,: 0 0 1
[a4=: 0.5 0 0, 0 _0.5 0,: 0 0 1
[a5=: 1 0 0, 0 1 0,: _1 _1 1

t2=: +/ .* & a2
t3=: +/ .* & a3
t4=: +/ .* & a4
t5=: +/ .* & a5

NB. =========================================================
Lab Section
For example, try:

   (RED;BLUE) gdpolygon gddraw (,}:"1 p);(,}:"1 t2 p)
)

NB. =========================================================
Lab Section Transformation Exercises
Let t6 be the transformation resulting from a4 +/ . * a5.

Display t6 p and compare the result of t4@t5 p with t6 p.

How do these compare to the result of the transformation, t7, associated with a5 +/ . * a4?
)
t6=: +/ . * & (a4 +/ . * a5)
t7=: +/ . * & (a5 +/ . * a4)

NB. =========================================================
Lab Section
Define a function rot such that rot y gives the rotation matrix (in homogeneous coordinates) of a rotation through an angle y.

Create a transformation t8 that rotates by an angle of p/3.
Apply it to p and view the result.

Compare t8@t5 and t5@t8.
)

NB. =========================================================
Lab Section
Next consider transformations associated with the identity matrix and inverse matrices.

Plot and view the results of those transformations on p.
)
PREPARE
a3=: 2 0 0, 0 1 0,: 0 0 1
a5=: 1 0 0, 0 1 0,: _1 _1 1
a8=:  0 0 1 ,~ ((cos,sin),:-@sin,cos) 1p1%3
PREPARE
]a9=: =i.3             NB. 3 by 3 identity matrix
]a10=: %. a3           NB. matrix inverse to a3
]a11=: %. a5           NB. another inverse
]a12=: %. a8           NB. another inverse

NB. =========================================================
Lab Section Experiment: Some Collages of Transformations
Define p to be a polygon in homogeneous coordinates giving the unit square with one vertex at the origin.

Define three transformations, t0, t1, t2 as follows:
)
[p=: (#:0 1 3 2),.1

[a0=: 0.5 0 0, 0 0.5 0,: 0 0 1
[a1=: 0.5 0 0, 0 0.5 0,: 0.5 0 1
[a2=: 0.5 0 0, 0 0.5 0,: 0 0.5 1

t0=: +/ . * & a0
t1=: +/ . * & a1
t2=: +/ . * & a2

NB. =========================================================
Lab Section
We want to create a collage of the transformations on the polygon. That is, we want the list of the three polygons that result from applying t0, t1 and t2 to p.

This is easy to implement using gerunds and gerund evaluate as below.
)
collage=: t0`t1`t2`:0

collage p

NB. =========================================================
Lab Section
Experiment with powers of collage, for example:

  collage ^: 2 p

Now draw each matrix of the result.
)

NB. =========================================================
Lab Section
For example:
)
RED gdpolygon01 gddraw ,"2 }:"1 collage ^: 6 p

NB. =========================================================
Lab Section
Repeat the experiment with the same transformations but use a five-pointed star (or some other polygon) as the starting polygon.

How does this affect the long term pattern?
)

NB. =========================================================
Lab Section Experiment: Iteration of Random Transformations
It turns out that object graphics cannot handle many thousands of objects. Thus when creating an image with many thousands of pixels, a pixel at a time, we do well to avoid object based graphics. Since arrays are easy to create in J we will create a large array giving a virtual screen. This can be viewed as a bitmap.
)

NB. =========================================================
Lab Section
Again define three transformations, t0, t1, t2 corresponding to the matrices you used in Section 2.9 to create a finite version of the Sierpinski triangle.

Now ?@3: is a function that, when applied to anything, randomly picks a number from 0 1 2. The function rt, given below, randomly applies one of the transformations to the input data and rt^:(i.5) applies a random transformation to the given input 4 times and catenates the results to the matrix.
)
PREPARE
[a0=: 0.5 0 0, 0 0.5 0,: 0 0 1
[a1=: 0.5 0 0, 0 0.5 0,: 0.5 0 1
[a2=: 0.5 0 0, 0 0.5 0,: 0 0.5 1

t0=: +/ . * & a0
t1=: +/ . * & a1
t2=: +/ . * & a2
PREPARE
rt=: t0`t1`t2@.(?@3:)

rt 0.9 0.9 1

rt^:(i.5) 0.9 0.9 1

NB. =========================================================
Lab Section
To view the result of rt in a bitmap, scale the result to the size of the bitmap, round to nearest integer and use as indices to insert 1 into a matrix of zeroes. Here we use a bitmap of size 500x500.

For example:
)
a=: }:"1  rt^:(i.125) 0.9 0.9 1
viewmat |.|: 1 (;/roundint 499 * a) } 500 500$0

NB. =========================================================
Lab Section
The last expression to view a set of indices will be used repeatedly, so we define a verb viewndx for it:
)
PREPARE
a=: }:"1  rt^:(i.125) 0.9 0.9 1
PREPARE
viewndx=: 3 : 'viewmat |.|: 1 (;/roundint 499 * y.) } 500 500$0'
viewndx a

NB. =========================================================
Lab Section
No pattern is discernable with only 125 values.

However, crank this up and the fractal starts to appear.

Repeat with 10000 values.
)
viewndx }:"1  rt^:(i.1000) 0.9 0.9 1

NB. =========================================================
Lab Section
Next try a different set of transformations.

Modify the above (including rt) to handle the four transformations associated with the following matrices.
)
[a3=: 0 0 0, 0 0.16 0,: 0.25 0 1
[a4=: 0.85 _0.04 0, 0.04 0.85 0,: 0.0375 0.17 1
[a5=: 0.2 0.23 0, _0.26 0.22 0,: 0.2 0.1025 1
[a6=: _0.15 0.26 0, 0.28 0.24 0,: 0.2875 _0.021 1

t3=: +/ . * & a3
t4=: +/ . * & a4
t5=: +/ . * & a5
t6=: +/ . * & a6

rt=: t3`t4`t5`t6 @.(?@4:)

NB. =========================================================
Lab Section
)
viewndx }:"1  rt^:(i.20000) 0.5 0.5 1

NB. =========================================================
Lab Section
Notice that the image looks like a fern, but the tip is very faint.

Modify rt so that it uses the same transformations, but the transformation associated with a4 is repeated 6 times (so a random number up to 9 is required). You should get a much better image.
)

NB. =========================================================
Lab Section
For example:
)
rt=: t3`t4`t5`t6 @.(?@9: { 0 1 1 1 1 1 1 2 3"_)
viewndx }:"1  rt^:(i.20000) 0.5 0.5 1
