bugGNU Octave - Bugs: bug #61961, `griddata` doesn't emit a warning...

 
 

bug #61961: `griddata` doesn't emit a warning for duplicate data points

Submitter:  Razdvizh <razdvizh>
Submitted:  Sat 29 Jan 2022 09:27:54 AM UTC
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Missed Error or Warning
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
* Mandatory Fields

Add a New Comment (Rich Markup)
   

Comment Type & Canned Response:
       

 

( Jump to the original submission )

Mon 31 Jan 2022 11:20:41 PM UTC, comment #16: 

comment #15:

> I'm not sure I understand what you are trying to do. But maybe take a look at `demo trimesh` if you'd like to plot a surface on an irregular grid.

Thank you. tetramesh works


Razdvizh <razdvizh>
Mon 31 Jan 2022 04:45:16 PM UTC, comment #15: 

I'm not sure I understand what you are trying to do. But maybe take a look at `demo trimesh` if you'd like to plot a surface on an irregular grid.

Re-titling this for the missing warning of `griddata` which might have clarified that that function cannot be used for that task. Like its name suggests, it's meant to be used for plotting on a (regular) grid.

Markus M├╝tzel <mmuetzel>
Project Member
Mon 31 Jan 2022 04:08:07 PM UTC, comment #14: 

comment #13:

> I think the approach in `griddata`, based on triangulation, is not adapted to interpolation on closed surfaces. Matlab doesn't support it either and returns a meaningful message: "Duplicate data points have been detected and removed - corresponding values have been averaged".
>
> Can we close this report as the failure is expected, or should we make it a feature request for a warning when duplicate data points are present?


Close, add note to doc. Do you know other approaches, software perhaps?

Razdvizh <razdvizh>
Mon 31 Jan 2022 03:43:49 PM UTC, comment #13: 

I think the approach in `griddata`, based on triangulation, is not adapted to interpolation on closed surfaces. Matlab doesn't support it either and returns a meaningful message: "Duplicate data points have been detected and removed - corresponding values have been averaged".

Can we close this report as the failure is expected, or should we make it a feature request for a warning when duplicate data points are present?

Pantxo Diribarne <pantxo>
Project Member
Mon 31 Jan 2022 03:29:41 AM UTC, comment #12: 

Same number of vertex, need more general way
cylinder([1:5 5 4 4:-2:-2])


Razdvizh <razdvizh>
Mon 31 Jan 2022 03:21:21 AM UTC, comment #11: 

In general if

f(x,y,z) = C

then solve for z such that

z = C - f2(x,y)

To plot, choose x,y in the 2-D plane that are valid using functions like meshgrid, ndgrid, or something of your own choosing.  Calculate z for those x and y and finally plot with mesh (x,y,z).

Rik <rik5>
Project Administrator
Mon 31 Jan 2022 02:54:34 AM UTC, comment #10: 

I would look at cylinder.m again.  For example, a right circular cone can be had by invoking the function as

cylinder ([1 0])

file:cone.png


Rik <rik5>
Project Administrator
Sun 30 Jan 2022 09:42:59 PM UTC, comment #9: 

For cylinder it's simple
For Hyperboloid z-cut differs by size, so you need always the same number points. Is there anywhere example with cone with different size of points on borders?

Razdvizh <razdvizh>
Sun 30 Jan 2022 09:26:28 PM UTC, comment #8: 

comment #6:

> > Of course it can, it's scattered data. Is there another method which works well with duplicate z data?
>
> This makes no sense.  The translation of the English language statement to a concrete mathematical statement is that
>

> V1 = f(x1,y1,z1)
> V2 = f(x1,y1,z1)

Hmm, so the problem griddata is design for explicit function
z=f(x,y)
I need implicit f(x,y,z)=const
I looked at cylinder.m
Duplicate some data and griddata is not needed at all?
 for comment #3
