bugGNU Octave - Bugs: bug #39370, contourc adds duplicates to the...

 
 

bug #39370: contourc adds duplicates to the start and end of the 1st level

Submitter:  Ben Abbott <bpabbott>
Submitted:  Sun 30 Jun 2013 09:30:19 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Confirmed Assigned to:  None
Originator Name:  Ben Abbott Open/Closed:  * Open
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 26 Jul 2022 04:50:17 PM UTC, comment #15: 

interesting note on eps.  I saw that part of mark_facets. figured it was just looking for sign-flips, and for that purpose it considered anything 0 +/- eps to be "0" but pushed it to eps to force a consistent sign.

Nicholas Jankowski <nrjank>
Group Member
Tue 26 Jul 2022 04:31:17 PM UTC, comment #14: 

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.

Rik <rik5>
Group administrator
Tue 26 Jul 2022 02:57:34 PM UTC, comment #13: 

looking again at matlab's old documentation, they had a 'contour algorithm' explanation


"Keep checking the edges of each square entered to determine the exit edge until the line(cx,cy) closes on its initial point or exits the grid."

https://lost-contact.mit.edu/afs/inf.ed.ac.uk/group/teaching/matlab-help/R2013b/matlab/creating_plots/the-contouring-algorithm.html

Wouldn't surprise me if there was an explicit 'if within eps, make the points exact' step in their algorithm. that would likely help merge intersecting/branching contours, but now that I try to think about it I'm not sure if continuous functions can have branches/intersections.

Nicholas Jankowski <nrjank>
Group Member
Tue 26 Jul 2022 02:34:25 PM UTC, comment #12: 

thanks for the quick followup.  I realize that there's no difference when plotting, but I have seen some use cases where people use contourc to generate a level contour for computation not plotting, and I assume the idea with the revised contourc is to clean up the redundancy for things like that. may also explain why it keeps things like the lev=0 case in the original set that gets thrown out for plotting.

checking in matlab for something similar to your sombrero plot (using peaks instead), it looks like their algorithm in both contour and contourc keep both endpoints for closed contours, and they are 'exact'. so, yes, octave applying a blind 'unique' trim postprocess may be a problem there as you'd lose the closing point to the contour. In general keeping them with O(eps) differences are likely irrelevant for plotting, but i guess it could create a leaky contour if it was used to define a boundary in a numerical analysis that assumed it didn't need to correct for it. can the contour algo produce intersecting/branching contours with similar issues or are they necessarily single paths?


>> [xx,yy,zz] = peaks();
>> a = contour(xx,yy,zz,[6 7 8])
a =
  Columns 1 through 14
    6.0000   -0.2668   -0.2500   -0.1250         0    0.1250    0.2244    0.2500    0.3750    0.4600    0.5000    0.5417    0.5488    0.5000
   31.0000    1.2500    1.2433    1.2169    1.2127    1.2266    1.2500    1.2580    1.3176    1.3750    1.4296    1.5000    1.6250    1.7446
  Columns 15 through 28
    0.4974    0.3750    0.3695    0.2500    0.1250         0   -0.1250   -0.2500   -0.3714   -0.3750   -0.4930   -0.5000   -0.5382   -0.5276
    1.7500    1.8711    1.8750    1.9349    1.9708    1.9829    1.9719    1.9366    1.8750    1.8725    1.7500    1.7337    1.6250    1.5000
  Columns 29 through 42
   -0.5000   -0.4536   -0.3750   -0.2668    7.0000   -0.2515   -0.2500   -0.1250         0    0.1250    0.2202    0.2500    0.3591    0.3750
    1.4476    1.3750    1.3095    1.2500   23.0000    1.3750    1.3741    1.3335    1.3246    1.3425    1.3750    1.3934    1.5000    1.6176
  Columns 43 through 56
    0.3758    0.3750    0.2953    0.2500    0.1250         0   -0.1250   -0.2500   -0.3000   -0.3750   -0.3789   -0.3750   -0.3664   -0.2515
    1.6250    1.6266    1.7500    1.7831    1.8383    1.8570    1.8405    1.7870    1.7500    1.6333    1.6250    1.5788    1.5000    1.3750
  Columns 57 through 62
    8.0000   -0.0849         0    0.0670         0   -0.0849
    5.0000    1.6250    1.5054    1.6250    1.6483    1.6250
