I'm trying to plot the Minkowski sum of a square and a rhombus.
# a unit square
sq <- square(1)
# a rhombus
line1 <- psp(0,1,1,0,window=owin())
line2 <- psp(0,0,1,1,window=owin())
rhombus <- line1 %(+)% line2
# the Minkowski sum of the rhombus and the square
ply <- rhombus %(+)% sq
The MinkowskiSum() function (or the %(+)% operator) returns an object with attribute $bdry that contains 2 different sets of boundary points.
ply$bdry # it has two sets of boundary points
After plotting it out, it seems like the first element in ply$bdry is the correct set of boundary points for ply.
par(mfrow=c(1,2))
plot(ply, axes=TRUE, border="blue", lwd=2, hatch=TRUE, hatchargs=list(col="blue"))
plot(rhombus, add=TRUE, axes=TRUE, lwd=2, hatch=TRUE, hatchargs=list(texture=5))
plot(sq, add=TRUE, border="red", lwd=2)
plot(ply, axes=TRUE, border="white")
polygon(ply$bdry[[1]]$x, y = ply$bdry[[1]]$y, border="lightpink", col="lavenderblush") # this is the correct set of boundary points of the Minkowski sum
polygon(ply$bdry[[2]]$x, y = ply$bdry[[2]]$y, border="blue", col="lavender") # this is not the correct set of boundary points of the Minkowski sum
I've read through the manual and the source code but still couldn't figure out:
Any help is appreciated!
You don't need to inspect the internal structure of the resulting object to decide how to plot it. Just use plot (which will be dispatched to plot.owin).
The result ply is a spatial window (object of class "owin"). Internally this is represented as a list with many components. One of the components is $bdry which gives the coordinates of the vertices of the polygon or polygons that represent the object. A spatial window may consist of multiple separate polygons, or polygons with holes; in these cases, $bdry is a list with one entry for each connected polygon representing a piece of the outer boundary or the boundary of a hole.
In your example, the result of the operation consists of two polygons. The first one is the one you would expect. The second one is a very tiny hole. You can see this using print.owin or summary.owin:
summary(ply)
Window: polygonal boundary
2 separate polygons (1 hole)
vertices area relative.area
polygon 1 8 7.00000e+00 1.00e+00
polygon 2 (hole) 3 -9.75782e-19 -1.39e-19
enclosing rectangle: [0, 3] x [0, 3] units
Window area = 7 square units
The tiny hole is the result of numerical rounding error. Polygon calculations in spatstat are performed mainly using the polyclip package which works in integer arithmetic. Input polygons in continuous space are discretised to integer coordinates with a default resolution of 1e-09 units, using the midpoint of the surrounding rectangle as the origin. The polyclip calculation results in a hole of size about 1e-09 by 1e-09. Use the ... arguments to MinkowskiSum to change the discretisation.