Thursday, March 16, 2017

Calculating Hilbert Curve Coordinates

You may have noticed from some of my previous blog posts that I've fallen down a rabbit hole learning about Gray codes and Hilbert curves.  They weren't just random topics.  I wrote those articles in order to understand the topic of this one, how to easily convert a Hilbert Index to Hilbert Coordinates via a method described in a paper by John Skilling.

Skilling, J. (2004), "Programming the Hilbert Curve", AIP Conference Proceedings.

Regrettably I don't yet understand why it works to a level that I'm satisfied with.  The paper describes the process and gives a summary of previous work by Butz and Lawder but doesn't really give a clear explanation of why the transform works.  It may be perfectly clear to someone who works in the field, but I don't get it... yet.

What I am capable of doing though is summarising the information here with some graphics to explain the concept to the next person that gets nerd sniped by this.  I'll start by recommending my last few blog posts as they will give you an idea as to my train of thought.

A basic concept to understand is how binary numbers are treated as coordinates.  In the animation below, the index is 6 bits long and counts up from 000000 to 111111.  This can be used to cover every coordinate on an 8 x 8 two dimensional grid if every second bit belongs to the coordinates of one dimension while the others belong to the other dimension.  In this case the red bits belong to the x axis while the blue bits belong to the y axis.  Everything I'll describe in this post is for two dimensions, but it works for any number of dimensions.  For example an 8x8x8 grid could be covered with a 9 bit index where bits 0, 3, 6 belong to one axis, 1, 4, 7 to another, and 2, 5 ,8 to another.

Binary Indexing Animation
Binary indexing in 2 dimensions
In the binary indexing animation you can see that there are diagonal lines.  These occur because more than one bit of the index can change at once affecting both dimensions.

The first step in the Skilling Transform is to take the binary reflected Gray code of the index.  This means that only one bit at a time can change.  Therefore the value of only one axis can change at a time leaving only vertical and horizontal movements.

Gray Code Animation
Gray code indexing in 2 dimensions

The addressing scheme detailed above is formally described in this section of the paper.  I put this here just so we're all on the same page about the variable names.

$$$p$$$ is the number of bits in each axis
$$$n$$$ is the number of dimensions

Addressing Scheme
Addressing description

Now onto the actual transform.  It's minimal and easy to compute.  It's not a recursive algorithm and takes a bit string that specifies the Hilbert index as an input and performs a series of swaps and inversions to return a bit string of the same length that describes the Hilbert coordinates.

Skilling transform

To see if I could understand the process better I wanted to visualise the transform.  The animation below show the transformation from Gray code indexing to the Hilbert curve.  You can see how the algorithm starts the transformation on the smaller features first and works up to larger features at each step.

Hilbert Curve Animation
Gray code to Hilbert curve via the Skilling transform in 2 dimensions.

Hopefully the worked example below will help people understand the process a little better.  The index is converted to its Gray code and the algorithm processes bits in a reverse order.  The grey boxes refer to the bit $$$r$$$ shown in the algorithm above.  The bits that are classified as low bits at each stage are underlined.
Hilbert Index to Coordinates
Worked example

I can't really offer much of a conclusion here.  The understanding I've gleaned so far is minimal and superficial at best.  I 'll leave this one for a while and come back to it in the future when I have more time.

I have a love hate relationship with problems like this.  I find them fascinating, but at the same time they remind me how much I miss formal study, and I get bummed out.  Kinda sucks having no one to bounce ideas off either.  Having said that though, I've finally generated some animations that I'm proud of.  You can find the undocumented scratch code here.

If you come across this page and have a better understanding of why this works I'd be glad to hear from you.

Saturday, March 4, 2017

Converting Binary to Gray Code with XOR

Over the last week I've been researching Gray Codes and how to generate them. I understand the basic concept but I wanted a deeper understanding.  Most of the explanations out there are a little dense so I'm going to try to make the available information a little bit more accessible to everyone. So let's start with a basic explanation of what a Gray code is.

A Gray Code is a special way of counting, where adjacent entries differ by only one digit. In practice this almost always refers to a binary counting sequence and more specifically a special code called the Binary Reflected Gray Code (BRGC). From here on, this is the Gray code that I'll be talking about.