>> b = contourc(xx(1,:),yy(:,1),zz,[6 7 8])
b =
  Columns 1 through 14
    6.0000   -0.5000   -0.4930   -0.3750   -0.3714   -0.2500   -0.1250         0    0.1250    0.2500    0.3695    0.3750    0.4974    0.5000
   31.0000    1.7337    1.7500    1.8725    1.8750    1.9366    1.9719    1.9829    1.9708    1.9349    1.8750    1.8711    1.7500    1.7446
  Columns 15 through 28
    0.5488    0.5417    0.5000    0.4600    0.3750    0.2500    0.2244    0.1250         0   -0.1250   -0.2500   -0.2668   -0.3750   -0.4536
    1.6250    1.5000    1.4296    1.3750    1.3176    1.2580    1.2500    1.2266    1.2127    1.2169    1.2433    1.2500    1.3095    1.3750
  Columns 29 through 42
   -0.5000   -0.5276   -0.5382   -0.5000    7.0000   -0.3750   -0.3000   -0.2500   -0.1250         0    0.1250    0.2500    0.2953    0.3750
    1.4476    1.5000    1.6250    1.7337   23.0000    1.6333    1.7500    1.7870    1.8405    1.8570    1.8383    1.7831    1.7500    1.6266
  Columns 43 through 56
    0.3758    0.3750    0.3591    0.2500    0.2202    0.1250         0   -0.1250   -0.2500   -0.2515   -0.3664   -0.3750   -0.3789   -0.3750
    1.6250    1.6176    1.5000    1.3934    1.3750    1.3425    1.3246    1.3335    1.3741    1.3750    1.5000    1.5788    1.6250    1.6333
  Columns 57 through 62
    8.0000         0    0.0670         0   -0.0849         0
    5.0000    1.6483    1.6250    1.5054    1.6250    1.6483

>> a(:,2)-a(:,32)
ans =
     0
     0
>> a(:,56)-a(:,34)
ans =
     0
     0
>> a(:,58)-a(:,62)
ans =
     0
     0
>> b(:,2)-b(:,32)
ans =
     0
     0
>> b(:,56)-b(:,34)
ans =
     0
     0
>> b(:,58)-b(:,62)
ans =
     0
     0


Nicholas Jankowski <nrjank>
Group Member
Tue 26 Jul 2022 01:39:43 AM UTC, comment #11: 

Just for fun, I put tic() and toc() around the uniquifying code and then did


[xx,yy,zz] = sombrero ();
c = contourc (xx,yy,zz, 10);


Running this a few times the extra code adds 10-13 milliseconds to the run time which isn't much.  I then increased the number of levels to 50 to stress the code, but it still took only 25 milliseconds.

I don't think performance is the issue, but I do think finite precision is.  For example,


[xx,yy,zz] = sombrero ();
c = contourc (xx,yy,zz, [0.95, 0.98, 0.99])
c =

 Columns 1 through 7:

    0.9500   -0.3618         0    0.3618    0.4000    0.5226    0.4000
   13.0000   -0.4000   -0.5226   -0.4000   -0.3618         0    0.3618

 Columns 8 through 14:

    0.3618         0   -0.3618   -0.4000   -0.5226   -0.4000   -0.3618
    0.4000    0.5226    0.4000    0.3618         0   -0.3618   -0.4000

 Columns 15 through 21:

    0.9800   -0.3024         0    0.3024         0   -0.3024    0.9900
    5.0000         0   -0.3024         0    0.3024         0    5.0000

 Columns 22 through 26:

   -0.1512         0    0.1512         0   -0.1512
         0   -0.1512         0    0.1512         0


For the contour at 0.95, the first and last entries are within an eps of each other.


c(:,2) - c(:,14)
ans =

  -5.5511e-17
            0


Perhaps this speaks to working on this in the C++ code.

Rik <rik5>
Group administrator
Tue 26 Jul 2022 01:28:03 AM UTC, comment #10: 

Good debugging.

I asked a question in comment #4 back in 2014 about whether there was any actual effect from having duplicate points in the contour lines.  I don't think there is which means this may not be that important to fix this.

Quoting the example in comment #8:


>> x = 0:2;
>> y = x;
>> z = x' * y;
>> c = contourc (x, y, z, 2:3)

c =
    2.0000    1.0000    2.0000    3.0000    1.5000    2.0000
    2.0000    2.0000    1.0000    2.0000    2.0000    1.5000

>> b = contour (x, y, z, 2:3)

b =
  Columns 1 through 7
    2.0000    1.0000    1.0000    2.0000    2.0000    3.0000    1.5000
    4.0000    2.0000    2.0000    1.0000    1.0000    2.0000    2.0000
  Column 8
    2.0000
    1.5000


You have to understand the contourc format, but the first contour is the points


(1, 2) (1, 2) (2, 1) (2, 1)


These points, in Octave, are used to directly create a patch object so the duplication is unnecessary, but immaterial.

I also think this means we don't need to follow Matlab exactly.  If the uniquified matrix from contourc will work to produce a plot then we might as well just settle on both routines producing the same matrix which is the smaller one.

