bugGNU Octave - Bugs: bug #46946, isosurface throws errors where...

 
 

bug #46946: isosurface throws errors where Matlab does not

Submitter:  Markus Mützel <mmuetzel>
Submitted:  Thu 21 Jan 2016 01:33:34 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  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 05 Jul 2016 12:02:16 PM UTC, comment #21: 

Yeah I saw that too, but too late :-)

I was misled by the fact that after a few isosurface manipulations the test did work.
But now I realize that happens because then __calc_isovalue_from_data__ has been called from isosurface and is already in the symbol table.
Normally, funcs in private/ subdirs can't be called directly (they aren't in the path), but once called from the parent subdir and parsed they do have a symbol table slot. Once the symbol table has a specific entry, that entry doesn't need to be searched on disk anymore. I have no idea about scope etc. but that must be the way it works.

I'd rather leave the test in. Someday, somehow we'll have a fix for testing private and nested funcs and subfuncs.

BTW I like the fact that isosurface has been improved quite a bit. The results look very nice. Thank you very much for that!

Philip Nienhuis <philipnienhuis>
Group Member
Tue 05 Jul 2016 10:18:13 AM UTC, comment #20: 

Philip, thank you for pushing the patch.

I misunderstood the impact of the test of a private function. I ran "test isosurface" from within Octave probably after using the function which did pass for me. I did not try "make check".

I see that Kai already converted the test to an xtest. I would also be OK with removing that test completely.

Markus Mützel <mmuetzel>
Group administrator
Mon 04 Jul 2016 08:35:59 PM UTC, comment #19: 

Thanks, I pushed your latest cset here:
http://hg.savannah.gnu.org/hgweb/octave/rev/634fbedbfb5b

with one mod: Mützel => Muetzel  I hope you don't mind :-)
My local copy of mercurial doesn't like umlauts :-(

I do get a FAIL & test isosurface fails, both the same and due to _calc_isovalue_from_data_ not being found. However:

(1) it is present in plot/draw/private/

(2) after running the examples it suddenly does work (PASSes 31 of 31 tests)

... so I guess it must be my machine's fault, or my cross-built dev version's fault.



Philip Nienhuis <philipnienhuis>
Group Member
Sun 03 Jul 2016 09:50:10 PM UTC, comment #18: 

Thank you, Philip, for your feedback.

In the attached patch, the test in _calc_isovalue_from_data_ is moved to isosurface.m.

The first example in the docstring give a isosurface for random data. So it is quite possible that there are some holes between the patches. Whether this example is very instructive is debatable. If you prefer, I could change that example to a sphere or anything else.

As light is now available in core Octave, I un-commented the respective lines of the second example in the docstring.

