contact: amantzafmath.uoa.gr & grad0961di.uoa.gr

Implementing python code for the computation/approximation of the decomposition of a convex lattice polygon P as P=A+B, where "+" denotes Minkowski addition.

Tools:

References:

- I.Z. Emiris, E. Tsigaridas, Minkowski decomposition of convex lattice polygons, "Algebraic geometry and geometric modeling", 2005
- S. Gao and A.G.B. Lauder, Decomposition of polytopes and polynomials, "Discrete and computational geometry", 2001
- A. Tuzikov, H.J.A.M. Heijmans, Minkowski decomposition of convex polygons into their symmetric and asymmetric parts, "Pattern Recognition Letters", 1998

Problem MINKOWSKI statement:

- Input: A sequence
of Vectors_2 defing the normal
fan of a convex lattice polygon P with n vertices

(or a sequence of Points_2 defining the vertices of P) - Output: All distinct sequences of Vectors_2 defining fans of summands A, where P=A+B

The normal fan of P is the set of vectors perpendicular to the edges u_i = p_{i+1} - p_i , i=0,...,n and of length equal to ||u_i||, where p_i 's are the vertices of P and p_n=p_0 by convention.

The polygon P is fully defined by a vertex and the set of it's edges. For the purpose of decomposition translations are not important, so the actual size of input is determined by the bitsize of the edges (or the fan) of P rather than the bitsize of the initial polygon vertices.

A subvector of a polygon edge u_i=(u_ix, u_iy) is of the form je_i, where j is an integer 0<j<=d_i and d_i=gcd(u_ix, u_iy). The vector e_i= u/d_i is the primitive edge of u.

The sum of all integers j for all edges u_i equals the number of integer points on the perimeter of the polygon (also known by Pick's formula) and is related to the complexity of the problem.

a) Define the primitive edge sequence. That is, for each primitive edge e_i=(e_x, e_y) consider the O(log d_i) vectors that are multiples of e_i with scalar factors all the powers of 2^0, 2, 2^2,... such that 2^r <= d_i/2, where d_i=gcd(u_ix, u_iy). Note that the sum of these integers does not exceed d_i (by geometric series sum). We dublicate the nessesary powers so that all these integers (resp. vectors) sum to d_i (resp. u_i).

Observe that these vectors suffice to represent any subvector of every edge of P, since every integer 1,...,d_i has a representation in the binary arithmetic system using a subset of the scalars defined before. In other words, the primitive edge sequence is a decomposition of the edges of P that generates (by addition) all subvectors of edges of P.

b) Give the set of all these vectors, that is of cardinality O(n log D) , D=max{d_i} as an input to the 2D-subset sum problem:

Input: A set of integer vectors in the plane

Output: One (or more) proper subset of the vectors that sums to (0,0).

Note that the input we created is polynomial in the bitsize of P, which is the input of the original MINKOWSKI problem.

c) How can we solve the 2D-subset sum problem? The problem can be easily reduced to 1D-subset sum, by considering for every (u_i,v_i) an integer u_i + Lv_i, where L is a sufficiently large integer so that it seperates the vanishing of a sum of u_i's from the vanishing of the corresponding v_i's.

The integer L should be as small as possible, since it affects the pseudopolynomial complexity factor of the subset sum problem.

The usual (1D) subset sum problem is:

Input: A set of integers

Output: One (or more) proper subset of the integers that sums to 0.

We initially attempted a direct algorithm for solving the 2D-subset sum problem. This used a 3D array of Vectors_2 and has been proven inefficient compared to reduction to the usual subset sum.

d) Solve the 1D-subset sum problem by dynamic programming and transform the solution back to MINKOWSKI output.

These ideas lead to a pseudopolynomial algorithm for MINKOWSKI of complexity O(n LV log D) where V=max{v_i}.

The notion of an approximate solution for MINKOWSKI is to be explored.

- A general graphical interface using Python-CGAL-visual for the graphical input/output of polygons and decompositions. This includes an integer grid along with functions for input of integer convex polygons, functions for drawing a lattice polygon on the grid (having as input either the edge sequence or the fan of the polygon) and other graphics-related routines.

- A set of auxiliary routines that pass from edge representation to fan representation, compute the complementary summand by Minkowski substraction, compute the primitive edge sequence and several other auxiliary tasks.

