[Octave-bug-tracker] [bug #39370] contourc adds duplicates to the start and end of the 1st level

[Rik]

Follow-up Comment #14, bug #39370 (project octave):

I think this may be simpler than thought, although it probably does require a
flag.

Going back to this example


>> x = 0:2;
>> y = x;
>> z = x' * y
z =
     0     0     0
     0     1     2
     0     2     4
>> b = contour (x, y, z, 0:4)
b =
  Columns 1 through 14
    1.0000    0.5000    1.0000    1.0000    2.0000    2.0000    1.0000   
1.0000    2.0000    2.0000    3.0000    1.5000    2.0000    4.0000
    4.0000    2.0000    1.0000    1.0000    0.5000    4.0000    2.0000   
2.0000    1.0000    1.0000    2.0000    2.0000    1.5000    2.0000
  Columns 15 through 16
    2.0000    2.0000
    2.0000    2.0000
>> c = contourc (x, y, z, 0:4)
c =
  Columns 1 through 14
         0         0         0    1.0000    2.0000    1.0000    0.5000   
1.0000    2.0000    2.0000    1.0000    2.0000    3.0000    1.5000
    4.0000    2.0000    1.0000         0         0    3.0000    2.0000   
1.0000    0.5000    2.0000    2.0000    1.0000    2.0000    2.0000
  Columns 15 through 17
    2.0000    4.0000    2.0000
    1.5000    1.0000    2.0000


There are two differences: 1) Matlab contourc includes the 0 level, and 2)
Matlab contourc eliminates duplicate points.

I think the first difference is because of the use of '<' versus '<=' to
include or exclude points in Octave.  The request for a level of 0 when Z just
equals zero seems to be different.  A way to see this is to change the level
in Octave from 0 to eps.


c = contourc (x,y,z, [eps, 1:4])
c =

 Columns 1 through 7:

   2.2204e-16   2.2204e-16   2.2204e-16   1.0000e+00   2.0000e+00   1.0000e+00
  5.0000e-01
   4.0000e+00   2.0000e+00   1.0000e+00   2.2204e-16   1.1102e-16   4.0000e+00
  2.0000e+00

 Columns 8 through 14:

   1.0000e+00   1.0000e+00   2.0000e+00   2.0000e+00   1.0000e+00   1.0000e+00
  2.0000e+00
   1.0000e+00   1.0000e+00   5.0000e-01   4.0000e+00   2.0000e+00   2.0000e+00
  1.0000e+00

 Columns 15 through 21:

   2.0000e+00   3.0000e+00   1.5000e+00   2.0000e+00   4.0000e+00   2.0000e+00
  2.0000e+00
   1.0000e+00   2.0000e+00   2.0000e+00   1.5000e+00   2.0000e+00   2.0000e+00
  2.0000e+00


If you de-dupe this then you get the Matlab result.

This makes me thing that the solution is going to be to go into the C++ code. 
It looks like the mark_facets() function in __contour__.cc might need
adjusting.  A code sample is


  for (unsigned int c = 0; c < nc; c++)
    for (unsigned int r = 0; r < nr; r++)
      {
        f[0] = Z(r, c) - lvl;
        f[1] = Z(r, c+1) - lvl;
        f[3] = Z(r+1, c) - lvl;
        f[2] = Z(r+1, c+1) - lvl;

        for (unsigned int i = 0; i < 4; i++)
          if (fabs(f[i]) < std::numeric_limits<double>::epsilon ())
            f[i] = std::numeric_limits<double>::epsilon ();

        if (f[1] * f[2] < 0)
          mark(r, c) += 2;

        if (f[0] * f[3] < 0)
          mark(r, c) += 8;
      }


Octave is trying to account for points that are numerically close to the
boundary, but I think it may be getting it wrong.  Recall that for this
example,


z =
     0     0     0
     0     1     2
     0     2     4


So for point (1,1) f[0] = 0 and f[1] = 0 and they are both replaced by eps. 
The test is then 


        if (f[1] * f[2] < 0)
          mark(r, c) += 2;


which will not be true because the multiplication results in eps ^ 2.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?39370>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/