bugGNU Octave - Bugs: bug #34604, Trivial delaunay throws qhull error

 
 

bug #34604: Trivial delaunay throws qhull error

Submitter:  Ben Abbott <bpabbott>
Submitted:  Thu 20 Oct 2011 03:21:38 PM UTC
   
 
Category:  None Severity:  4 - Important
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Ben Abbott Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 25 Oct 2011 05:40:41 PM UTC, comment #11: 

Rik,

Thanks for taking care of this. I'm still working on my patch idea.


Ben Abbott <bpabbott>
Group Member
Tue 25 Oct 2011 05:28:45 PM UTC, comment #10: 

Ben, I checked in a patch that fixes your problem by adding 'Qz' to the list of Qhull options for low-dimensional data.  From all that I understand about Qhull this will not a concern.  Just to be sure, I ran multiple sets of random vectors through Qhull with and without the Qz option and the Delaunay triangulation is always the same.

The changeset is here (http://hg.savannah.gnu.org/hgweb/octave/rev/7ff0bdc3dc4c)

Rik <rik5>
Group administrator
Fri 21 Oct 2011 06:32:10 PM UTC, comment #9: 

Rik,

I've been looking at improving the compatibility of the patch() command. My thought was to mesh patches into triangles to allow them to have interpolated colors. I was going to use delaunay() to mesh convex patches.

I don't expect to finish this effort this week.

Please include the addition of "Qz" to the defaults when you push your changeset.

I don't know enough about the internals of qhull to have an opinion on whether or not not "Qz" can cause problems in some cases. Thus, I'm inclined to be conservative and only include it when precision error result (as jwe suggested). Is that possible?

Ben Abbott <bpabbott>
Group Member
Fri 21 Oct 2011 05:44:31 PM UTC, comment #8: 

Is there any possible harm in adding Qz by default, or should it only be added if there is a failure without it?

John W. Eaton <jwe>
Group administrator
Fri 21 Oct 2011 05:32:35 PM UTC, comment #7: 

Yes to adding 'Qz' to the defaults.  I am working on what has grown to be an enormous patch for all the geometry functions.  Octave was not cleaning up memory correctly which was why 'clear all' was needed.  I've changed a bunch of other things to clean up the code as well.  I expect by the end of the weekend to check something in.

Rik <rik5>
Group administrator
Fri 21 Oct 2011 03:34:02 AM UTC, comment #6: 

Rik,

What do you think, should the "Qz" option be added by default?

Ben Abbott <bpabbott>
Group Member
Thu 20 Oct 2011 08:23:45 PM UTC, comment #5: 

The error message does appear to be accurate.  The center point of the square is exactly on a facet line (the diagonal from the corners of the square).  The warning message suggests using two options--both of which work for me.

The first option is 'QJ' which "Joggles" the input points.  This does eliminate the exact co-planarity.  However, its use isn't particularly recommended because it reduces the accuracy of the triangulation.

The second option is 'Qz' which adds a point at infinity.  This option seems to work well.


y = [0 1 1 0];
x = [0 0 1 1];
tri = delaunay (x,y, 'Qz')
tri =

   4   2   3
   4   2   1


Perhaps this option should always be added to low-dimensional data?  The relevant code is in src/__delaunayn__.cc:89.


  // default options
  if (dim <= 3)
    options = "Qt Qbb Qc";
  else
    options = "Qt Qbb Qc Qx";


Interestingly, I tried the same trial figure in three dimensions, i.e., a unit cube and delaunay3 fails in exactly the same way.  By adding 'Qz' I can get delaunay3 to pass as well. 

Rik <rik5>
Group administrator
Thu 20 Oct 2011 07:10:07 PM UTC, comment #4: 

opps, I had a typo that swallowed the rest of my comment.


y = [0 1 1 0];
x = [0 0 1 1];
tri = delaunay ([x 0.5], [y 0.5])
tri =

  5   2   1
  5   2   3
  4   5   1
  4   5   3


Ben Abbott <bpabbott>
Group Member
Thu 20 Oct 2011 07:08:50 PM UTC, comment #3: 

In case it might give some hint as to what is wrong, I Just noticed that adding a point to the interior resolves the problem.


y = [0 1 1 0];
x = [0 0 1 1];
tri = delaunay ([x 0.5], [y 0.5])
tri =

   5   2   1
   5   2   3
   4   5   1
   4   5   3
-verbatiim-

Ben Abbott <bpabbott>
Group Member
Thu 20 Oct 2011 06:58:22 PM UTC, comment #2: 

Thanks Rik.

I've changed the OS field to "any".

Ben Abbott <bpabbott>
Group Member
Thu 20 Oct 2011 06:55:11 PM UTC, comment #1: 

The same error message appears for a Linux build (Kubuntu 10.04) and qhull version 2009.1-1.  The same problems with further delaunay operations are also confirmed.  This is with a development tip from 10/20/2011.

Rik <rik5>
Group administrator
Thu 20 Oct 2011 03:21:38 PM UTC, original submission:  


y = [0 1 1 0];
x = [0 0 1 1];
tri = delaunay (x, y)
QH6154 qhull precision error: initial facet 1 is coplanar with the interior point
ERRONEOUS FACET:
- f1
    - flags: bottom simplicial upperDelaunay flipped
    - normal:   -0.4082  -0.4082   0.8165
    - offset:         -0
    - vertices: p1(v2) p2(v1) p0(v0)
    - neighboring facets: f2 f3 f4

While executing:  | qhull d Qt Qbb Qc
Options selected for Qhull 2010.1 2010/01/14:
  run-id 2019060291  delaunay  Qtriangulate  Qbbound-last  Qcoplanar-keep
  _pre-merge  _zero-centrum  Qinterior-keep  Pgood  _max-width  1
  Error-roundoff 1.4e-15  _one-merge 9.7e-15  Visible-distance 2.8e-15
  U-coplanar-distance 2.8e-15  Width-outside 5.5e-15  _wide-facet 1.7e-14

The input to qhull appears to be less than 3 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p3(v3):     1     0   0.5
- p1(v2):     0     1   0.5
- p2(v1):     1     1     1
- p0(v0):     0     0     0

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 1.4e-15.  The center point, facets and distances
to the center point are as follows:

center point      0.5      0.5      0.5

facet p1 p2 p0 distance=    0
facet p3 p2 p0 distance=    0
facet p3 p1 p0 distance=    0
facet p3 p1 p2 distance=    0

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:         0         1  difference=    1
  1:         0         1  difference=    1
  2:         0         1  difference=    1

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 1.4e-15.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.


This is a Delaunay triangulation and the input is co-circular or co-spherical:
  - use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)
  - or use 'QJ' to joggle the input and avoid co-circular data