These code have practical uses in many areas.  The one that I'm most familiar with is in rotary encoders where it's advantageous that only one digit changes at a time.  In a normal binary counting sequence you may have multiple bits change at once.  For example from $$$0111$$$ to $$$1000$$$. The problem with this is that if all the bits don't change at exactly the same moment, the value $$$1010$$$ may be read from the sensor.  If Gray codes are used, as you progress through the series only one bit changes at a time.  This leaves no ambiguity as to what the value is.

Let's start with a formal definition of a Gray code and try to explain things in easier to understand terms.

$$G_0=\epsilon \\ G_{k+1}=0G_k,1G_k^R$$

The series for zero bits is shown as $$$G_0=\epsilon$$$ and indicates an empty series.  The next equation, $$$G_{k+1}=0G_k,1G_k^R$$$ indicates that the series for a Gray code series of $$$k+1$$$ bits is created by taking the the Gray code of $$$k$$$ bits and joining a reversed copy to its end, with zeros added before entries in the first half of the series, and ones added before entries in the second half of the series.  This is how the series ensures that consecutive entries only differ by one digit.  This may be clearer in the diagram below.

Constructing $$$G_{k+1}$$$ from $$$G_{k}$$$

If $$$G_k$$$ is a Gray code where consecutive entries differ by one digit, then appending a mirror of itself to its end won't change that, except for where the two series join, as those two entries will be the same. Adding zero to the first half of the new series and one to the second half of the new series won't change it either but it does mean that the entries where the series join now differ by one digit. This means that the new series $$$G_{k+1}$$$ is now also a Gray code.  It also proves that the entries in the output are unique.  Starting with $$$G_0$$$, we know that all its entries are unique.  If we then assume that all the elements of $$$G_k$$$ are unique, then the elements of $$$G_{k+1}$$$ must also be unique as $$$G_k$$$ with zeros added before it and a reverse copy of $$$G_k$$$ with ones added before it create a new unique series. It's then provable by induction that the each entry in a Gray code is unique.

This gives us the instructions needed to construct a Gray code of any length.  Start with the minimal Gray Code that contains one empty element and keep doubling its length when adding a new digit.  

Gray Code
Constructing $$$G_3$$$ from $$$G_2$$$
Getting a picture in your mind of what a Gray code looks like will be helpful for the next step. For some situations the above representation can help, but in others, the vertical table representation with the index $$$n$$$ of each element is better.
Gray Code
Gray code indices in decimal

Now that I've explained the basics of what a Gray code is, I want to detail a proof of a well known formula to generate them with an xor operation.

$$H[n] = n \oplus (n>>1)$$

Where $$$H[n]$$$ is the $$$n$$$th element of the Gray code. $$$H$$$ is used at the moment as it's a hypothesis that this formula is equivalent to $$$G[n]$$$ described above.

I found a well thought out explanation of this proof in a Quora post but it could do with some further explanation.  So I'll repeat the proof here with some commentary to help those new to the subject. To be clear all of these operations are performed on the binary representation of the numbers involved.

First, let's define a helper function.

$$F[n,k]=(2^k-1) \oplus n$$

$$$F[n,k]$$$ is a function that inverts all the bits of a binary string $$$n$$$ of length $$$k$$$. $$$2^k$$$ is a binary string of length $$$k+1$$$ that is all zeros except for the most significant bit. Subtracting one from this will create a binary string of length $k$ that is all ones.  When the string $n$ is xored with this all of the bits are inverted. This is a function that essentially reverses the direction of  a binary count. For example, when $$$k=3$$$, and for a value of $$$n=001$$$, $$$F[001,3] = 110$$$.  Where $$$001$$$ is the second entry when counting in binary, $$$110$$$ is the second last.  It's a mathematical representation of reflecting/mirroring a binary count.

Now we have that out of the way we will prove the following

$$H[x] \oplus H[F[x,k]]=2^{k-1}$$

What does that mean?  It's saying that the output of the hypothesis function xored with the hypothesis function with the input ordering reversed will give a string of length $$$k$$$ that is all zeros except for the most significant bit.  This is a fancy way of saying that the output of $$$H$$$ is symmetrically except for the most significant bit.  This is a property of the Gray code that we are trying to prove. From the beginning of the series each element of a Gray code is equal to the element the same distance from the end except for the most significant bit.