(file #37696)

Markus Mützel <mmuetzel>
Group administrator
Sun 03 Jul 2016 11:16:40 AM UTC, comment #17: 

#6 seems to work fine; the first example gives a nice 3D plot with some holes between patches (intended ?).

"test isosurface" PASSes 30 of 30 tests.
However, "make check" gives a FAIL for

:
  scripts/plot/draw/private/__calc_isovalue_from_data__.m ..... PASS      0/1
                                                                  FAIL    1


...and that is because it is in a private subdir which isn't in the path. Once in Octave, cd'ing to plot/draw/private:

>> cd (strrep (which ("isosurface"), "isosurface.m", "private"))
>> test __calc_isovalue_from_data__
PASSES 1 out of 1 test


I suppose it is better to move that test into isosurface.m itself.
FYI there's a bug report about tests in private subdirs: bug #38776. that may be of interest to you as regards


__calc_isovalue__




Philip Nienhuis <philipnienhuis>
Group Member
Sat 02 Jul 2016 08:30:25 PM UTC, comment #16: 

Well I was just ready to make another cross-build (the 5th and last today I hope) so I'll back out the previous and apply the latest patch. Thanks for your prompt action!

There's no _unite_shared_vertices_.m in .../plot/draw/private/

btw, Release -> dev, Status -> Ready for test

Philip Nienhuis <philipnienhuis>
Group Member
Sat 02 Jul 2016 08:23:03 PM UTC, comment #15: 

Philip, thank you for trying to test this patch.
Is this file present in scripts/plot/draw/private? The function from the error message should be one of the two files added by this patch.

These two files had Windows line endings. The attached patch fixes this. Additionally, I added these two files to the corresponding make-files.
Could you please try with the attached patch "isosurface_6.patch"?

(file #37674)

Markus Mützel <mmuetzel>
Group administrator
Sat 02 Jul 2016 05:54:52 PM UTC, comment #14: 

On Win7 I tried the first example in the texinfo header and got this:

>> [x, y, z] = meshgrid (1:5, 1:5, 1:5);
>>           val = rand (5, 5, 5);
>>           isosurface (x, y, z, val, .5)
error: '__unite_shared_vertices__' undefined near line 159 column 33
error: called from
    isosurface at line 159 column 31
>>


(Using file #36998)

Philip Nienhuis <philipnienhuis>
Group Member
Sat 02 Jul 2016 12:01:03 PM UTC, comment #13: 

Lachlan, sorry for the confusion I might have caused. The last comment was ment to approve that Marco was right with his comments on an earlier version of this patch.

The latest version file #36998 takes a different approach and should not suffer from that error.

Markus Mützel <mmuetzel>
Group administrator
Sat 02 Jul 2016 01:51:36 AM UTC, comment #12: 

Markus, did you intend to attache a patch to comment #11?  It seems that Marco pointed out a problem with the latest version of this patch that I can see.

Retagging as "In progress".

Lachlan Andrew <lachlan>
Wed 27 Apr 2016 01:00:41 PM UTC, comment #11: 

You are right. I did not consider that in the regime where things get interesting here the difference between the two points is in the order of eps of the values...

dist = (a-b)^2
error = sqrt( (2*(a-b) * eps(a))^2 + (2*(a-b) * eps(b))^2 )

Assuming close points where eps(a) = eps(b) and with (a-b) = eps(a), this reduces to

error = 2 * eps(a)^2


I hope the new version that compares differences instead of distances is correct now...

Markus Mützel <mmuetzel>
Group administrator
Wed 27 Apr 2016 12:06:05 PM UTC, comment #10: 

Dear Markus,

I will try your new patch asap. About eps/eps^2, I just tested


fvc.vertices=[2/3,2/3,2/3;1-1/3,1-1/3,1-1/3];
dx = fvc.vertices(:,1) .- fvc.vertices(:,1)';
dy = fvc.vertices(:,2) .- fvc.vertices(:,2)';
dz = fvc.vertices(:,3) .- fvc.vertices(:,3)';
mutual_dist = dx.^2 + dy.^2 + dz.^2


and observed an eps^2. But I see that in your new patch you just consider differences along each coordinate and then eps is fine.

Marco Caliari <caliari>
Group Member
Tue 26 Apr 2016 05:14:24 PM UTC, comment #9: 

Dear Marco,

thank you for having a look at this patch. I am happy it could help you.

You raised some valid points. It is not correct to use the same eps for all distances. Like you have pointed out, it should be different for each two points. I did not consider sets of vertices that span over several orders of magnitudes.

I like your approach of calculating the mutual distances between all points. It should save memory and be faster than the trivial approach I took.

However, I think it is not correct to compare to eps^2. Usually, the error always increases when a value is squared. Probably, it is better to assume that the uncertainties (floating point precission) are independent in each coordinate and use the laws of error propagation (e.g. here [1]). Simplifying by the assumption that the uncertainties are similar in each coordinate, that would read (employing your approach for calculating the mutual distances):

IP = fvc.vertices * fvc.vertices';
norms = diag (IP);
mutual_dist = norms + norms' - 2*IP;
is_very_close = mutual_dist < 2 * sqrt (3) * eps (max (norms, norms'));


I had forgotten the squareroot in the patch.
I hope this is right. Please, correct me if I am wrong.

But that call to eps takes enormous time on my machine ("demo surface 1" takes 22 seconds on my system). Making the simplification "eps (x) = x * eps" performs much faster (0.02 seconds):

IP = fvc.vertices * fvc.vertices';
norms = diag (IP);
mutual_dist = norms + norms' - 2*IP;
is_very_close = mutual_dist < 2 * sqrt (3) * eps * max (norms, norms');


But this still uses a lot of memory. The following example does not run for me (out of memory):

[x,y,z] = meshgrid (-2:0.05:2, -2:0.05:2, -2:0.5:2);
val = x.^2 + y.^2 + z.^2;
isosurface2 (x, y, z, val, 1);


I ran into a similar issue for reducepatch (patch #8912). Using _unite_shared_vertices_ from that patch, the above example runs in below 1.2s for me (not very fast but maybe acceptable). "demo surface 1" takes 0.03 seconds (still ok imho).

Attached please find an updated patch with those changes. I incorporated your change in lines 51&52 of _unite_shared_vertices_.m.

[1] https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Simplification


(file #36998)

Markus Mützel <mmuetzel>
Group administrator
Tue 26 Apr 2016 11:21:48 AM UTC, comment #8: 

Dear Markus,

I recently encountered your original bug 3., thanks for fixing it.
I have some comments about your latest patch. First of all, since you compute the square of the distances, I think you should compare with eps^2. Moreover, max (abs (fvc.vertices(:))) is not correct: try the following set


fvc.vertices = [1,1,1;1e20,1e20,1e20;1.1,1.1,1.1;10,10,10]


You will get first, third and fourth point "very_close". The distance between two points should be relative to the maximum 2-norm of those two points. I tried the following code


x = [fvc.vertices(:,1)';fvc.vertices(:,2)';fvc.vertices(:,3)'];
IP = x' * x;
norms = diag(IP);
mutual_dist = norms+norms'-2*IP;
is_very_close = mutual_dist < SPF * eps^2 * max(norms',norms);


where SPF is some possible factor and it seems to work nice. What do you think?

Marco Caliari <caliari>
Group Member
Wed 17 Feb 2016 08:53:41 AM UTC, comment #7: 

No problem if it takes a little longer. Thus, the patch can "mature" a little:
I realized that "unique" does not really do the job for recognizing shared vertices. This is because very close doubles are treated as different values:

a = [2/3 0 0; 2/3 1 0; 1-1/3 0 0];
unique(a, "rows")

Result:

ans =
   0.66667   0.00000   0.00000
   0.66667   1.00000   0.00000
   0.66667   0.00000   0.00000

This is because

>> a(1,:)-a(3,:)
ans =
  -1.1102e-016  0.0000e+000  0.0000e+000

>> eps(2/3)
ans =   1.1102e-016


At first I thought this was a bug. But Matlab behaves exactly the same. And it is true that 2/3 != 1-1/3 in double precision.

Therefore, I changed the logic for detection of shared vertices by calculating the mutual distances between all vertices and comparing to "eps" (isosurface_4.patch). There is probably a much better way than my current implementation since it performs e.g. more than twice the calculations really needed. But I cannot think of one without using loops which would perform even worse in Octave.
So feel free to give advise.


(file #36360, file #36361)

Markus Mützel <mmuetzel>
Group administrator
Mon 08 Feb 2016 09:04:18 PM UTC, comment #6: 

Don't worry, your contributions will make it into Octave - some day :-) I'm sure they'll make it into 4.2.x (the next major release) - but that's still far away (end of this year?).
The main developers are quite busy with revamping the source tree and getting 4.0.1 released. A lot of work.

I've asked (in the patch tracker for reducevolume.m) where your contibutions should be put (in what subdir of/scripts in the dev tree, or -after installetion- what subdir of /usr/local/share/octave/<version>/m/ ). Once I get an answer I'll proceed further, but mind you, I have little spare time as well.

Philip Nienhuis <philipnienhuis>
Group Member
Mon 08 Feb 2016 06:04:30 PM UTC, comment #5: 

Sorry for posting yet another patch.

Matlab returns vertices shared by adjacent faces only once. At the moment isosurface returns some of the vertices several times. This corresponds to the behavior of Matlab when the argument "noshare" is passed to the function.
I changed the patch such that "unique" is called to find shared vertices unless you pass the additional argument "noshare".
The documentation and the tests are changed respectively.

Is there any chance that this patch will make it into the distribution some day?

(file #36278, file #36279, file #36280)

Markus Mützel <mmuetzel>
Group administrator
Tue 26 Jan 2016 10:25:30 AM UTC, comment #4: 

Thank you for spotting that. I copied&pasted most of the input argument checking from isocaps and must have overlooked those references.

Attached please find isosurface_2.diff with your comments addressed.

(file #36168)

Markus Mützel <mmuetzel>
Group administrator
Mon 25 Jan 2016 10:20:26 PM UTC, comment #3: 

A quick glance:
In L.198 I read in the error message "isocaps.m: ..." where this function is a subfunc of isosurface.m.

In l.518: I read "called from isocaps" but as you say in comment #2, also from isosurface ?

Philip Nienhuis <philipnienhuis>
Group Member
Mon 25 Jan 2016 09:35:15 AM UTC, comment #2: 

Changeset for this bug attached.

I adressed all points mentioned and adapted the help text and the tests.

Points 1.-3. indeed are only changes to the input argument checking.
Point 4. is added functionality. I added '__calc_isovalue_from_data__.m' for this (is the name matching the conventions?). The reason why it is in a separate file is that I would like to re-use it for patch #8872 ('isocaps.m').

While working with the function, I found another difference to Matlab: When calling isosurface without output arguments with different iso values repeatedly in Matlab, it adds a new isosurface to the same axes each time. In Octave, a new axes is created each time.
I removed the command

hax = newplot ();


to match the behaviour of Matlab.

Please, let me know if something needs to be changed.

(file #36160)

Markus Mützel <mmuetzel>
Group administrator
Sat 23 Jan 2016 09:57:10 PM UTC, comment #1: 

A changeset is always welcome.  A lot of this should be easy to fix since it is in an m-file rather than C++ code.  The problem appears to be very limited input argument checking in isosurface.m, before it goes ahead and calls _marching_cube_ and patch.

Rik <rik5>
Group administrator
Thu 21 Jan 2016 01:33:34 PM UTC, original submission:  

1. 'verbose'

In Matlab, when you run this line

isosurface (rand (10, 10, 10), .5, 'verbose')

a figure with the isosurface is shown together with the following lines in the command window:

ISOSURFACE: Computing triangles and vertices...
         V-done
.........
ISOSURFACE: number of vertices=1335   number of triangles=2287


When you run the same lines in Octave, the following error is thrown:

error: COLORS must be a matrix of dim 3 and of same size as C
error: called from
    __marching_cube__ at line 121 column 7
    isosurface at line 157 column 49


Even if the usefulness of the additional output is not really obvious to me, 'isosurface' should not throw an error when this switch is passed.

Btw, Octave does accept the additional parameter "verbose" (ignoring it), when being called with 6+ arguments, e.g.:

[x, y, z] = meshgrid (1:10, 1:10, 1:10);
a = rand (10, 10, 10);
c = rand (10, 10, 10);
isosurface (x, y, z, a, .5, c, "verbose")



2. Empty x, y, z

In Matlab, this line works without error:

isosurface ([], [], [], rand (10, 10, 10), .5)

Otave throws the following error:

error: __marching_cube__: XX, YY, ZZ, C must be matrices of dim 3
error: called from
    __marching_cube__ at line 104 column 5
    isosurface at line 160 column 29


3. x, y, z being vectors

In Matlab, this line works without error:

isosurface (1:10, 1:10, 1:10, rand (10, 10, 10), .5)

Otave throws the following error:

error: __marching_cube__: XX, YY, ZZ, C must be matrices of dim 3
error: called from
    __marching_cube__ at line 104 column 5
    isosurface at line 160 column 29


4. call without iso value

In Matlab, these lines work:

isosurface (rand (10, 10, 10))
[x, y, z] = meshgrid (1:10, 1:10, 1:10);
isosurface (x, y, z, rand (10, 10, 10))

In Octave, the first line throws "error: Invalid call to isosurface" and displays the help text. The last line throws:

error: __marching_cube__: ISO must be scalar value
error: called from
    __marching_cube__ at line 116 column 5
    isosurface at line 160 column 29



If you like I could take a deeper look at these bugs and probably submit a patch.

Markus Mützel <mmuetzel>
Group administrator

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #37696:  isosurface_7.patch added by mmuetzel (28KiB - text/x-diff)
file #37674:  isosurface_6.patch added by mmuetzel (26KiB - text/x-diff)
file #36998:  isosurface_5.patch added by mmuetzel (25KiB - text/x-diff)
file #36360:  isosurface_4.patch added by mmuetzel (20KiB - application/octet-stream)
file #36361:  isosurface.m added by mmuetzel (19KiB - text/x-objcsrc)
file #36278:  isosurface_3.patch added by mmuetzel (21KiB - application/octet-stream)
file #36279:  __calc_isovalue_from_data__.m added by mmuetzel (2KiB - text/x-objective-c)
file #36280:  isosurface.m added by mmuetzel (17KiB - text/x-objective-c)
file #36168:  isosurface_2.diff added by mmuetzel (20KiB - application/octet-stream)
file #36160:  isosurface.patch added by mmuetzel (20KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by philipnienhuis
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mmuetzel (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 20 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-07-04 philipnienhuis StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2016-07-03 mmuetzel Attached File- Added isosurface_7.patch, #37696
    2016-07-02 philipnienhuis StatusPatch Submitted Ready For Test
        Release4.0.0 dev
    2016-07-02 mmuetzel Attached File- Added isosurface_6.patch, #37674
    2016-07-02 lachlan StatusIn Progress Patch Submitted
    2016-07-02 lachlan StatusPatch Submitted In Progress
    2016-06-04 mtmiller Operating SystemMicrosoft Windows Any
    2016-06-04 mtmiller Item GroupSegfault, Bus Error, etc. Incorrect Result
        StatusNone Patch Submitted
    2016-04-26 mmuetzel Attached File- Added isosurface_5.patch, #36998
    2016-02-17 mmuetzel Attached File- Added isosurface_4.patch, #36360
        Attached File- Added isosurface.m, #36361
    2016-02-08 philipnienhuis Carbon-Copy- Added philipnienhuis
    2016-02-08 mmuetzel Attached File- Added isosurface_3.patch, #36278
        Attached File- Added _calc_isovalue_from_data_.m, #36279
        Attached File- Added isosurface.m, #36280
    2016-01-26 mmuetzel Attached File- Added isosurface_2.diff, #36168
    2016-01-25 mmuetzel Attached File- Added isosurface.patch, #36160

    Back to the top

    Powered by Savane 3.14-79a4.
    Corresponding source code