Image processing doesn’t always have to do with scientific computing. Image processing techniques can easily be applied to create artistic filters capable of producing art that would be impossible or difficult for a human to recreate by hand. Previously, I wrote about an algorithm that can make an image look like an oil painting. In this article, I’ll be discussing an algorithm that can be used to roughly simulate what an image would look like if it were in stained glass.
Properties of stained glass
Stained glass has been around for centuries. This type of art form has several properties that we will try to simulate. Stained glass has regions of solid colors, and these regions are separated by black lines. For this algorithm, we will be using the most basic polygon, triangles. The general flow is to generate a set of points then create triangles from said points. For each triangle, we can calculate the average color. Then we can simply draw solid color triangles, and draw black lines for all the edges of the triangles and we’re done!
Delaunay triangulation
Creating the triangulation is the tricky part. How do we take a set of points and determine a good set of triangles from these points? For this project, I used Delaunay triangulation. Basically, the goal is to create triangles such that no point lies within the circumcircle of any triangle. This is a conceptually simple task, but can actually be quite annoying to solve quickly. Since I was more interested in the end result, I decided to use some code geompack. There is no need to reinvent the wheel here. Originally, I used the dtris2 function to calculate the Delaunay triangulation quickly. However, it never seemed to work quite right, and often gave incorrect triangulations. So I opted to use the painfully slow points_delaunay_naive_2d function instead. Most of my articles revolve around high performance computing, but in this case, I’ll settle for a naive O(N^4) algorithm. Ouch! I did, however, modify it slightly so that it uses OpenMP and completes in one pass rather than two. With an 8 core computer, it’s around 16 times faster which eliminates a little bit of the pain.
Filling the triangles
In the interests of finishing quickly, I also took a naive approach to determining the color value of each triangle. To do this, I initialize an array ppColor[n_triangles][3] to store the red, green and blue components for each triangle to zero. Then, for each pixel, for each triangle, test to see if that pixel is in the triangle. If it is, take the original RGB value of that pixel, and add it to the ppColor value for that triangle index. Below is the simple, unoptimized code. Again, this code is extremely unoptimized, but it is quite clear what is going on.
unsigned long **ppColor = new unsigned long*[nTriangles]; // could've been a flattened array instead, but no matter. int *pPixelCount = new int[nTriangles]; memset(pPixelCount, 0, sizeof(int)*nTriangles); for (int i=0; i < nTriangles; i++) { ppColor[i] = new unsigned long[3]; // r,g,b ppColor[i][0] = ppColor[i][1] = ppColor[i][2] = 0; } // for each pixel, for each triangle, see which triangle it's in. #pragma omp parallel for shared(ppColor, pPixelCount) for (int y=0; y < height; y++) { for (int x=0; x < width; x++) { point p; p.x = (float)x; p.y = (float)y; for (int i=0; i < nTriangles; i++) { int v1 = pTriVert[i*3+0]; int v2 = pTriVert[i*3+1]; int v3 = pTriVert[i*3+2]; point p0, p1, p2; p0.x = pDPoints[v1*2]; p0.y = pDPoints[v1*2+1]; p1.x = pDPoints[v2*2]; p1.y = pDPoints[v2*2+1]; p2.x = pDPoints[v3*2]; p2.y = pDPoints[v3*2+1]; int matched = IsPointInTriangle(p0, p1, p2, p); if (matched) { #pragma omp critical { Color color; // GetPixel is slow and using LockBits would be much faster, but I'm really going for simplicity here... pBitmap->GetPixel(x,y,&color); ppColor[i][0] += color.GetR(); ppColor[i][1] += color.GetG(); ppColor[i][2] += color.GetB(); pPixelCount[i]++; } i=nTriangles; // to exit; } } } } #pragma omp parallel for shared(ppColor) for (int i=0; i < nTriangles; i++) { if (pPixelCount[i]) { ppColor[i][0] /= pPixelCount[i]; ppColor[i][1] /= pPixelCount[i]; ppColor[i][2] /= pPixelCount[i]; } else { // extremely narrow and spread out triangles are possible where no pixels are fully in bounds. // just set it to black since it'll look like a line anyway ppColor[i][0] = ppColor[i][1] = ppColor[i][2] = 0; } }
Once the totals are added up, the color can be divided by total number of pixels in each triangle, this giving the average color for that triangle. Once that is done, drawing the triangles and lines on the borders of the triangles is trivial.
Picking the points
Picking the points to use occurs in two phases. One approach would be to generate a bunch of random points. However, it would become important to check for duplicate points. Instead of adding this burden, I chose another path. For each pixel, there is a small probability of becoming a point. Doing this approach prevents duplicate points from being formed. For artistic reasons, I maintain different probabilities for image borders. Also, I purposefully ensured that each corner of the image is a point. Below is a first pass at the algorithm.
vector <sPoint> originalPoints; for (int i=0; i < height; i++) { for (int j=0; j < width; j++) { int curRand = rand()%10000; float fRand = (float)curRand / 10000.0f; float density = imageDensity; if (j == 0 || i==0 || j == width-1 || i == height-1) { density = borderDensity; } if((j==0 && i==0) || (j==0 && i==height-1) || (j==width-1 && i==0) || (j==width-1 && i==height-1)) { density = 1000; // put vertices in corner } if (fRand <= density) { sPoint curPoint; curPoint.x = j; curPoint.y = i; originalPoints.push_back(curPoint); } } }
The image looks nice enough, but it doesn’t look good if two points are too close together. Therefore, I created a function to remove all the duplicate points. The image does look much better. Do note, however, that it isn’t perfect. It can cause some interesting artifacting at the bottom of the image. If this is undesirable, you can always modify it to protect the borders only against other border pixels, then do all other pixels normally. The image looks better without the crowded pixels!
vector <sPoint> CImageProcessor::RemoveClosePoints(vector <sPoint> originalPoints, float minDistance, int width, int height) { vector <sPoint> vettedPoints; int nPoints = originalPoints.size(); sPoint cornerP; cornerP.x=0; cornerP.y=0; vettedPoints.push_back(cornerP); cornerP.x=width-1; cornerP.y=0; vettedPoints.push_back(cornerP); cornerP.x=width-1; cornerP.y=height-1; vettedPoints.push_back(cornerP); cornerP.x=0; cornerP.y=height-1; vettedPoints.push_back(cornerP); for (int i=0; i < nPoints; i++) { sPoint p = originalPoints.at(i); if ((p.x == 0 && p.y == 0) || (p.x==width-1 && p.y==0) | (p.x == 0 && p.y == height-1) || (p.x==width-1 && p.y== height-1)) { continue; } else { int tooClose = 0; for (int j=0; j < vettedPoints.size(); j++) { sPoint vetP = vettedPoints.at(j); double dX = vetP.x - p.x; double dY = vetP.y - p.y; double dist = sqrt(dX*dX + dY*dY); if (dist < minDistance) { tooClose=1; j=vettedPoints.size(); } } if (tooClose == 0) vettedPoints.push_back(p); } } return vettedPoints; }
You can download the full source code here. It has everything except for the Delaunay triangulation code.