Matlab gives,


x = [0 0 1 1];
y = [0 1 1 0];
tri = delaunay (x, y)

tri =

     3     2     4
     2     1     4


After this error, the dellunay test also fails.


test delaunay
QH6164 qhull internal error (qh_qhull): temporary sets not empty(2)

While executing:  | qhull d Qt Qbb Qc
Options selected for Qhull 2010.1 2010/01/14:
  run-id 2032169751  delaunay  Qtriangulate  Qbbound-last  Qcoplanar-keep
  _pre-merge  _zero-centrum  Qinterior-keep  Pgood  _max-width  2
  Error-roundoff 2.8e-15  _one-merge 1.9e-14  Visible-distance 5.5e-15
  U-coplanar-distance 5.5e-15  Width-outside 1.1e-14  _wide-facet 3.3e-14
Last point added to hull was p3.  Last merge was #1.

At error exit:

Delaunay triangulation by the convex hull of 5 points in 3-d:

  Number of input sites: 5
  Number of Delaunay regions: 0
  Number of non-simplicial Delaunay regions: 1

Statistics for:  | qhull d Qt Qbb Qc

  Number of points processed: 5
  Number of hyperplanes created: 7
  Number of facets in hull: 5
  Number of distance tests for qhull: 6
  Number of distance tests for merging: 29
  Number of distance tests for checking: 30
  Number of merged facets: 1

  ***** testif HAVE_QHULL
 x = [-1, 0, 1, 0, 0];
 y = [0, 1, 0, -1, 0];
 assert (sortrows (sort (delaunay (x, y), 2)), [1,2,5;1,4,5;2,3,5;3,4,5])
!!!!! test failed
__delaunayn__: qhull failed


If I run "clear all" then the test passes.

This may be a problem local to building on MacOS, with MacPorts, or something I am doing. I'd appreciate the effort if someone else would try the example.

Ben Abbott <bpabbott>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by bpabbott (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2011-10-25 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2011-10-21 rik5 StatusConfirmed In Progress
    2011-10-20 bpabbott Operating SystemMac OS Any
    2011-10-20 rik5 StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code