Fix the bug discovered by Victor Ostromoukhov. Due to the bug some qu…

…aternions which lead to bijective digitized rotations were not certified as such. This was because the old condition to check if the sytem Ax = y is feasible was not sufficient for the integer programming problem.

QuaternionCertification.m

@@ -1,4 +1,4 @@
-(* :Title: QuaternionCertification, Version 1.0, 2016*)
+(* :Title: QuaternionCertification, Version 1.1, 2017*)
 
 (* :Author: Kacper Pluta, [email protected]
 	    Laboratoire d'Informatique Gaspard-Monge - LIGM, A3SI, France 
@@ -9,22 +9,28 @@
   A package which implements algorithm which allows to verify if a given Lipschitz quaternion induce
   bijective digitized rotations. The algorithm was introduced in:
 
-@inproceedings{Pluta:CTIC:2016,
-	author = {Pluta, K. and Romon, P. and Kenmochi, Y. and Passat, N.},
-	booktitle = {{CTIC}},
-	publisher = {Springer},
-	series = {{Lecture Notes in Computer Science}},
-	title = {{Bijectivity certification of 3D digitized rotations}},
-	year = {Forthcoming 2016}
+@Inbook{Pluta:CTIC:2016,
+  author="Pluta, Kacper and Romon, Pascal and Kenmochi, Yukiko and Passat, Nicolas",
+  editor="Bac, Alexandra and Mari, Jean-Luc",
+  title="Bijectivity Certification of 3D Digitized Rotations",
+  bookTitle="Computational Topology in Image Context: 6th International Workshop, CTIC 2016,
+             Marseille, France, June 15-17, 2016, Proceedings", 
+  year="2016",
+  publisher="Springer International Publishing",
+  address="Cham",
+  pages="30--41",
+  isbn="978-3-319-39441-1",
+  doi="10.1007/978-3-319-39441-1_4",
+  url="http://dx.doi.org/10.1007/978-3-319-39441-1_4"
 }
 
 *)
 
 (* :Context: QuaternionCertification` *)
 
-(* :Package Version: 1.0 *)
+(* :Package Version: 1.1 *)
 
-(* : Copyright: Kacper Pluta, 2016
+(* : Copyright: Kacper Pluta, 2017
   All rights reserved.
 
  Redistribution and use in source and binary forms, with or without
@@ -58,7 +64,7 @@
 *)
 
 
-(* :Mathematica Version: 10.0.0 *)
+(* :Mathematica Version: 11.0.0 *)
 
 BeginPackage["QuaternionCertification`"]
 
@@ -105,25 +111,75 @@
 
 (*
 
-Procedure: QuaternionRgionMultiplication
+Procedure: SystemA
 
 Summary:
-  This procedure builds a system [A|B] (see paper cited in :Summary: section of this package).
+  This procedure builds the matrix A from the system (see paper cited in :Summary: section of this
+  package).
 
 Parameters:
   a, b, c, d -- a quaternion represented by four integers.
-  bb1, bb2, bb3, bb4 -- a four dimensional integer point.
 
 Output:
-  The augment matrix [A|B].
+  The matrix A.
 
 *)
-SystemAB[a_, b_, c_, d_, bb1_, bb2_, bb3_, bb4_] := 
+SystemA[a_, b_, c_, d_] := 
   Return[ {
-   { b, c, d, -b, -c, -d, bb1 }, { a, -d, c, -a, -d, c, bb2 },
-   { d, a, -b, d, -a, -b, bb3 }, { -c, b, a, -c, b, -a, bb4 } 
+   { b, c, d, -b, -c, -d }, { a, -d, c, -a, -d, c },
+   { d, a, -b, d, -a, -b }, { -c, b, a, -c, b, -a } 
   } ];
-(* end of SystemAB *)
+(* end of SystemA *)
+
+
+(*
+
+Procedure: CheckMultiplicity
+
+Summary:
+  This procedure checks of elements of the second list are multiples of the respective elements of
+  the first list
+
+Parameters:
+  Two lists of integers which are of the same length.
+
+Output:
+  True if all the elements of the second lists are integer multiples of the elements of the first
+  list.
+
+*)
+CheckMultiplicity[ListA_, ListB_] := Module[{result},
+  If[Length[ListA] != Length[ListB], Return[False],
+    result = True;
+    Do[ result = result && IntegerQ[ ListB[[x]] / ListA[[x]] ], {x, 1, Length[ListA]} ];
+    Return[result];
+  ];
+];
+(* end of CheckMultiplicity *)
+
+
+(*
+
+Procedure: HasIntegerSolutions
+
+Summary:
+  This procedure checks if y is a solution to the system Ax = y (see the paper cited in :Summary:)
+  To check the feasibility of the system Ax = y we use a similar condition to the one give by
+  A.  Schrijver in Theory of linear and integer programming (page 51)
+
+Parameters:
+  A Smith Normal Form of the matrix A and a four dimensional integer point y
+
+Output:
+  True if Ax = y has integer solutions.
+
+*)
+HasIntegerSolutions[SNF_, y_] := Module[{yp, diag},
+  diag = Diagonal[ SNF[[2]] ];
+  yp = SNF[[1]].y;
+  If[ yp[[-1]] != 0, Return[False], Return[ CheckMultiplicity[ diag[[;;-2]], yp[[;;-2]] ] ] ];
+];
+(* end of HasIntegerSolutions *)
 
 
 (*
@@ -140,29 +196,28 @@
   True if q leads to 3D bijective digitized rotation and False otherwise.
 
 *)
-CertifyQuaternion[q_] := Module[ {p, pp, cube1, cube2, bounds, y, h, i, j, k, pivot47},
+CertifyQuaternion[q_] := Module[ {p, pp, cube1, cube2, bounds, y, h, i, j, k, SNF},
   If[MemberQ[Map[IntegerQ, q], False], Throw["This is not a Lipschitz quaternion."]];
   p = q / GCD[ q[[1]], q[[2]], q[[3]], q[[4]] ];
   pp = { p[[1]], -p[[2]], -p[[3]], -p[[4]] };
   cube1 = TransformedRegion[OpenUnitCube[], QuaternionRegionMultiplication[ pp, { Indexed[#, 1],
   Indexed[#, 2], Indexed[#, 3], Indexed[#, 4] } ] &];
   cube2 = TransformedRegion[OpenUnitCube[], QuaternionRegionMultiplication[ { Indexed[#, 1],
   Indexed[#, 2], Indexed[#, 3], Indexed[#, 4] }, pp] & ]; 
-  bounds = Round[RegionBounds[cube1]];  
+  bounds = Round[RegionBounds[cube1]];
+  SNF = SmithDecomposition[ SystemA[ p[[1]], p[[2]], p[[3]], p[[4]] ] ];
+  (*Compute GCDs for minors of the all valid dimensions*)
   Do[
      y = {h, i, j, k};
      (*If zero the y times q goes to R^3*)
      If[ h * p[[1]] - i * p[[2]] - j * p[[3]] - k * p[[4]] != 0, Continue[]];
-     pivot47 = HermiteDecomposition[ SystemAB[p[[1]], p[[2]], p[[3]], p[[4]], y[[1]], y[[2]], y[[3]],
-                                              y[[4]]]][[2]][[4]][[7]]; 
-     (*System has solution iff the last pivot is zero.*)
-     If[y \[NotElement] cube1 || pivot47 != 0, Continue[]];
+     If[y \[NotElement] cube1 || !HasIntegerSolutions[ SNF, y ], Continue[]];
      If[y \[NotElement] cube2, Goto[end]];
    , {h, bounds[[1]][[1]], bounds[[1]][[2]]},
    {i, bounds[[2]][[1]], bounds[[2]][[2]]},
    {j, bounds[[3]][[1]], bounds[[3]][[2]]},
    {k, bounds[[4]][[1]], bounds[[4]][[2]]}
- ];
+  ];
   (* Calling just Return[] in the If[] will take us at the end where is the second Return and output
     will be True. *) 
   Return[True];