xx=coord(:,1);xx(:,2)=xx;
yy=coord(:,2);yy(:,2)=yy;
zz(1:2*columns(t),1)=-1;zz(1:2*columns(t),2)=1;
 mesh(xx,yy,zz);

Razdvizh <razdvizh>
Sun 30 Jan 2022 07:55:12 PM UTC, comment #7: 

From the discussion so far, the following approaches might be more useful than griddata for this specific application:

  • convhulln (will only use outermost vertices, discarding interior ones in case of duplicated coordinates)
  • kriging toolbox (pkg load stk), useful to handle data with missing or duplicate points or inconsistent repetitions.
  • manually fitting a surface, such as a plane or polynomial z = f(x,y) and plotting that.

Do these options offer any help?

Arun Giridhar <arungiridhar>
Project Member
Sun 30 Jan 2022 07:08:30 PM UTC, comment #6: 

> Of course it can, it's scattered data. Is there another method which works well with duplicate z data?


This makes no sense.  The translation of the English language statement to a concrete mathematical statement is that

V1 = f(x1,y1,z1)
V2 = f(x1,y1,z1)

What function produces different values for identical inputs and is also useful?

As Pantxo pointed out, griddata called as

zz = griddata (coord(:,1),coord(:,2),coord(:,3),xx,yy);

is doing interpolation over a 2-D surface.  If you are trying to do interpolation over a 3-D volume then the calling form is

VI = griddata (X, Y, Z, V, XI, YI, ZI)

This would be appropriate, for example, for temperature data where you might have data sampled at the same height (same Z value) but on the left-hand side of the room (low X value) and on the right-hand side of the room (high X value).

For example of a cylinder with mesh walls see cylinder.m function included with Octave.

Rik <rik5>
Project Administrator
Sun 30 Jan 2022 03:22:00 PM UTC, comment #5: 

@pantxo: Do you know algorithm to produce mesh between closed curves?
Mesh for cylinder lying on walls is kinda wrong


Razdvizh <razdvizh>
Sun 30 Jan 2022 09:50:45 AM UTC, comment #4: 

@razdvitch: When I read the doc for griddata, I see that the function first interpolates your data with a function f defined as z = f (x,y) and then uses evaluates that function on your regular (xx,yy) grid. AFAIU this approach requires that there is an unique value z for a given coordinate (x,y). A cylinder has an infinity of z coordinates for a given (x,y) so I think griddata is a no go. 

Pantxo Diribarne <pantxo>
Project Member
Sun 30 Jan 2022 07:39:19 AM UTC, comment #3: 

comment #1:

> Is there an absolute minimum viable example here?


clear all;clf

t=-2*pi:.1:2*pi;
coord(1:columns(t),1)=cos(t);
coord(1:columns(t),2)=sin(t);
coord(1:columns(t),3)=-1;

coord(1+columns(t):2*columns(t),1)=cos(t);
coord(1+columns(t):2*columns(t),2)=sin(t);
coord(1+columns(t):2*columns(t),3)=1;

limits = [min([coord(:,1)]), max([coord(:,1)]) ,...
min([coord(:,2)]), max([coord(:,2)]),...
min([coord(:,3)]), max([coord(:,3)])];
xlabel("x");ylabel("y");zlabel("z");

view(3);

xlabel("x");ylabel("y");zlabel("z");
plot3(coord(:,1),coord(:,2),coord(:,3));
axis ("equal",limits);

hold on;

 [xx,yy] = meshgrid (linspace (min(coord(:,1)),max(coord(:,1)),200));
  zz = griddata (coord(:,1),coord(:,2),coord(:,3),xx,yy);
 mesh(xx,yy,zz);
 hidden('off');


Razdvizh <razdvizh>
Sun 30 Jan 2022 07:10:08 AM UTC, comment #2: 

comment #1:

> It's not clear what "vertical vertices" means.  The initial data you want to interpolate can't have duplicate Z-values for a given X and Y point.  How would the interpolation know which Z-value you wanted?

Of course it can, it's scattered data. Is there another method which works well with duplicate z data?