As $$$F[x,k]$$$ is an inverted version of $$$x$$$ the most significant bit of one or the other, but not both, will be equal to one. Let's assume that it's $$$x$$$. It doesn't matter which one is chosen, the proof works either way.  We can now write

$$x=2^{k-1} \oplus y \quad where \quad y<2^{k-1}$$

It then follows that $$$F[x,k] = F[y,k-1]$$$.  This says that inverting $$$x$$$ which we have defined to start with one, will now start with zero and is equal to the inverted value of the new shorter bit string $$$y$$$.  For example if $$$x=1100101$$$ then $$$y = 100101$$$.  So $$$F[1100101,7]=0011010$$$ and $$$F[100101,6]=011010$$$.

Now onto the proof.  It's all maths at this point, you can ignore conceptual ideas here and just work through the process.

\begin{align}H[x] \oplus H[F[x,k]]&=(2^{k-1} \oplus y) \oplus (2^{k-2} \oplus (y>>1)) \oplus H[F[y,k-1]] \\ &=2^{k-1} \oplus 2^{k-2} \oplus H[y] \oplus H[F[y,k-1]] \\ &= 2^{k-1} \oplus 2^{k-2} \oplus 2^{k-2} \\ &= 2^{k-1}\end{align}

Now we've proven this, for strings of length $$$k$$$ using strings of length $$$k-1$$$ it just needs to be proved for strings of length 1 ie $$$k=1$$$ for the inductive proof to be complete.

\begin{align} H[0] \oplus H[1] &= 0 \oplus 0 \oplus 1 \oplus 0 \\&=1\end{align}

All of that above just proves that $$$H[n] = n \oplus (n>>1)$$$ results in an output series that is mirrored with the most significant bit inverted.  Now we need to prove the following

$$H[2^k \oplus x] = 2^k \oplus H[F[x,k]] \quad where \quad x<2^k$$

This states that if we add a one to the input to the hypothesis function, it's equivalent to inverting the input and then adding a one to the output.  It's hard to see, but this describes the process of going from a Gray code of $$$k-1$$$ bits to a Gray code of $$$k$$$ bits.

\begin{align} H[2^k \oplus x] &= 2^k \oplus x \oplus 2^{k-1} \oplus (x>>1) \\ &= 2^k \oplus 2^{k-1} \oplus H[x] \\ &=2^k \oplus H[F[x,k]]\end{align}

As we did before, this needs to be proved for $$$k=1$$$ for the induction proof to work.

\begin{align}H[10 \oplus 1] &= 10 \oplus H[0] \\H[11] &= 10 \oplus 0 \oplus 0  \\ 11 \oplus 01 &= 10 \\ 10=10 \end{align}

OK, that was a confusing rabbit hole we just fell down, but it proves that $$$H[n]=G[n]$$$ for all $$$n$$$.  I'm not entirely sure we need both parts of that proof, but it can't hurt to prove it twice. 😅 

We can now be sure that $$$G[n] = n \oplus (n>>1)$$$  How simple is that!  A Gray code can be generated by xoring its index with itself but right shifted by one.  This means that each digit of the Gray code is the xor of two consecutive digits in the input.

Generating a Gray code from an index - input in green output in purple
In the image above it's kind of like differentiating the bit stream.  If two consecutive bits are different it generates a one in the output.  So if we go from index to Gray code with something like a differentiation you might think we can go the other way with integration.  You sure can.  I won't prove it but you can look it up.

Now all this seems a little simple and useless apart from the use case I gave at the start, but this goes into some deep concepts.  We're looking at mathematics over a finite field, in this case $$$GF(2)$$$. Now here's the really cool part.  If you use each digit of a Gray code as a coordinate and plot it, you visit each vertex of a square for a two bit Gray code.  For a three bit Gray code you visit each vertex of a cube travelling along the edges. More generally for a k bit code you travel along the edges of a k-dimensional hypercube visiting all vertices.  This is called a Hamilton cycle.

Sorry if I'm ranting but I find all this fascinating, but it hasn't solidified in my mind yet.

If you've made it this far you deserve a joke.

Q. What's the best way to visit all the vertices of a hypercube?
A. On a Hamilton bicycle.

Hey give me a break.  I didn't say it'd be funny.

In case the equations don't render properly, here is a PDF of this page.
A PDF of the original solution on Quora.