At least for starters, I would probably do the culling in contourc.m.  It is an m-file, and potentially slow, but the number of contours is not that large (default is 10).

This would be approximately right Octave code


diff -r 7d4cf04665e6 scripts/plot/draw/contourc.m
--- a/scripts/plot/draw/contourc.m        Wed Jul 20 16:37:58 2022 +0200
+++ b/scripts/plot/draw/contourc.m        Mon Jul 25 18:25:21 2022 -0700
@@ -163,6 +163,28 @@ function [c, lev] = contourc (varargin)
     endwhile
   endif

+  ## Post-process to remove duplicate contour points
+  ## FIXME: Possibly better to do this in C++ code __contour__
+  i = j = 0;
+  c2 = zeros (size (c));
+  while (i < columns (c))
+    i++;  j++;
+
+    lvl = c(1,i);
+    n = c(2,i);
+
+    ctmp = c(:, i+(1:n));
+    ctmp = unique (ctmp.', "rows", "stable").';
+    n2 = columns (ctmp);
+
+    c2(1,j) = lvl;
+    c2(2,j) = n2;
+    c2(:,j+(1:n2)) = ctmp;
+
+    i += n;  j += n2;
+  endwhile
+  c = c2(:,1:j);
+
 endfunction


Note that this won't catch points that differ by a few eps.  One would need to use uniquetol for that, but it doesn't appear to have a "stable" option.


(file #53472)

Rik <rik5>
Group administrator
Mon 25 Jul 2022 08:25:08 PM UTC, comment #9: 

ok, according to the discussion linked below, the change in contourc occurred sometime around matlab 7.4 (2009a).

https://www.mathworks.com/matlabcentral/answers/95263-why-do-contour-and-contourc-produce-different-results-in-matlab-7-8-r2009a

This seems to confirm the problem isn't that contourc adds duplicates, as that's expected by contour, it's that octave's contourc doesn't trim the dupes from the matrix.

a compatibility hack/fix for contourc may just involve applying a unique check and trim on the output. unfortunately contour still expects the full contourc output unmodified. perhaps adding a 'legacy' flag to contourc that defaults to false with the new trimmed output. that would enable continued use of the code by contour and other functions without breaking them.

Regarding any other issues with the two functions - I also notice  that the outputs seems to have the same number of points for larger cases, the ordering is sometimes but not always different (may have the same points but one sorts with increasing in y while the other sorts decreasing), and contour seems to throw out some levels that contourc keeps.

slightly expanding on the simple 3pt example with other levels:

matlab 2022a:

>> 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


note contourc keeps the 0 level while contour doesn't, and contourc trims the dupes from each level while contour leaves them in.  octave currently matches the output of contour, with no 0 level and dupes left in.  (so maybe the 'legacy' flag would also have to somehow prevent clipping some data levels. haven't finished stepping through to see where the 0 level currently gets trimmed out.

Nicholas Jankowski <nrjank>
Group Member
Mon 25 Jul 2022 02:11:15 PM UTC, comment #8: 

looking at this a bit, it seems that octave is just pushing out the same output that comes from contour, while matlab truncates the additional points:

Octave 8.0.0:


>> x = 0:2;
>> y = x;
>> z = x' * y;
>> a = contourc (x, y, z, 2:3)
a =

   2.0000   1.0000   1.0000   2.0000   2.0000   3.0000   1.5000   2.0000
   4.0000   2.0000   2.0000   1.0000   1.0000   2.0000   2.0000   1.5000

>> b = contour (x, y, z, 2:3)
b =

   2.0000   1.0000   1.0000   2.0000   2.0000   3.0000   1.5000   2.0000
   4.0000   2.0000   2.0000   1.0000   1.0000   2.0000   2.0000   1.5000


whereas with matlab:


>> x = 0:2;
y = x;
z = x' * y;
c = contourc (x, y, z, 2:3)
c =
    2.0000    1.0000    2.0000    3.0000    1.5000    2.0000
    2.0000    2.0000    1.0000    2.0000    2.0000    1.5000
>> b = contour (x, y, z, 2:3)
b =
  Columns 1 through 7
    2.0000    1.0000    1.0000    2.0000    2.0000    3.0000    1.5000
    4.0000    2.0000    2.0000    1.0000    1.0000    2.0000    2.0000
  Column 8
    2.0000
    1.5000


so i don't think the existance of the duplication is necessarily a bug.  the fact that the duplication isn't pared down in contourc might be.

currently the matlab help for contourc states:

The matrix returned by contourc might be inconsistent with the results from the contour, contourf, and contour3 functions. To get the matrix used by one of those functions, call that function with an output argument.
https://www.mathworks.com/help/matlab/ref/contourc.html


whereas if you find some old documentation you'll find statements like:

Description
contourc calculates the contour matrix C used by the M-file contour to draw actual contour plots.
...
http://www.math.clemson.edu/~warner/M860/Matlab/contour.html


None of that explains WHY contour in both programs, and originally both contourc's produce the duplicate points, but at least it seems it was historically consistent. at some point between the 1994 docs linked above and this 2013 report, mathworks decided contourc should produce a non-duplicative set. if it's simple and consistent, it may not take much to have the contourc m-file just clip out those values.

Nicholas Jankowski <nrjank>
Group Member
Fri 22 Jul 2022 09:01:48 PM UTC, comment #7: 

bug still present in dev branch of Octave 8.0.0

Nicholas Jankowski <nrjank>
Group Member
Tue 08 Dec 2020 03:39:19 AM UTC, comment #6: 

as of Octave v.6.1.0, the examples from comment #0:

The help contourc help example now reads:


x = 0:2;
y = x;
z = x' * y;
c = contourc (x, y, z, 2:3)
 => c =
       2.0000   1.0000   1.0000   2.0000   2.0000   3.0000   1.5000   2.0000
       4.0000   2.0000   2.0000   1.0000   1.0000   2.0000   2.0000   1.5000


which matches wat Octave 6.1.0 produces:


>> x = 0:2;
>> y = x;
>> z = x' * y;
>> a = contourc (x, y, z, 2:3)
a =

   2.0000   1.0000   1.0000   2.0000   2.0000   3.0000   1.5000   2.0000
   4.0000   2.0000   2.0000   1.0000   1.0000   2.0000   2.0000   1.5000


Matlab 2020b produces:

>> x = 0:2;
>> y = x;
>> z = x' * y;
>> contourc (x, y, z, 2:3)
ans =
    2.0000    1.0000    2.0000    3.0000    1.5000    2.0000
    2.0000    2.0000    1.0000    2.0000    2.0000    1.5000


Nicholas Jankowski <nrjank>
Group Member
Mon 21 Nov 2016 10:09:13 PM UTC, comment #5: 

This issue is still present in Octave 4.2.0.

Hartmut <hardy>
Wed 31 Dec 2014 06:02:02 AM UTC, comment #4: 

At first this looks bad, but I think this is mostly a performance hit.  The duplicate points shouldn't affect the drawing of the contour, although I can't guarantee that.

I took a look at _contour_.cc and the code is very old and I couldn't figure out what was going on easily.

Rik <rik5>
Group administrator
Sat 06 Jul 2013 03:45:58 PM UTC, comment #3: 

Ok. It is a bug in _contourc_.cc

The first level list has duplicates for the beginning and end.

x = [1 1 2 2];
y = [2 2 1 1];

is the same as;

x = [1 2];
y = [2 1];


Ben Abbott <bpabbott>
Group Member
Sat 06 Jul 2013 03:08:53 PM UTC, comment #2: 

To be honest, I'm not sure if this is a bug or a documentation error.  I'll try to find the time to dig deeper in the next few days.

Ben Abbott <bpabbott>
Group Member
Wed 03 Jul 2013 03:56:05 PM UTC, comment #1: 

Is the documentation on the example incorrect or is the contourc code incorrect?  It seems like Octave has been generating a 2x8 matrix for the example since at least version 3.2.4.  Meanwhile, Matlab and the documentation suggest a 2x6 matrix is correct.

Maybe we should re-title the report.

Rik <rik5>
Group administrator
Sun 30 Jun 2013 09:30:19 AM UTC, original submission:  

"help contourc" displays the example below.


          x = 0:2;
          y = x;
          z = x' * y;
          contourc (x, y, z, 2:3)
             =>   2.0000   2.0000   1.0000   3.0000   1.5000   2.0000
                  2.0000   1.0000   2.0000   2.0000   2.0000   1.5000


The sources produce this result.


a = contourc (x, y, z, 2:3)
a =

   2.0000   1.0000   1.0000   2.0000   2.0000   3.0000   1.5000   2.0000
   4.0000   2.0000   2.0000   1.0000   1.0000   2.0000   2.0000   1.5000


Matlab doesn't match either result.


ans =

    2.0000    1.0000    2.0000    3.0000    1.5000    2.0000
    2.0000    2.0000    1.0000    2.0000    2.0000    1.5000


Ben Abbott <bpabbott>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #53472:  contourc.diff added by rik5 (754B - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -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
    2022-07-26 rik5 Attached File- Added contourc.diff, #53472
    2014-12-31 rik5 StatusNone Confirmed
    2013-10-29 mtmiller CategoryDocumentation Octave Function
        Item GroupDocumentation Matlab Compatibility
    2013-07-06 bpabbott Summarycontourc example is incorrect contourc adds duplicates to the start and end of the 1st level

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code