Razdvizh <razdvizh>
Sat 29 Jan 2022 11:45:37 PM UTC, comment #1: 

It's not clear what "vertical vertices" means.  The initial data you want to interpolate can't have duplicate Z-values for a given X and Y point.  How would the interpolation know which Z-value you wanted?

Is there an absolute minimum viable example here?

I tried this

[xx, yy] = meshgrid (0:1);
zz = xx;
mesh (xx,yy,zz);  # for visualization, this is an inclined plane
## Try interpolating a Z-value at an existing vertex
z = griddata (xx, yy, zz, 1, 1)

which works in that the answer z is exactly 1 as expected.

Rik <rik5>
Project Administrator
Sat 29 Jan 2022 09:27:54 AM UTC, original submission:  

With some slope you could get decent result, but if slope became little or none you've got blanks

Example is for cylinder.
See also https://octave.discourse.group/t/strange-result-with-griddata/2147

clear all;clf
f2=fopen("c5.bin","rb+");coord=fread(f2,[47716,3],"double");fclose(f2);
f2=fopen("c6.bin","rb+");coordnum=fread(f2,[152,1],"double");fclose(f2);

limits = [min([coord(:,1)]), max([coord(:,1)]) ,...
min([coord(:,2)]), max([coord(:,2)]),...
min([coord(:,3)]), max([coord(:,3)])];
xlabel("x");ylabel("y");zlabel("z");

invs=@(n) line("xdata",coord(1+coordnum(n):coordnum(n+1),1), ...
"ydata",coord(1+coordnum(n):coordnum(n+1),2), ...
"zdata",coord(1+coordnum(n):coordnum(n+1),3),"color", "b", "marker", "none");

axis ("equal",limits);
view(3);
xlabel("x");ylabel("y");zlabel("z");

g_step=20;
Nr=1:g_step:rows(coordnum);
hl1=arrayfun(invs, Nr);
hold on;

 [xx,yy] = meshgrid (linspace (min(coord(:,1)),max(coord(:,1)),200));
  zz = griddata (coord(:,1),coord(:,2),coord(:,3),xx,yy);
 mesh(xx,yy,zz);
 hidden('off');

Razdvizh <razdvizh>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #52785:  s11.png added by razdvizh (54KiB - image/png)
file #52786:  s12.png added by razdvizh (26KiB - image/png)
file #52780:  s10.png added by razdvizh (33KiB - image/png)
file #52779:  cone.png added by rik5 (63KiB - image/png)
file #52776:  s9.png added by razdvizh (144KiB - image/png)
file #52775:  s8.png added by razdvizh (21KiB - image/png)
file #52772:  s6.png added by razdvizh (51KiB - image/png)
file #52773:  c6.bin added by razdvizh (1KiB - application/octet-stream)
file #52774:  c5.bin added by razdvizh (1MiB - 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 mmuetzel (Posted a comment)
  • -email is unavailable- added by arungiridhar (Posted a comment)
  • -email is unavailable- added by pantxo (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by razdvizh (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 project members can vote.

     

     

     

     

    Follow 13 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-01-31 razdvizh Attached File- Added s11.png, #52785
        Attached File- Added s12.png, #52786
    2022-01-31 mmuetzel Item GroupNone Missed Error or Warning
        StatusNone Confirmed
        Release5.2.0 dev
        SummaryGriddata cannot produce mesh on vertical vertices `griddata` doesn't emit a warning for duplicate data points
    2022-01-31 razdvizh Attached File- Added s10.png, #52780
    2022-01-31 rik5 Attached File- Added cone.png, #52779
    2022-01-30 razdvizh Attached File- Added s9.png, #52776
    2022-01-30 razdvizh Attached File- Added s8.png, #52775
    2022-01-29 razdvizh Attached File- Added s6.png, #52772
        Attached File- Added c6.bin, #52773
        Attached File- Added c5.bin, #52774

    Back to the top


    Powered by Savane 3.9