Generating sloping 3D models from Chinese characters
Tags: python, openscad, 3dprinting, opensource
Published: 2022-03-14 19:01:22

In this weekend project that eventually took a few months I combined three passions: Chinese, programming and 3D printing.

I was wondering whether it’d be possible to generate a 3D model of a Chinese character with the height representing the order of writing.

# Results

Some examples of what you can generate:

{% embed https://www.youtube.com/embed/tn4bemVyHxs?loop=1&playlist=tn4bemVyHxs %}

# How it works

## Reading the outlines of Chinese characters

I’ll be using the character 好 (hăo) for all the examples here. It means “good” and is the second part of 你好 (nĭ hăo) which means “hello” in Chinese. 好 looks as follows in a common font:

![好](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/s12x7n8kxvtzm3uv2xjl.png)

Each Chinese character is written in multiple strokes, and the order in which you draw the strokes is important to get the character visually right. In 好, there are 6 strokes.

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/d4pj97n8vxgka47jzl3t.gif">
  <figcaption>The 6 strokes of 好. Source: [Hanzi Writer](https://hanziwriter.org/demo.html)</figcaption>
</figure>

There’s a great database of Chinese strokes available over at [Make me a Hanzi](https://github.com/skishore/makemeahanzi). Every character is defined in a JSON, for example: 

```jsx
{
  "character": "好",
  "strokes": ["M 330 202 Q 361 175 399 134 Q 415 119 424 118 Q 433 118 439 128 Q 446 138 442 170 Q 435 206 361 247 L 319 270 Q 292 286 258 304 Q 237 314 240 335 Q 261 393 281 453 L 293 492 Q 317 568 337 644 Q 347 690 366 715 Q 379 737 373 750 Q 360 769 313 797 Q 294 810 276 801 Q 263 794 273 778 Q 303 733 247 486 L 236 442 Q 218 373 195 336 Q 185 314 206 296 Q 254 268 294 233 L 330 202 Z","M 294 233 Q 287 226 281 217 Q 250 180 196 143 Q 183 134 165 124 Q 149 114 133 104 Q 120 95 131 92 Q 212 86 327 199 Q 328 200 330 202 L 361 247 Q 406 322 421 385 Q 449 488 463 510 Q 473 526 458 537 Q 416 576 387 569 Q 374 565 378 550 Q 387 531 387 507 L 385 481 Q 384 469 382 455 Q 375 376 319 270 L 294 233 Z","M 387 507 Q 341 501 293 492 L 247 486 Q 183 479 115 468 Q 94 465 61 471 Q 48 471 45 462 Q 41 450 49 441 Q 68 422 96 400 Q 106 396 118 402 Q 190 436 236 442 L 281 453 Q 320 463 362 474 Q 372 478 385 481 C 414 489 417 511 387 507 Z","M 671 521 Q 788 635 822 648 Q 843 655 835 672 Q 831 688 760 725 Q 739 735 716 725 Q 661 703 575 676 Q 553 669 498 669 Q 473 669 482 648 Q 491 635 511 623 Q 544 605 578 627 Q 597 636 691 676 Q 706 682 719 673 Q 732 664 726 649 Q 693 595 655 531 C 640 505 649 500 671 521 Z","M 717 430 Q 702 497 671 521 L 655 531 Q 648 535 640 538 Q 618 547 608 540 Q 595 533 608 519 Q 645 491 653 444 Q 656 434 659 421 L 668 384 Q 701 204 658 103 Q 643 76 607 83 Q 576 89 548 94 Q 536 97 542 85 Q 546 78 564 65 Q 604 31 618 5 Q 628 -14 645 -11 Q 660 -10 687 17 Q 775 107 726 391 L 717 430 Z","M 726 391 Q 783 397 947 397 Q 966 398 971 406 Q 977 416 960 430 Q 909 467 848 454 Q 793 445 717 430 L 659 421 Q 562 409 452 393 Q 431 392 447 375 Q 460 362 478 357 Q 497 351 514 356 Q 586 375 668 384 L 726 391 Z"],
  "medians": [[[282,788],[307,769],[327,733],[264,465],[216,321],[235,298],[386,194],[411,166],[424,133]],[[390,556],[417,530],[424,516],[422,504],[387,361],[338,255],[304,207],[260,165],[206,127],[137,97]],[[59,457],[107,434],[373,491],[380,501]],[[493,656],[517,646],[550,644],[680,692],[706,699],[743,696],[771,669],[765,657],[677,546],[674,535],[663,536]],[[613,530],[637,519],[659,499],[674,474],[687,432],[711,289],[709,166],[692,92],[672,59],[648,41],[551,85]],[[449,384],[504,377],[860,427],[906,426],[960,412]]]
}
```

The outlines of strokes of this character (好) are defined in the `strokes` array. They are stored in an SVG path format.

The medians of each stroke are defined in the `medians` array.

There are 6 strokes, so each of these arrays contains 6 items.

All points are on a 1024x1024 grid.

Let’s read this data using Svgpathtools, a Python library:

```python
import parse_path from svgpathtools
stroke_path = parse_path(character_data["strokes"][0])
```

Svgpathtools gives us back a Path object. To draw these strokes with Matplotlib, we need to sample some points on that stroke:

```python
points = []
# for every segment in the path...
for segment in stroke_path:
    # take 4 samples
    for i in [0, 0.25, 0.5, 0.75]
        point = segment.point(i)
        ps.append(point)
points.append(points[0])
xs_stroke, ys_stroke = zip(*points)
plt(xs_stroke, ys_stroke, "g-")
```

This will plot the outlines of the strokes. We can also plot the medians:

```python
xs_medians, ys_medians = zip(*character_data["medians"][0])
plt(xs_medians, ys_medians, "ro", markersize=2)
```

And after some subplotting and aligning, we get this:

![Plot 好 1.png](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/498zi2vfa4v8vrd31rfw.png)

Great! Now how do we turn this into 3D objects?

Well, we can try to use SolidPython, a library that outputs OpenSCAD, a language for creating 3D objects.

```python
stroke_polygon = polygon(stroke_points)
stroke_obj = linear_extrude(10)(stroke_polygon)
```

If we do this for every stroke in the character, we get this output:

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/g0hqtnnrtrnd16we5dok.png">
  <figcaption>Every stroke in 好, linearly extruded</figcaption>
</figure>

Not bad. But we want to generate actual 3D objects here, not boring flat shapes.

We’ll need a way to split up each stroke into roughly equally sized parts.

## Cutting the character into parts

So we’ll need equidistant points for our next steps, and one point for each part of the stroke.

Luckily Svgpathtools makes it easy for us to sample the position of a point on a path given a number between 0 and 1:

```python
# create Svgpathtools Lines out of the medians
medians_lines = [Line(p1, p2) for p1, p2 in pairwise(medians)]
# create a Path out of those lines
medians_path = Path(*medians_lines)
# sample the Path
median_points = [medians_path.point(i) for i in np.linspace(0, 1, parts_count)]
```

If we now also plot those points in blue, we get this:

![Plot 好 2](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/96o033wipfrgt1dqzchn.png)

Now we now roughly how our strokes should be divided. But how do we actually calculate the shapes? We can use Voronoi diagrams for this.

![Plot 好 3](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/pn9uisigy9luus6no9yn.png) 

Because we want closed Voronoi regions (not regions stretching to infinity), we add 4 extra points around the shape, at a distance so they can’t interfere with the Voronoi diagram inside the center 1024x1024 square.

Zoomed out, that looks like this:

![Plot 好 4](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/2w2d5xnjbq71iz1923fq.png)

If we now intersect those polygons with the character, linearly extrude each part and move it according to the distance in the stroke, we get this:

![Linearly extruded character](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/2wopogn5fax7qcq0zbvp.png)

```python
obj += up(z)(
    linear_extrude(thickness)(
        intersection()(voronoi_polygon, strokes_polygon)
    )
)
```

Not terrible.

If we add some colors according to the position of the part in the stroke it looks like this:

![Colored linearly extruded character](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/b2otofkywo1v41tml5ao.png)

Now you might think we can increase the smoothness of the part by cutting it into more parts, but that won’t work in corners. 

We’ll do something different: we skew each part according to the slope.

## Shearing each part to generate a smooth surface

A skew in 3D is usually called a shear, and it can be represented with a 3D matrix called a shear matrix.

According to Wikipedia:

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/3qqlf8kt7jt8wycsfxfw.png">
  <figcaption>Source: <a target="_blank" href="https://en.wikipedia.org/wiki/Shear_matrix">https://en.wikipedia.org/wiki/Shear_matrix</a></figcaption>
</figure>

I don’t really understand what that means, but I know some matrix algebra. I can see that the 4th dimension of the input will influence the 1st dimension of the output.

If you reason about it, you can conclude: if only the 1st dimension is influenced, then the shear will be parallel to the axis of the 1st dimension: the x axis.

In our case, we need the shear to be parallel to the z dimension. The {% katex inline %}x{% endkatex %} and {% katex inline %}y{% endkatex %} inputs will be able to influence the {% katex inline %}z{% endkatex %} output.

{% katex %}
\begin{bmatrix}
1 & 0 & 0 & 0 \\\\
0 & 1 & 0 & 0 \\\\
a & b & 1 & 0 \\\\
0 & 0 & 0 & 1
\end{bmatrix}
{% endkatex %}

In this case, the output will be skewed as follows: {% katex inline %}z' = a*x + b*y + z{% endkatex %}

So we’ll need to calculate {% katex inline %}a{% endkatex %} and {% katex inline %}b{% endkatex %}.

That can’t be too hard. We already know the x, y and z coordinates of the center of each part of the stroke. The x and y coordinates of those are in blue on the charts above. The z coordinate we assign ourselves.

In our case, we want x and y to influence z. So the shear will be parallel to the z axis. Actually, the axis we want to shear is the axis between {% katex inline %}p_1{% endkatex %} and {% katex inline %}p_2{% endkatex %}.

We could calculate {% katex inline %}a{% endkatex %} and {% katex inline %}b{% endkatex %}, but we can also just rotate the object around the z axis so that the axis we want to shear is just the x axis.

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/sr0r0dlcvusjg48j17vl.png">
  <figcaption>Top view of some parts of a stroke</figcaption>
</figure> 

A rotation matrix looks as follows:

{% katex %}
\begin{bmatrix}
cos(\theta) & -sin(\theta) & 0 & 0 \\\\
sin(\theta) & cos(\theta) & 0 & 0 \\\\
0 & 0 & 1 & 0 \\\\
0 & 0 & 0 & 1
\end{bmatrix}
{% endkatex %}

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ms6ghui0fqo59nmcbg7r.png">
  <figcaption>Top view: We need to rotate each part by -θ so it aligns with the x axis</figcaption>
</figure>

```python
# calculate theta by using the arctangent
theta = np.arctan2(y2 - y1, x2 - x1)
# construct the rotation matrix, using -theta as the angle
rot_mat = np.matrix(
    (
        (cos(-theta), -sin(-theta), 0, 0),
        (sin(-theta), cos(-theta), 0, 0),
        (0, 0, 1, 0),
        (0, 0, 0, 1),
    )
).reshape((4, 4))
```

Because we aligned the shear axis with the x axis, {% katex inline %}b = 0{% endkatex %}. The shear matrix will now look as follows:

{% katex %}
\begin{bmatrix}
1 & 0 & 0 & 0 \\\\
0 & 1 & 0 & 0 \\\\
a & 0 & 1 & 0 \\\\
0 & 0 & 0 & 1
\end{bmatrix}
{% endkatex %}

{% katex inline %}a{% endkatex %} is the amount that {% katex inline %}z{% endkatex %} influences {% katex inline %}x{% endkatex %}. ({% katex inline %}x' = x + a*z{% endkatex %}) So we can reason: {% katex inline %}a{% endkatex %} needs to be proportional to the angle of the slope of the stroke. We’ll call this angle {% katex inline %}\alpha{% endkatex %}.

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/p86rie7pcrz7ze3hwaie.png">
  <figcaption>Side view: We need to shear each part by α</figcaption>
</figure>

We know the location of the points ({% katex inline %}p_1{% endkatex %} and {% katex inline %}p_2{% endkatex %}), so we can calculate the tangent as follows:

{% katex %}
distance_{xy} = \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}
\newline
tangent_\alpha = (z_2-z_1) / distance_{xy}
{% endkatex %}

```python
distance_xy = sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
tangent_alpha = (z2 - z1) / distance_xy
```

We actually don’t need to calculate {% katex inline %}\alpha{% endkatex %}! The tangent of {% katex inline %}\alpha{% endkatex %} is all we need to construct the shear matrix. This makes sense, the tangent is the ratio of {% katex inline %}z{% endkatex %} proportional to {% katex inline %}x{% endkatex %} and {% katex inline %}y{% endkatex %}.

Remember from SOH-CAH-TOA: Tangent = Opposite/Adjacent.

{% katex %}
\begin{bmatrix}
1 & 0 & 0 & 0 \\\\
0 & 1 & 0 & 0 \\\\
tangent_\alpha & 0 & 1 & 0 \\\\
0 & 0 & 0 & 1
\end{bmatrix}
{% endkatex %}

```python
shear_mat = np.matrix(
    (
        (1, 0, 0, 0),
        (0, 1, 0, 0),
        (tangent_alpha, 0, 1, 0),
        (0, 0, 0, 1)
    )
).reshape((4, 4))
```

Now we know the correct shear matrix.

But there’s something important we’re missing: our parts are not located at the origin, so the shear will also move the whole part depending on the distance to the origin. To solve this, we move the part to the origin first.

```python
translate_mat = np.matrix(
    (
        (1, 0, 0, -x1),
        (0, 1, 0, -y1),
        (0, 0, 1, 0),
        (0, 0, 0, 1)
    )
).reshape((4, 4))
```
 
<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/n63w6m1hao6kn8ap8fw8.png">
  <figcaption>Top view of part rotated at the origin</figcaption>
</figure>

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ww0ekhevqylo0subs8ue.png">
  <figcaption>Side view of part before shearing at the origin</figcaption>
</figure>

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/jpi445jza8094zittk57.png">
  <figcaption>Side view of part after shearing at the origin</figcaption>
</figure>

In total, the matrix will look like this:

```python
mat = np.identity(4)
# Apply translation
mat = np.matmul(translate_mat, mat)
# Apply rotation
mat = np.matmul(rot_mat, mat)
# Apply shear
mat = np.matmul(shear_mat, mat)
# Undo rotation
mat = np.matmul(np.linalg.inv(rot_mat), mat)
# Undo translation
mat = np.matmul(np.linalg.inv(translate_mat), mat)

total_obj += up(z)(multmatrix(np.asarray(mat))(obj))
```

We end up with something like this:

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/4qeo29gkmlw1eiu9tuix.png">
  <figcaption>The sheared result</figcaption>
</figure>

Pretty good!

![Character with sheared parts](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ruppunhycx2pfyp5fh3d.png)

Animated the shear:

{% embed https://www.youtube.com/embed/rASgv5nZfso?loop=1&playlist=rASgv5nZfso %}

As you can see, the shear in the corners and at the endpoints is not ideal yet. I partly solved this by smoothing the path with the points, and by taking the moving average of the shear.

# Smoothing

As you can see the smoothness in the corners is not great. I think it’s impossible to get perfect smoothness when using this part-splitting method, but we can make some improvements.

The problem is partly that the median points create Voronoi regions that are rather sharp in the corners. The sharper the corners in the median points, the sharper the region.

I added an algorithm to remove the sharpness of these corners.

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/1ysc93nd1ptcvg2ihmbj.png">
  <figcaption>Voronoi regions without path smoothing</figcaption>
</figure>

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/clskk301150z83i7ku0h.png">
  <figcaption>Voronoi regions with path smoothing</figcaption>
</figure>
 
I didn’t know any smoothing algorithms that could do this so I invented my own:

![Smoothing algorithm](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/zqwlw3jhtgtkhrf8ygt7.png)

- Create an empty list of output points
- Add the first input point to the list of output points
- For every 3-pair of input points:
    - Calculate the corner between the 3 points
    - Calculate the bisecting angle of the corner
    - Add a new point to the list of output points, with requirements:
        - The point is on the bisecting line
        - The point is on the inside of the corner
        - The sharper the corner, the higher the distance from the input point
- Add the last input point to the list of output points
- Re-interpolate points so they’re equidistant

Repeat this algorithm n times.

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/ja9yguw6yiakkpf7yjar.png">
  <figcaption>First iteration of the algorithm. The magnitude is in orange, the resulting point is in magenta.</figcaption>
</figure>

{% embed https://www.youtube.com/embed/nMTWrJPksAI?loop=1&playlist=nMTWrJPksAI %}
<figure>
  <figcaption>First step: 0 smoothing steps, last step: 10 smoothing steps</figcaption>
</figure>

{% embed https://www.youtube.com/embed/5YXa18Nz8fI?loop=1&playlist=5YXa18Nz8fI %}
<figure>
  <figcaption>First step: 0 smoothing steps, last step: 15 smoothing steps</figcaption>
</figure>

This algorithm converges to a straight line.

I also added a moving average of the shears, this helps with smoothing out the corners even more.

{% embed https://www.youtube.com/embed/nMTWrJPksAI?loop=1&playlist=nMTWrJPksAI %}
<figure>
  <figcaption>First step: remember 0 previous parts, last step: remember 10 previous parts</figcaption>
</figure>

Both are significant improvements to the smoothness of the part, and result in very usable 3D models of the parts.

These features also allow us to continuously increase the number of parts, as the corner issue has now been solved. We can increase the number of parts per stroke and achieve an almost smooth surface!

{% embed https://www.youtube.com/embed/HFSQVxjATzA?loop=1&playlist=HFSQVxjATzA %}
<figure>
  <figcaption>The end result</figcaption>
</figure>

## Flattening the part

We can also add a mode to extend each stroke to the bottom:

<figure>
  <img src="https://dev-to-uploads.s3.amazonaws.com/uploads/articles/nayb3lkbvsgizy99jhoz.png">
  <figcaption>To bottom mode</figcaption>
</figure>

But for printing a ring or something mostly flat, this is much too tall. It would be better if we could “flatten” it before stretching it to the bottom. We can achieve this by rotating the top surface so it’s more or less flat.

I called this untilted mode.

![Untilting](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/r7k31ujoy5euklb3f4l5.png)

It works by doing Principal Component Analysis on the point cloud of the center point of each part. We can then undo the tilt by matrix-multiplying the eigenvectors of the PCA with the object.

For example:

```python
from sklearn.decomposition import PCA

pca = PCA()
pca.fit(medians_3d)
eigenvectors = pca.components_

mat = [(*values, 0) for values in eigenvectors] # pad with 0

obj = multmatrix(mat)(obj)
```

If we would just run this, we can see that the top part is flat, as we want. But the character is rotated and skewed as well.

![Before unskewing](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/o1tfxk42xlb84yy8xuk2.png) 

{% embed https://www.youtube.com/embed/mzu3xy6xKBo?loop=1&playlist=mzu3xy6xKBo %}
<figure>
  <figcaption>Untilted mode before unskewing</figcaption>
</figure>

We need to add some extra code to force the transformation to keep vertical lines vertical.

```python
# make sure z doesnt have an effect on x and y
eigenvectors[0][2] = 0
eigenvectors[1][2] = 0
```

The eigenvectors could also be inverted for some reason I’m not sure of.

```python
# make sure all eigenvectors are in the same direction as the current axis
for i in range(3):
    # e.g. for Z, if 0, 0, 1 would map to *, *, < 0, then invert it
    if eigenvectors[i][i] < 0:
        eigenvectors[i] = -eigenvectors[i]
```

We also need to remove possible rotation around the z axis and scaling of the xy plane.

```python
# make sure there's no rotation around z or scaling of the xy plane
eigenvectors[0][0] = 1
eigenvectors[0][1] = 0
eigenvectors[1][0] = 0
eigenvectors[1][1] = 1
```

After these checks, the character is as we want:

![Skew correction](https://dev-to-uploads.s3.amazonaws.com/uploads/articles/6qoo4b6fqioijj2ss1i5.png) 

# Future

I think it should be possible to generate perfectly smooth surfaces by transforming lattices on a mesh. This might be possible with a Blender plugin of some sorts. I think it would even be possible to reuse the smoothing algorithm, whether or not the Voronoi diagrams could be used remains to be seen. For now, this simple script will do for most 3D printing needs.

OpenSCAD rendering and rendering to STL is a bottleneck as well. I think it should be possible to port this to a web based application with sliders for the parameters, so you would get instant feedback instead of having to wait for the rendering.

Porting this to Javascript using OpenJSCAD, or to a compiled language and then to WebAssembly is also a possibility.

# Interested?

If you want to try generating your own part, check out [the repo](https://github.com/hansottowirtz/3d-hanzi-generator), and if you have any ideas or suggestions, open an issue or send me a message on [Twitter](https://twitter.com/hansottowirtz).

{% github https://github.com/hansottowirtz/3d-hanzi-generator %}