{"id":505,"date":"2022-09-02T04:03:36","date_gmt":"2022-09-02T04:03:36","guid":{"rendered":"https:\/\/timallanwheeler.com\/blog\/?p=505"},"modified":"2022-09-02T04:04:59","modified_gmt":"2022-09-02T04:04:59","slug":"embedded-contours-with-marching-squares","status":"publish","type":"post","link":"https:\/\/timallanwheeler.com\/blog\/2022\/09\/02\/embedded-contours-with-marching-squares\/","title":{"rendered":"Embedded Contours with Marching Squares"},"content":{"rendered":"\n

I recently started working on new chapters for Algorithms for Optimization<\/a>, starting with a new chapter on quadratic programs. We already have a chapter on linear programs, but quadratic programs are sufficiently different and are ubiquitous enough that giving them their own chapter felt warranted.<\/p>\n\n\n

\n
\"\"
A simple example quadratic program.<\/figcaption><\/figure>\n<\/div>\n\n\n

The chapter starts with a general formulation for quadratic programs and transforms the problem a few times into simpler and simpler formulations, until we get a nonnegative least squares problem<\/a>. This can be solved, and then we can use that solution to back our way up the transform stack back to a solution of our original solution.<\/p>\n\n\n\n

In texts, these problem transforms are typically solely the realm of mathematical manipulations. While math is of course necessary, I also tried hard to depict what these various problem representations looked like to help build an intuitive understanding of what some of the transformations were. In particular, there was one figure where I wanted to show how the equality constraints of a general-form quadratic program could be removed when formulating a lower-dimensional least squares program:<\/p>\n\n\n\n

\"\"
A 3d general-form quadratic program with an equality constraint on the left, and a 2d least squares program on the right.<\/figcaption><\/figure>\n\n\n\n

This post is about this figure, and more specifically, about the blue contour lines on the left side.<\/p>\n\n\n\n

What is this Figure?<\/h2>\n\n\n\n

The figure on left shows the following quadratic program:<\/p>\n\n\n

\n
\"\"
general-form quadratic program<\/figcaption><\/figure>\n<\/div>\n\n\n

As we can see, the quadratic program is constrained to lie within a square planar region \\(0 \\leq x \\leq 1\\), \\(0 \\leq x_3 \\leq 1\\), and \\(x_2 = x_3\\). The feasible set is shown in blue.<\/p>\n\n\n\n

The right figure on the right shows the least-squares program that we get when we factor out the equality constraint. No need to get into the details here, it isn’t really important for this post, but we can use our one equality equation to remove one variable in our original problem, producing a 2d problem without any equality constraints:<\/p>\n\n\n

\n
\"\"
least squares program<\/figcaption><\/figure>\n<\/div>\n\n\n

The original feasible set is rotated and stretched, and the objective function contours are rotated and stretched, but other than that we have a similar topology and a pretty standard problem that we can solve.<\/p>\n\n\n\n

How do we Make the Figures?<\/h2>\n\n\n\n

I make the plots with PGFPlots.jl<\/a>, which uses the pgfplots<\/a> LaTeX package under the hood. It is a nice synergy of Julia coding and LaTeX within the LaTeX textbook.<\/p>\n\n\n\n

Right Plot<\/h3>\n\n\n\n

The right plot is pretty straight forward. It has a contour plot for the objective function:<\/p>\n\n\n\n

Plots.Contour(x -> norm(A_lsi*x - b_lsi, 2),\n              xdomain, ydomain, xbins=600, ybins=600)<\/pre>\n\n\n\n
some over + under line plots and a fill-between for the feasible region:<\/p>\n\n\n\n
x1_domain_feasible = (0.0, (0.5+\u221a2)\/2)\npush!(plots,\n      Plots.Linear(x1->max(-x1, x1 - sqrt(2)), x1_domain_feasible,\n                   style=\"name path=A, draw=none, mark=none\"))\npush!(plots,\n      Plots.Linear(x1->min(x1, 1\/2-x1), x1_domain_feasible,\n                   style=\"name path=B, draw=none, mark=none\"))\npush!(plots,\n      Plots.Command(\"\\\\addplot[pastelBlue!50] fill between[of=A and B]\"))<\/pre>\n\n\n\n
and a single-point scatter plot plus a label for the solution:<\/p>\n\n\n\n
push!(plots,\n      Plots.Scatter([x_val[1]], [x_val[2]],\n      style=\"mark=*, mark size=1.25, mark options={solid, draw=black, fill=black}\"))\npush!(plots,\n      Plots.Command(\"\\\\node at (axis cs:$(x_val[1]),$(x_val[2])) [anchor=west] {\\\\small \\\\textcolor{black}{\\$y_2^\\\\star\\$}}\"))<\/pre>\n\n\n\n
These plots get added to an Axis<\/code> and we’ve got our plot.<\/p>\n\n\n\n
Left Plot<\/h3>\n\n\n\n
The left plot is in 3D. If we do away with the blue contour lines, its actually even simpler than the right plot. We simply have a square patch, a single-point scatter plot, and our label:<\/p>\n\n\n\n
Axis([\n   Plots.Command(\"\\\\addplot3[patch,patch type=rectangle, faceted color=none, color=pastelBlue!40] coordinates {(0,0,0) (1,0,0) (1,1,1) (0,1,1)}\"),\n   Plots.Linear3([x_val[1]], [x_val[2]], [x_val[3]], mark=\"*\", style=\"mark size=1.1, only marks, mark options={draw=black, fill=black}\"),\n   Plots.Command(\"\\\\node at (axis cs:$(x_val[1]),$(x_val[2]),$(x_val[3])) [anchor=west] {\\\\small \\\\textcolor{black}{\\$x^\\\\star\\$}}\"),\n  ], \n  xlabel=L\"x_1\", ylabel=L\"x_2\", zlabel=L\"x_3\", width=\"1.2*$width\", \n  style=\"view={45}{20}, axis equal\",\n  legendPos=\"outer north east\"\n)<\/pre>\n\n\n\n
The question is – how do we create that contour plot, and how do we get it to lie on that blue patch?<\/p>\n\n\n
\n
\"\"<\/figure>\n<\/div>\n\n\n
How do you Draw Contours?<\/h2>\n\n\n\n
Contour plots don’t seem difficult to plot, until you have to sit down and plot them. They are just lines of constant elevation – in this cases places where \\(\\|Ax – b\\| = c\\) for various constants \\(c\\). But, unlike line plots, where you can just increase \\(x\\) left to right and plot successive \\((x, f(x))\\) pairs, for contour plots you don’t know where your line is going to be ahead of time. And that’s assuming you only have one line per contour value – you might have multiple:<\/p>\n\n\n
\n
\"\"
A contour plot with multiple contour lines per level.<\/figcaption><\/figure>\n<\/div>\n\n\n
So, uh, how do we do this?<\/p>\n\n\n\n
The first thing I tried was trying to write my function out the traditional way – as \\((x_1, f(x_1))\\) pairs. Unfortunately, this doesn’t work out, because the contours aren’t 1:1.<\/p>\n\n\n\n
The second thing I briefly considered was a search method. If I fix my starting point, say, \\(x_1 = 0, x_2 = 0.4\\), I can get the contour level \\(c = \\|Ax – b\\| \\). Then, I could search toward the right in an arc to find where the next contour line was, maybe bisect over an arc one step distance away and proceed that way:<\/p>\n\n\n
\n
\"\"
An iterative search procedure for following a contour line, proceeding left to right.<\/figcaption><\/figure>\n<\/div>\n\n\n
I think this approach would have worked for this figure, since there were never multiple contour lines for the same contour value, but it would have been more finicky, especially around the boundaries, and it would not have generalized to other figures.<\/p>\n\n\n\n
Marching Squares<\/h2>\n\n\n\n
I happened upon a talk by Casey Muratori on YouTube, “Papers I have Loved<\/a>“, a few weeks ago. One of the papers Casey talked about was Marching Cubes, A High Resolution 3D Surface Construction Algorithm. <\/em>This 1987 paper by W. Lorenson and H. Cline presented a very simply way to generate surfaces (in 3D or 2D) in places of constant value: \\(f(x) = c\\).<\/p>\n\n\n\n
Casey talked about it in the context of rendering metaballs<\/em>, but when I found myself presented with the contour problem, it was clear that marching cubes (or squares, rather, since its 2d) was relevant here too.<\/p>\n\n\n\n
The algorithm is a simple idea elegantly applied.<\/p>\n\n\n\n
The contour plot we want to make is on a 2d surface in a 3d plot, but let’s just consider the 2d surface for now. Below we show that feasible square along with our objective function, \\(f(x) = \\|Ax – b\\|\\):<\/p>\n\n\n
\n
\"\"
An image plot showing our object function. We want to find contour lines.<\/figcaption><\/figure>\n<\/div>\n\n\n
Our objective is to produce contour lines, the same sort as are produced when we make a PGFPlots contour plot:<\/p>\n\n\n
\n
\"\"<\/figure>\n<\/div>\n\n\n
We already know that there might be multiple connected contour lines for a given level, so a local search method won’t be sufficient. Instead, marching squares uses a global approach combined with local refinement.<\/p>\n\n\n\n
We begin by evaluating a grid of points over our plot:<\/p>\n\n\n
\n
\"\"<\/figure>\n<\/div>\n\n\n
We’ll start with a small grid, but for higher-quality images we’d make it finer.<\/p>\n\n\n\n
Let’s focus on the contour line \\(f(x) = 3\\). We now determine, for each vertex in our mesh, whether its value is above or below our contour target of 3. (I’ll plot the contour line too for good measure):<\/p>\n\n\n
\n
\"\"
Blue points are where the objective function <= 3, and yellow points are where it is > 3.<\/figcaption><\/figure>\n<\/div>\n\n\n
Marching squares will now tile the image based on each 4-vertex square, and the binary values we just evaluated. Any squares with all vertices equal can simply be ignored – the contour line doesn’t pass through them. This leaves the following tiles:<\/p>\n\n\n
\n
\"\"
The white tiles are those that contain our contour line(s).<\/figcaption><\/figure>\n<\/div>\n\n\n
Now, each remaining tile comes in multiple possible “flavors” – the various ways that its vertices can be above or below the objective value. There are \\(2^4\\) possible tiles, two of which we’ve already discarded, leaving 14 options. The Wikipedia article<\/a> has a great  set of images showing all of the tile combinations, but here we will concern ourselves with two more kinds: tiles with one vertex below or one vertex above, and tiles with two neighboring vertices above and two neighboring vertices below:<\/p>\n\n\n
\n
\"\"
The white squares have one vertex above or one vertex below. The gray squares have two neighboring vertices above and two neighboring vertices below. <\/figcaption><\/figure>\n<\/div>\n\n\n
Let’s start with the first type. These squares have two edges where the function value changes. We know, since we’re assuming our objective function is continuous, that our contour line must pass through those edges. As such, if we run a bisection search along each of those edges, we can quickly identify those crossing points:<\/p>\n\n\n
\n
\"\"
The red points are where bisection search has found a crossing point.<\/figcaption><\/figure>\n<\/div>\n\n\n
We then get our contour plot, for this particular contour value, by drawing all of the line segments between these crossing points:<\/p>\n\n\n
\n
\"\"<\/figure>\n<\/div>\n\n\n
We can get higher resolution contour lines by using a finer grid:<\/p>\n\n\n
\n
\"\"
The same contour line with double the resolution.<\/figcaption><\/figure>\n<\/div>\n\n\n
To get our full contour plot, we simply run marching squares on multiple contour levels. To get the original plot, I run this algorithm to get all of the little tile line segments, then I render them as separate Linear3 line segments on top of the 3d patch. Pretty cool!<\/p>\n\n\n\n
\"\"
The created figure<\/figcaption><\/figure>\n\n\n\n
What does the Code Look Like?<\/h2>\n\n\n\n
The code for marching squares is pretty elegant.<\/p>\n\n\n\n
We start off with defining our grid:<\/p>\n\n\n\n
nbins = 41\nx1s = collect(range(0.0, stop=1.0, length=nbins))\nx2s = collect(range(0.0, stop=1.0, length=nbins))\nis_above = falses(length(x1s), length(x2s))<\/pre>\n\n\n\n
Next we define a bisection method for quickly interpolating between two vertices:<\/p>\n\n\n\n
function bisect(\n    u_target::Float64, x1_lo::Float64, x1_hi::Float64, \n    x2_lo::Float64, x2_hi::Float64, k_max::Int=8)\n    f_lo = norm(A*[x1_lo, x2_lo, x2_lo]-b)\n    f_hi = norm(A*[x1_hi, x2_hi, x2_hi]-b)\n    \n    if f_lo > f_hi # need to swap low and high\n        x1_lo, x1_hi, x2_lo, x2_hi = x1_hi, x1_lo, x2_hi, x2_lo\n    end\n    \n    x1_mid = (x1_hi + x1_lo)\/2\n    x2_mid = (x2_hi + x2_lo)\/2\n    for k in 1 : k_max\n        f_mid = norm(A*[x1_mid, x2_mid, x2_mid] - b)\n        if f_mid > u_target\n            x1_hi, x2_hi, f_hi = x1_mid, x2_mid, f_mid\n        else\n            x1_lo, x2_lo, f_lo = x1_mid, x2_mid, f_mid\n        end\n        x1_mid = (x1_hi + x1_lo)\/2\n        x2_mid = (x2_hi + x2_lo)\/2\n    end\n    \n    return (x1_mid, x2_mid)\nend<\/pre>\n\n\n\n
In this case I’ve baked in the objective function, but you could of course generalize this and pass one in.<\/p>\n\n\n\n
We then iterate over our contour levels. For each contour level (u_target), we evaluate the grid and  then evaluate all of our square patches and add little segment plots as necessary:<\/p>\n\n\n\n
for (i,x1) in enumerate(x1s)\n    for (j,x2) in enumerate(x2s)\n        is_above[i,j] = norm(A*[x1, x2, x2] - b) > u_target\n    end\nend\nfor i in 2:length(x1s)\n    for j in 2:length(x2s)\n        # bottom-right, bottom-left, top-right, top-left\n        bitmask = 0x00\n        bitmask |= (is_above[i,j-1] << 3)\n        bitmask |= (is_above[i-1,j-1] << 2)\n        bitmask |= (is_above[i,j] << 1)\n        bitmask |= (is_above[i-1,j])\n        plot = true\n        x1_lo1 = x1_hi1 = x2_lo1 = x2_hi1 = 0.0\n        x1_lo2 = x1_hi2 = x2_lo2 = x2_hi2 = 0.0\n        if bitmask == 0b1100 || bitmask == 0b0011\n            # horizontal cut\n            x1_lo1, x1_hi1, x2_lo1, x2_hi1 = x1s[i-1], x1s[i-1], x2s[j-1], x2s[j]\n            x1_lo2, x1_hi2, x2_lo2, x2_hi2 = x1s[i], x1s[i], x2s[j-1], x2s[j]\n        elseif bitmask == 0b0101 || bitmask == 0b1010\n            # vertical cut\n            x1_lo1, x1_hi1, x2_lo1, x2_hi1 = x1s[i-1], x1s[i], x2s[j-1], x2s[j-1]\n            x1_lo2, x1_hi2, x2_lo2, x2_hi2 = x1s[i-1], x1s[i], x2s[j], x2s[j]\n        elseif bitmask == 0b1000 || bitmask == 0b0111\n            # cut across bottom-right\n            x1_lo1, x1_hi1, x2_lo1, x2_hi1 = x1s[i-1], x1s[i], x2s[j-1], x2s[j-1]\n            x1_lo2, x1_hi2, x2_lo2, x2_hi2 = x1s[i], x1s[i], x2s[j-1], x2s[j]\n        elseif bitmask == 0b0100 || bitmask == 0b1011\n            # cut across bottom-left\n            x1_lo1, x1_hi1, x2_lo1, x2_hi1 = x1s[i-1], x1s[i], x2s[j-1], x2s[j-1]\n            x1_lo2, x1_hi2, x2_lo2, x2_hi2 = x1s[i-1], x1s[i-1], x2s[j-1], x2s[j]\n        elseif bitmask == 0b0010 || bitmask == 0b1101\n            # cut across top-right\n            x1_lo1, x1_hi1, x2_lo1, x2_hi1 = x1s[i-1], x1s[i], x2s[j], x2s[j]\n            x1_lo2, x1_hi2, x2_lo2, x2_hi2 = x1s[i], x1s[i], x2s[j-1], x2s[j]\n        elseif bitmask == 0b0001 || bitmask == 0b1110\n            # cut across top-left\n            x1_lo1, x1_hi1, x2_lo1, x2_hi1 = x1s[i-1], x1s[i], x2s[j], x2s[j]\n            x1_lo2, x1_hi2, x2_lo2, x2_hi2 = x1s[i-1], x1s[i-1], x2s[j-1], x2s[j]\n        else\n            plot = false\n        end\n        if plot\n            x1a, x2a = bisect(u_target, x1_lo1, x1_hi1, x2_lo1, x2_hi1)\n            x1b, x2b = bisect(u_target, x1_lo2, x1_hi2, x2_lo2, x2_hi2)\n            push!(contour_segments, \n                Plots.Linear3([x1a, x1b], [x2a, x2b], [x2a, x2b],\n                               style=\"solid, pastelBlue, mark=none\")\n            )\n        end\n    end\nend<\/pre>\n\n\n\n
In this code, all we ever do is figure out which two line segments to bisect over, and if we get then, run the bisection and plot the ensuing line. Its actually quite neat and compact.<\/p>\n\n\n\n
Conclusion<\/h2>\n\n\n\n
So there you have it - a home-rolled implementation for contour plotting turns out to be useful for embedding 2d contours in a 3d plot.<\/p>\n\n\n\n
Note that the implementation here wasn't truly general<\/mark><\/span><\/span>Which is fine - it was sufficient to generate the figure. As long as we are aware of where we cut corners, we are okay.<\/span><\/span><\/span><\/span>. In this case, we didn't handle some tile types, namely those where opposing corners have the same aboveness, but neighboring corners have differing aboveness. In those cases, we need to bisect across all four segments, and we typically need to remove ambiguity by sampling in the tile's center to determine which way the segments cross.<\/p>\n\n\n\n
I find it really satisfying to get to apply little bits of knowledge like this from time to time. Its part of computer programming that is endlessly rewarding. I have often been amazed by how a technique that I learned about in one context surprised me by being relevant in something else later on.<\/p>\n","protected":false},"excerpt":{"rendered":"
I recently started working on new chapters for Algorithms for Optimization, starting with a new chapter on quadratic programs. We already have a chapter on linear programs, but quadratic programs are sufficiently different and are ubiquitous enough that giving them their own chapter felt warranted. The chapter starts with a general formulation for quadratic programs […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-505","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/posts\/505","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/comments?post=505"}],"version-history":[{"count":31,"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/posts\/505\/revisions"}],"predecessor-version":[{"id":560,"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/posts\/505\/revisions\/560"}],"wp:attachment":[{"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/media?parent=505"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/categories?post=505"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/timallanwheeler.com\/blog\/wp-json\/wp\/v2\/tags?post=505"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}