- A function for the enumeration of all decompositions by exhaustive search.

- A dynamic programming algorithm for the direct solution of the 2D-subset sum problem (this has been proven inefficient compared to reduction to the usual subset sum).

- A dynamic programming algorithm for the usual subset sum problem along with routines for the reduction of MINKOWSKI input to subset sum input as well as transformation of the subset sum output to MINKOWSKI output.

- A test routine that generates the polygons used in our experiments below, executes the algorithms on them and returns results and timings.

Two-summands:

One can find the complementary summand for each solution by Minkowski substraction:

Three-summand: Complementary:

Four-summands:

Complementaries

Our Dynamic programming implementation does not return the totality of possible summands (there are exponentially many of them after all), but a subset of them, ranging from small (i.e. line summands, if they exist) to large ones. All solutions could be explored by a careful backtracking, but our intention is to discover a small amount of "interesting" decompositions rather than output an exponentially large set. In this example Dynamic programming outputs 3 decompositions: One containing line summand, one with a four-summand and one with a three-summand.

An interesting task would be to define a compact representation of all possible summands. For example instead of outputting 4 line summands, we could denote that the two parallel edges define 4 decompositions by outputting the integer range [1,4] along with pointers to these two edges. We expect that this technique is generalizable and may lead to new algorithms for MINKOWSKI.

We consider scalar multiples of these polygons, by a scalar factor k, in order to create large instances (and also hard: with a plethora of integer points on the perimeter).

For every polygon and for k ranging from 20 to 400 we list the number of decompositions discovered ("#") and the time needed ("t") in seconds on a Intel Core2Duo 2GHz machine with 1GB of main memory running Ubuntu 7.10 and Python 2.5.

Dynamic Programming |
||||||||||||||||

k | # | t | # | t | # | t | # | t | # | t | # | t | # | t | # | t |

20 | 5 | 0.011 | 11 | 0.016 | 11 | 0.060 | 11 | 0.029 | 11 | 0.073 | 17 | 0.133 | 17 | 0.091 | 23 | 0.273 |

40 | 6 | 0.056 | 13 | 0.072 | 13 | 0.278 | 13 | 0.136 | 13 | 0.339 | 20 | 0.688 | 20 | 0.414 | 27 | 1.290 |

60 | 8 | 0.171 | 17 | 0.210 | 17 | 0.864 | 17 | 0.427 | 17 | 0.991 | 26 | 1.903 | 26 | 1.300 | 35 | 3.812 |

80 | 7 | 0.284 | 15 | 0.358 | 15 | 1.272 | 15 | 0.671 | 15 | 1.670 | 23 | 2.889 | 23 | 1.951 | 31 | 5.803 |

100 | 8 | 0.406 | 17 | 0.534 | 17 | 2.232 | 17 | 1.092 | 17 | 2.747 | 26 | 5.102 | 26 | 3.300 | 35 | 10.004 |

120 | 9 | 0.709 | 19 | 0.910 | 19 | 3.750 | 19 | 1.950 | 19 | 4.631 | 29 | 8.478 | 29 | 5.369 | 39 | 16.914 |

140 | 9 | 0.936 | 19 | 1.219 | 19 | 5.087 | 19 | 2.457 | 19 | 6.228 | 29 | 11.607 | 29 | 7.720 | 39 | 22.329 |

160 | 8 | 1.036 | 17 | 1.364 | 17 | 6.212 | 17 | 2.887 | 17 | 7.766 | 26 | 14.274 | 26 | 9.006 | 35 | 26.784 |

180 | 10 | 1.620 | 21 | 2.238 | 21 | 9.398 | 21 | 4.429 | 21 | 11.725 | 32 | 20.913 | 32 | 13.515 | 43 | 40.139 |

200 | 9 | 1.938 | 19 | 2.427 | 19 | 11.031 | 19 | 5.289 | 19 | 13.398 | 29 | 23.769 | 29 | 15.142 | 39 | 45.751 |

220 | 11 | 2.765 | 23 | 3.595 | 23 | 14.869 | 23 | 7.185 | 23 | 18.192 | 35 | 32.912 | 35 | 21.688 | 47 | 65.636 |

