Re-posted from: https://cancandan.github.io/julia/graphics/cuda/2022/05/07/triangles.html
Suppose we want to draw a batch of images, where each image is made up of randomly positioned and colored triangles, blending, so it will look like this:
and then find the euclidean distance of each such image to a given target image. Why on earth am I doing this? Well this turns into an interesting optimization problem of finding the closest triangle image and also an excuse to practise Julia. The inspiration is from this repo.
Towards framing this as an optimization problem, we will represent a triangle as a vector of size 10, made up of floating point numbers between 0 and 1. Four numbers of this vector are for the color of the triangle; r,g,b and alpha, and for the three vertices of the triangle we need 6 numbers, (x1,y1), (x2,y2), (x3,y3). Hence, if we want to draw M images, each image having N triangles, we need a matrix of size (10,N,M), which will be our parameters matrix. I want to randomly create such a matrix and min-max scale it along the triangle dimension, by which I mean, for each image, I first find the minimum and maximum of a triangle parameter among the N triangles, and then subtract from the parameter this minimum and then divide the result by the difference between the maximum and the minimum. I want to end up with an array of size (3,w,h,M) for the images, where w is width and h is height, and an array of size M for the distances. Let’s see how fast we can do this.
First order of business is setting this up, note that I am scaling the numbers to the given width and height:
The first thing that comes to mind is to use a 2d graphics library for drawing, and since the Cairo lib is available, let’s try that, the function below is drawing the triangles on a blank white Cairo canvas, and copying it to the img array at the end:
Now let’s draw each image in this fashion:
Benchmarking this with @btime withCairo();
I see 428.101 ms (3157 allocations: 205.09 MiB)
.
The cross product method
Now something really cool. Move your mouse inside and outside of the triangle below, you will see a bar chart depicting the magnitude and direction of the 3rd components of the cross products of vectors AB with AP (reds), BC with BP (greens) and CA with CP (blues). Observe that all those bars point same direction only inside the triangle!
let sketch = function(p) {
p.setup = function(){ p.createCanvas(400, 400); p.background(1); };
// p.draw = function() { // p.background(220); // p.triangle(30, 75, 58, 20, 86, 75); // }; // };
p.drawArrow = function(base, vec, myColor) { p.push(); p.stroke(myColor); p.strokeWeight(3); p.fill(myColor); p.translate(base.x, base.y); p.line(0, 0, vec.x, vec.y); p.rotate(vec.heading()); let arrowSize = 7; p.translate(vec.mag() - arrowSize, 0); p.triangle(0, arrowSize / 2, 0, -arrowSize / 2, arrowSize, 0); p.pop(); };
p.draw = function() { p.background(220); [x1,y1,x2,y2,x3,y3]=[30, 375, 58, 20, 386, 75]; // triangle(x1,y1,x2,y2,x3,y3); r = p.color(255, 0, 0) g = p.color(0, 255, 0) b = p.color(0, 0, 255)
p.strokeWeight(4); p.stroke(r); // p.line(x1, y1, x2, y2); p.drawArrow(p.createVector(x1,y1), p.createVector(x2-x1,y2-y1), r);
p.stroke(g); // p.line(x2, y2, x3, y3); p.drawArrow(p.createVector(x2,y2), p.createVector(x3-x2,y3-y2), g);
p.stroke(b); // p.line(x3, y3, x1, y1); p.drawArrow(p.createVector(x3,y3), p.createVector(x1-x3,y1-y3), b);
p.stroke(p.color(0, 0, 0));
p.strokeWeight(1);
p.text("A", x1-10, y1+10); p.text("B", x2-20, y2); p.text("C", x3, y3-10); p.text("P", p.mouseX-20, p.mouseY);
if (p.mouseX <= 400 && p.mouseX >= 0 && p.mouseY <= 400 && p.mouseY >= 0) {
// p.line(x1, y1, p.mouseX, p.mouseY); // p.line(x2, y2, p.mouseX, p.mouseY); // p.line(x3, y3, p.mouseX, p.mouseY); p.drawArrow(p.createVector(x1,y1), p.createVector(p.mouseX-x1,p.mouseY-y1), r); p.drawArrow(p.createVector(x2,y2), p.createVector(p.mouseX-x2,p.mouseY-y2), g); p.drawArrow(p.createVector(x3,y3), p.createVector(p.mouseX-x3,p.mouseY-y3), b);
// p.push(); cr1=(x2-x1)*(p.mouseY-y1)-(y2-y1)*(p.mouseX-x1) cr2=(x3-x2)*(p.mouseY-y2)-(y3-y2)*(p.mouseX-x2) cr3=(x1-x3)*(p.mouseY-y3)-(y1-y3)*(p.mouseX-x3) p.fill(r); p.rect(300, 300, 20, -cr1/1000);
p.translate(30,0); p.fill(g); p.rect(300, 300, 20, -cr2/1000);
p.translate(30,0); p.fill(b); p.rect(300, 300, 20, -cr3/1000);
p.translate(-60,0); p.text("3rd component of \n Cross Products =>", 180, 280); // p.pop(); }
}; };
new p5(sketch, 'container');
Whats great about this is that cross products are just multiplications and subtractions, perfect job to parallelize with a GPU.
What needs to be done is clear, for each of the M images, for each of the N triangles, our operation is to update a pixel color to blend with the current triangles color if that pixel is inside the triangle. We will parallelize this operation with a CUDA kernel, shown below:
We will need to pass our parameters and target array to the GPU, and then call the kernel with @cuda
. We can create white canvases with CUDA.ones
here, so no need to pass it.
Benchmarking this I see 120.315 ms (38689 allocations: 52.53 MiB)
Awesome, that’s about 4x speedup. Note that I benchmarked with --check-bounds=no
, which is a startup option that you pass when launching julia that disables the bounds checking. In the next post, I will talk about the very cool PGPE algorithm used in evojax to steer these images towards a target image. You can see one example of this here.
I am pretty new to Julia, please let me know if you have any comments, suggestions.