Wednesday, February 22, 2017

Hilbert Curve Generation With Lookup Tables

I've been researching something that requires a Hilbert curve and I thought I'd share how to generate the path and also move between the index and coordinates of points.

For those of you unfamiliar with Hilbert curves let me quickly explain what one is and why you would want to use one.  The curve covers every point in a square with side length 2^n by moving up, down, left and right, starting at one corner and ending at an adjacent corner.  You could accomplish the same requirements by just scanning back and forth across the rows, but the advantage the Hilbert curve has is that nearby points on the 2D grid are generally also near each other on the curve.  If you were to scan back and forth you get a lot of points that are directly beside each other in 2D space but very far apart on the curve.

So why would you want to use one?  They're great for turning 2D areas into 1D streams of data while maintaining locality.  This comes in handy in image processing.  Depending on what you want to do, this may make the processing a lot easier.

Generation of the curve can be done recursively by first selecting an initial shape type and then using tables of sub-types and locations to generate the next stage.  In the animation below the initial shape is type 1.  When you go to the sub-type table you see that type 1 is replaced with types 2, 1, 1, and 4 in that order.  The position of each sub-type is defined by the parent type.  As the parent is type 1 the locations of the sub-types from the location table are:

2 at (0, 0)
1 at (0, 1)
1 at (1, 1)
4 at (1, 0)

The location table can then be used to find the coordinates of the sub-types.  As binary coordinates are used, the x and y coordinates of the sub-type can be appended to its coordinates to find the new sub-coordinates.   For example, the coordinates of the type 4 sub-type at position (1, 0) above are (1, 1) (0, 1) (0, 0) (1, 0).  Appending these to the coordinates of (1, 0) you find the new sub coordinates (11, 01) (10, 01) (10, 00) (11, 00).  Trust me it's confusing at first but after a while it all makes complete sense.  I'm still trying to come up with appropriate terminology for the process.  These tables are from a paper

It's a little confusing at first but I recommend reading that to help understand the process.

Recursive Hilbert Curve Generation

What if you don't want to generate the full curve and just want to know where along the curve a point lays or where a point on the curve is in 2D space?  Let's use the below example to demonstrate this.

44th Point at location (111, 101)

The above image shows the 44th point of the third stage at (111, 101).  It's actually the 45th point but we start counting at  0.  We add an index row to the tables provided in the above paper to make the process easier.

Generation Tables

We start by converting the index to a binary number with 2 bits for every stage.  In this case it's the third stage, so we get 6 bits.  These bits are placed two at a time on rows under the index column.  We also define the initial shape as type 1.  From this, the sub-type and location for the first row can be found.  The sub-type becomes the type for the next row.  This process is repeated two more times until we have three locations.  To get the final location coordinates, just append the x and y coordinates downwards to get (111, 101).

Converting an Index to a location

To go from coordinate to index, split the location up by removing the highest bit from the x and y coordinate for each line.  As the start type is known, the index can be found.  This index is then used to find the sub-type for the next row.  This is repeated two more times.  The indices are then concatenated to get 101100 in binary which is 44.

Converting a location to an index

I hope this helped to explain the process.  Trust me, even if it's still confusing, having diagrams and a couple of hopefully easy to follow examples will help.

Saturday, February 11, 2017

Entry and Exit Points for Space Filling Paths on a Grid

I've been researching space filling curves recently and came across an interesting algorithm to generate curves for arbitrary areas. To be completely accurate, space filling curves cover the entire unit square in continuous space, whereas what I'm looking at is the discrete form that covers a grid. Below when a path or curve is mentioned, what is meant is a way to cover a rectangular area by moving one square at a time, either up, down, left, or right, without crossing over itself. The curve starts in one corner and ends in another.

Being almost completely new to this, I set about understanding the description of the algorithm. Although it was accompanied by code, details of the process were scant. The first part described rules for valid corner exit points given a starting point on grids with edges that are either odd or even. This is what I want to explore in this blog post.

A rectangular area can have an either odd or even width and an either odd or even height. There are 4 possible combinations of these edge dimensions and for each one, a set of directions were supplied. I didn't just want to take these at face value. Some sort of proof was needed.

Space Filling Curve Algorithm by Lutz Tautenhahn