240 | 10 | 2.871 | 21 | 3.808 | 21 | 16.155 | 21 | 7.813 | 21 | 19.798 | 32 | 36.619 | 32 | 23.625 | 43 | 71.561 |

260 | 9 | 3.080 | 19 | 4.072 | 19 | 17.159 | 19 | 8.363 | 19 | 21.089 | 29 | 38.207 | 29 | 25.282 | 39 | 76.175 |

280 | 10 | 3.903 | 21 | 5.159 | 21 | 22.091 | 21 | 10.669 | 21 | 26.980 | 32 | 48.893 | 32 | 32.218 | 43 | 97.786 |

300 | 11 | 4.968 | 23 | 6.494 | 23 | 27.730 | 23 | 13.519 | 23 | 35.191 | 35 | 61.290 | 35 | 40.402 | 47 | 135.049 |

320 | 9 | 6.306 | 19 | 6.879 | 19 | 26.163 | 19 | 12.704 | 19 | 32.093 | 29 | 57.906 | 29 | 38.283 | 39 | 127.395 |

340 | 11 | 6.344 | 23 | 10.958 | 23 | 45.373 | 23 | 17.638 | 23 | 43.604 | 35 | 90.129 | 35 | 54.920 | 47 | 175.801 |

360 | 11 | 7.220 | 23 | 9.363 | 23 | 39.985 | 23 | 19.497 | 23 | 58.086 | 35 | 88.505 | 35 | 58.273 | 47 | 227.275 |

380 | 13 | 9.661 | 27 | 12.239 | 27 | 61.804 | 27 | 25.500 | 27 | 75.913 | 41 | 126.344 | 41 | 84.147 | 55 | 380.184 |

400 | 10 | 8.519 | 21 | 10.589 | 21 | 45.117 | 21 | 21.997 | 21 | 55.252 | 32 | 110.953 | 32 | 68.678 | 43 | 273.083 |

In order to compare with exhaustive search (that returns all possible summands), we list the timings for k=1,...,7 (larger instances cannot be exhaustively searched within reasonable time):

Exhaustive Search |
||||||||||||||||

k | # | t | # | t | # | t | # | t | # | t | # | t | # | t | # | t |

1 | 1 | 0.000 | 3 | 0.000 | 3 | 0.000 | 3 | 0.000 | 3 | 0.001 | 7 | 0.002 | 9 | 0.002 | 23 | 0.007 |

2 | 2 | 0.000 | 8 | 0.001 | 8 | 0.001 | 8 | 0.001 | 12 | 0.004 | 34 | 0.014 | 44 | 0.012 | 236 | 0.121 |

3 | 3 | 0.001 | 15 | 0.004 | 15 | 0.004 | 15 | 0.004 | 27 | 0.132 | 99 | 0.065 | 135 | 0.065 | 1263 | 1.325 |

4 | 4 | 0.002 | 24 | 0.009 | 24 | 0.009 | 24 | 0.009 | 54 | 0.047 | 238 | 0.311 | 324 | 0.255 | 4724 | 7.601 |

5 | 5 | 0.003 | 35 | 0.018 | 35 | 0.018 | 35 | 0.018 | 91 | 0.119 | 479 | 0.869 | 665 | 0.761 | 13943 | 32.867 |

6 | 6 | 0.005 | 48 | 0.034 | 48 | 0.034 | 48 | 0.033 | 146 | 0.254 | 884 | 2.006 | 1224 | 2.043 | 34936 | 114.496 |

7 | 7 | 0.009 | 63 | 0.058 | 63 | 0.057 | 63 | 0.058 | 215 | 0.501 | 1487 | 4.398 | 2079 | 4.387 | 77503 | 352.390 |

Observe (from the fourth column of the above table or by an easy combinatorial argument) that a square polygon of side k has exactly k^2+2k -1 decompositions. This number is incremented by one in the table above, because the trivial decomposition is counted. In fact, exhaustive search will count each distinct decomposition twice (once for each summand), so there are ceil[ (k^2 +2k-1)/2 ] distinct decompositions of a square polygon. The "ceil" comes from the fact that when k in even there exists a "self-dual" decomposition (to two equal squares of side k/2) counted only once.

Similar formulas for the number of decompositons occur for the rest of the canonical polygons.