Geometry visualization issue after clip operation and rotation

[FilipeFcp]

After subtracting one geometry from another of the same size along a given dimension and rotating the object, a one-pixel artifact appears in the plot:

To Reproduce

 import tidy3d as td

box  = td.Box(center = (0,-9.01,10),
              size = (10,0.01,10))
cylinder = td.Cylinder(radius = 9,
                       center = (0,-9.01,0),
                       length=0.01,
                       axis = 1)
ax = (box-cylinder).plot(x=0)
ax.set_aspect('auto')
ax = (box-cylinder).rotated(np.pi/2,axis = 0).plot(x=0)
ax.set_aspect('auto') 

Expected behavior
Adding a small extra value for the length argument of the cylinder produces the expected behavior: length=0.01+ 0.01

[jewettaijfc]

Why this happens:

• The output geometry is different when the intersection of two shapes is rotated. It's not just a visual glitch. The second shape really has an (extremely thin) triangle protruding from it.
• This happens because we are rotating (transforming) the coordinates of the shapes before we compute their intersections. (Relevant code is here.)
• Consequently, floating-point roundoff error due to rotation sometimes causes tiny distortions in the overlapping regions of the two shapes. This prevents them from being subtracted correctly.
• One way to solve this problem would be to switch the order: perform the shape intersection calculations before the rotation in the plane.
• Unfortunately, there is no easy way to do that. (DETAILS: This is because shapely can only compute the intersection between 2D shapes (not 3D shapes). To make each geometric shape in our simulation 2D, it must first be projected onto the image plane. This is accomplished using the to_2D array. Unfortunately, this array also contains rotation information. There's no way to only project onto the plane without also rotating the object. So, if we are using shapely, then we have to apply the to_2D transformation (including rotations) first before computing the shape intersections.)

---

Note to self:

FEEL FREE TO SKIP. These are my own notes.

• When 3D volumetric shapes (eg: box, cylinder) are plotted in 2D, Tidy3D computes their intersection with the display plane (eg. YZ axes) (using the Geometry.intersections_plane() function, which is invoked by Geometry.plot()).
• The plots look different because the resulting shapes being displayed actually are different in the rotated and unrotated cases.

Unrotated case: (box-cylinder).plot(x=0)

• I checked the list of coordinates generated by that function (stored in the shapes_intersect variable).
• In the 1st plot example ((box-cylinder).plot(x=0)), the shapes_intersect variable stores the coordinates of a simple rectangle (with only 4 unique points at the rectangle's corners).

shapely.get_coordinates(shapes_intersect[0])=array([
       [-9.015, 15.   ],
       [-9.005, 15.   ],
       [-9.005,  9.   ],
       [-9.015,  9.   ],
       [-9.015, 15.   ]

Rotated case: (box-cylinder).rotated(np.pi/2,axis = 0).plot(x=0)

• In the 2nd plot example the shapes_intersect variable has a couple extra vertices, including one located at -5,-9.015, which is responsible for the infinitely thin line we are seeing that extends to the right.

shapely.get_coordinates(shapes_intersect[0])=array([
       [ -5.   ,  -9.015],  # <- problematic
       [-10.   ,  -9.015],
       [-15.   ,  -9.015],
       [-15.   ,  -9.01 ],
       [-15.   ,  -9.005],
       [-10.   ,  -9.005],
       [ -9.   ,  -9.005],
       [ -9.   ,  -9.015],
       [ -5.   ,  -9.015],  # <- problematic
       [ -5.   ,  -9.015],  # <- problematic
])

[jewettaijfc]

The thin line that we don't want protruding from the object is actually a very thin triangle.

Possible Solution

Discard vertices on the polygon boundary if their absence does not affect the area of the polygon.

Question

Am I allowed to do that? Are objects allowed to be arbitrarily thin? Is it okay to discard extremely thin triangles from shapely objects?

[tylerflex]

Hi @jewettaijfc I believe we disallow polygons with small areas (near numerical precision). But maybe the area here was large enough to make the cut.

https://github.com/flexcompute/tidy3d/blob/develop/tidy3d/components/geometry/polyslab.py#L49

@momchilmm @lucas-flexcompute do you know how we should handle this?

[jewettaijfc]

Hi Tyler!
To clarify: Discarding small polygons is not the same as discarding small triangles that are part of a larger polygon. That's what I'm asking for permission to do.

[jewettaijfc]

(In this case the width of the triangle was approximatley 1e-15, which I suspect is pretty close to 64-bit precision.)
Thanks for your fast reply, @tylerflex !

I am wondering if there is some way to represent infinitely thin charged sheets or current sheets in Tidy3D. (?) If so, if I start deleting thin triangles, I worry I might accidentally erase those.

[tylerflex]

ok looks like fp_eps here is actually 1e-7 or so. so something else is going on. Maybe the other people can comment, it's been a while since I worked with this code.

I am wondering if there is some way to represent infinitely thin charged sheets or current sheets in Tidy3D. (?) If so, if I start deleting thin triangles, I worry I might accidentally erase those.

We do support these, but they are evaluated I think as lines in that plot_shape method.

[jewettaijfc]

Great! Then I will go ahead and remove those super-thin triangles from each polygon.

I started a new PR which uses this approach.
(Warning: It is not yet complete.)

(Currently, I only bother to do this for clipped polygons. In other words, polygons that were created by ClipOperation.intersections_plane() or ClipOperation.intersections_tilted_plane(). If you want me do this consistently for all polygons that appear in plots, let me know.)