The proof below doesn't say anything about the existence of a path/curve that starts at one corner and ends at another covering all squares. All it does is show if a path is possible or not. Although it's relatively easy to prove that a path can be found for all of the scenarios below, we'll leave that for another time.

So what are the results?  In grids with two even edges, a path can be found from a corner to either adjacent corner, but not to the diagonal corner.  In grids with an even and odd side a path can be found to the diagonal corner and the adjacent corner on the even side but not the odd side.  Finally in a grid with two odd sides a path can be found to any other corner except if one of the odd edges is equal to 1.  In that case the path is only diagonal or along the odd side that isn't one.  That last bit isn't in the proof but it's obvious.

Allow directions in a space filling grid

Proving the directions is straight forward. If you imagine that the corner start square of a grid is covered, and every time you move, another square is covered, it becomes obvious that there are (wh-1) moves until the entire space is covered. The total number of moves is also equal to the number of up, down, left, and right moves. If you want to move right side of the grid the number of right moves minus the number of left moves has to be equal to (w-1). If you want to stay on the left, the number of left moves is equal to the number of right moves.  Similarly if you want to move to the bottom of the grid, the number of down moves minus the number of up moves is equal to (h-1). If you want to stay at the top, the number of down and up moves mus t be equal. When these equations are combined the directions above can be found.

Grid movement proof

It was redundant to prove both the Odd-Even and Even-Odd case as they are a diagonal mirroring of the same scenario.  It's also sufficient to only show the paths starting at the top left corner as paths starting in other corners are vertical and or horizontal mirrors of paths in the top left corner.

Although I didn't end up using the curve filling path in the first image, it helped me understand some of the subtleties of curve generation.

Wednesday, February 1, 2017

Voronoi Stippling

I'm trying to write my own software to convert normal raster images to a corresponding stippled version.  If you're unfamiliar with the term, stippling is a style of artwork that that uses are large number of small dots to represent a scene.  Areas with a high density of dots appear darker, while areas without dots are white.  Normally I like to have reasonably good understanding of a subject before I do a blog post about it, but in this case I'll give a quick outline of what I'm trying to achieve and then fill in the details in a later post.

So, how do you get an image like this,
Stippled Image

From something like this.
Original Image

To start with I'm trying to replicate the process in a paper by Adrian Secord.

Weighted Voronoi Stippling
Adrian Secord
2nd International Symposium on Non-Photorealistic Animation and Rendering (NPAR 2002)

The first step is to generate a number of points equal to the desired number of stippling points in the final image.  They can be randomly distributed over the image, or ideally via a rough approximation, close to the locations of the stippling dots in the output image.  As the process I'm about to describe is iterative, the closer the initial approximation is the final output, the faster the algorithm will converge.

After the points have been created, their Voronoi diagram is generated.  A Voronoi diagram produces a region for each input point that identifies all the locations that are closest to that dot than any other.  There are some tricks to making sure that the regions around the boundary are properly created but in general it's a simple call to scipy.spatial.Voronoi

That's great, but the process so far hasn't taken the input image into account(unless you use it to generate the initial point locations).  Now that a number of points and their corresponding regions are created, the weighted centroid of the area is calculated.  The centroid then becomes the location of the new point.  This process of creating Voronoi diagrams and calculating centroids is repeated until some convergence criteria is met and is better known as Lloyd's Algorithm.

In the image below, the black dots are the initial points and the white dots are the centroids of the Voronoi regions after an iteration.
Voronoi Diagram
Voronoi Diagram with centroids
Overlaying this diagram on the original image helps to illustrate the process a little better.  Obviously the image below is sparsely stippled and the object isn't recognisable with only 64 points.

Voronoi Diagram
Voronoi Iteration

I've generated the next diagram with I think 256 points to explain a subtle point.  The Voronoi region for each point is a convex polygon.  Therefore its centroid will always be located within itself.  As the number of points increases the size of the regions decreases.  This means that the process will take a long time to converge as the points can only move a small amount during each iteration.   That's why you need to place the initial points as close to their final location as possible.  Or do you?

I only discovered this problem after programming my solution and needed a quick way to overcome it.  If a large number of points take a long time to converge then a small number should converge quickly.  So I started with only two points.  After a few iterations, each point splits, like a cell undergoing Mitosis.  Each point splits into two new points diametrically opposed on a unit circle centred on the old point.  The direction of the split is randomly chosen.  (I have a hunch that continually splitting in the same direction will cause artifacts that will delay convergence)  These new points are then iterated on for a while and then split again.  The key is to get the points close to where they need to be early and then converge and add the detail later.
Voronoi Diagram
Dense Voronoi diagram
I have some ideas I want to try to make the image better.  There are obviously some linearity problems with the output.  White parts aren't as light as they are should be.  My code is a mess too.  It needs a clean up.  I'm having fun though.

Tuesday, January 17, 2017

Arrange Rectangles with Python to Design a Carton

Recently I've been interested in cardboard boxes and how to design and model them.  In my last blog post I described how to create a box in Fusion 360.  It works but it's a cumbersome, frustrating, and time consuming process.  What was needed was a way to go from concept to template quickly.  After looking around and not finding a tool for the job, I designed one.  It's called rectangleBuilder.

Laying out a template isn't too hard once you get the hang of it.  Draw all the rectangles that make up the box you intend to build and then dimension them.  Then define relationships between all the rectangles.  An example of a relationship can be seen in the plan below.  The bottom middle of rectangle B is in the same location as the top middle of rectangle A.

Box Template Plan

This can be done by defining a Rectangle class in python,  It is initialized with its width and height and by default is placed with its bottom left corner at (0,0).  The Rectangle object is also aware of the position of all its corners and edge midpoints.  When called with a code for one of these points it returns the xy coordinates in list form.  For example, calling an object with the argument 'TM' will return the coordinates of the top middle.  All of this shorthand is explained in the code.

Python Code
Define rectangle sizes

This is the part that I think is really cool.  Positioning the rectangles by defining geometric relationships.  In the code below the bottom middle of rectangle A is positioned at the origin.  Rectangle B is then positioned so that it's bottom middle is in the same position as the top middle of rectangle A.  The last line demonstrates the offset feature.  It essentially says, place rectangle M so that its bottom left corner is in the same position as the bottom right of rectangle B then shift it in the y direction by an amount (t+b).  The few lines of code above and below completely describe the template for the box.

Python Code
Position rectangles

This data can then be used to generate an SVG template in python to see if everything is as expected.  The outline to cut is in black and the bend lines are marked in red.

Box Template
Box template

That's all well and good but how do we transfer that to cardboard?  I don't have a plotter or laser cutter so it has to be done by hand.  After some experimenting, the below instructions are the easiest way to describe the rectangle edges so that they can be easily drawn.  The y coordinate and the start and end x coordinates of the horizontal lines are listed in increasing y order.

Stats and Horizontal lines

The same process is repeated for the vertical lines as well.

Vertical lines

Take the piece of cardboard that you plan to use and draw a base line in the direction of the x axis.  At each end of this line draw perpendicular guides in the direction of the y axis.  Along each of these lines mark all the y coordinates from the horizontal line list using the base line as the zero point.  You can then draw the line segments in the following manner.  77.0 , (84.0 - 124.5) means place a ruler on the 77.0 marks you would have drawn in the last step with the 0 mark on the left guide.  Then draw a line from 84.0 to 124.5 on the ruler.  By working through the list the design will gradually appear.  The vertical list is supplied just to be complete and isn't really needed as the ends of the horizontal lines can be joined to form the rectangles.

Edge guides

The layout looks like the SVG file so everything worked out OK.

Box Plan
Box to cut out

The bend lines can be perforated with a tool like the one below that's used for transferring dressmaking patterns.

Perforated Cardboard
Perforated Bend Lines

Cut out the outlines with a scalpel.

Flattened box

When folded together everything works quite well.  You can see in the image below what my intention was.  The top extends over to cover the entire top and although not necessary is a fun little design challenge.

Cardboard Box
Top overlap

As it's a tight fit, the box needs to be taped together to stop it springing open.

Cardboard Box
Assembled Box

Everything seems to fit nicely.  A little tweaking could make things better, but I'm happy with it.  In the end the fit all comes down to how accurately the template is drawn and bent.  Folding across the corrugations is easy to get right, but folding in the same direction as the corrugations is hard and I need to come up with a better way to do it.  When folding with the corrugation it's not uncommon for the fold to follow the middle of the corrugation.  This could be a couple mm away from from the line.
Get The Code!