bugGNU Octave - Bugs: bug #33539, griddata does not support 'v4'...


bug #33539: griddata does not support 'v4' method

Submitter:  Arnaud Delorme <arnodelorme>
Submitted:  Sat 11 Jun 2011 03:45:56 PM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  3 - Low Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
* Mandatory Fields

Add a New Comment (Rich Markup)

Comment Type & Canned Response:


( Jump to the original submission )

Mon 13 Apr 2020 04:12:46 PM UTC, comment #14: 

Checked in one more change to issue a "not implemented" message for "cubic" method (https://hg.savannah.gnu.org/hgweb/octave/rev/bb50407f9c60).

Rik <rik5>
Project Administrator
Mon 13 Apr 2020 02:24:46 AM UTC, comment #13: 

thanks Rik. Only issue with your error message on griddatan is that cubic hasn't actually been implemented on griddata yen.

Nicholas Jankowski <nrjank>
Project Member
Mon 13 Apr 2020 01:30:24 AM UTC, comment #12: 

I reviewd the patch and checked it in here https://hg.savannah.gnu.org/hgweb/octave/rev/bb929d5a34cb.

I also reworked the input validation and the documentation for griddatan in this changeset https://hg.savannah.gnu.org/hgweb/octave/rev/289882040316.

Marking as fixed and closing report.

Rik <rik5>
Project Administrator
Thu 09 Apr 2020 08:24:19 PM UTC, comment #11: 

forgot to mention patch also adjusts the error messages in griddatan for calling an invalid method add added some input validation BISTS.

Nicholas Jankowski <nrjank>
Project Member
Thu 09 Apr 2020 08:23:05 PM UTC, comment #10: 

thanks, the in-comments explanation of the v4 method is much better. I was honestly unaware of .', so that makes sense too. I included those edits, added a couple very simple tests to make sure linear, v4, and nearest all get tested, and removed the 'error if 'v4' called' test.  Also added the method note to the line about griddata compatibility in NEWS.

Patch attached. if all is well it should close out this bug. Other methods could get handled elsewhere if someone has a compatible implementation for them. There's already another bug for 'cubic' (bug #35178). Haven't seen one for 'natural' but could create one or expand the cubic one to include it to round out the set.

(file #48795, file #48796)

Nicholas Jankowski <nrjank>
Project Member
Wed 25 Mar 2020 09:25:37 PM UTC, comment #9: 

I added some comment proposal. Perhaps too much, so don't hesitate to simply take what you want.

Code looks well thought. I know we are talking about a real valued domain but you may consider replacing ' by .'

(file #48676)

Hg200 <hg200>
Mon 23 Mar 2020 01:29:33 PM UTC, comment #8: 

not sure when I'll get to finishing it up, here's a version of griddata with the v4 implementation in case anyone wants it ahead of a patch. still missing BISTs for that (and really any) input method.  if anyone has simple BISTS for simple or matlab compatible verification I'll add them.

(file #48655)

Nicholas Jankowski <nrjank>
Project Member
Tue 17 Mar 2020 07:59:34 PM UTC, comment #7: 

ok, i now have a v4 implementation using Mike's reference and one a bit more recent that elucidated the algorithm a bit for me.

noting that cubic and natural are still unimplemented, I'll plan to put one patch together to cover them all.

for now, here's the v4 codeblock:

elseif (strcmp (method, "v4"))
    ## Biharmonic Spline as per D. Sandwell 1987 and Deng & Tang  2011
    ## 2D greens function - G(x) = |x|^2 * (ln|x| -1)
       #wi = sum_j_N alpha_j * G(xi-xj)
       #   = sum_j_N alpha_j * [ |xi-xj|^2 * (ln|xi-xj| -1)]

    ##step1 - find weight function coefficients alpha(j) from point distances:
    r = sqrt ((x - x').^2 + (y - y').^2);
    D = (r.^2) .* (log (r) - 1);
    D(isnan (D)) = 0;  ##correct for 0*inf =NaN when x1-x2 = 0

    ## Solve linear eqn #W(Xi) = Zi = ALPHA * G -> PHI = inv(ALPHA)*Z = PHI \ Z
    alpha_j = D \ z;

    ##step2 - use alphas and greens functions to get interpolated points
    ## use dim3 projection for vectorized calculation to avoid loops.
    ## TODO: if this approach is too memory intensive, revert portion to loop
    x = permute (x, [3, 2, 1]);
    y = permute (y, [3, 2, 1]);
    alpha_j = permute (alpha_j, [3, 2, 1]);
    r_i = sqrt ((xi - x).^2 + (yi - y).^2);
    Di(isnan (Di)) = 0; ##correct for 0*inf =NaN when x1-x2 = 0

    zi = sum (Di .* alpha_j, 3);

max error against a matlab test is on the order of 1e-15, so I think I'll call that good until someone find an error. the whole function is missing BISTs too.

Nicholas Jankowski <nrjank>
Project Member
Tue 10 Mar 2020 10:34:26 PM UTC, comment #6: 

Awesome, thanks mike. any chance you have references for the other two.missing method algorithms? I spent a minutes looking for natural neighbor method but didn't found an easy license compatible le one. I think the other is a cubic spline?

Was gonna put these on my eventual todo list. But might add a 'short projects' description in the wiki first.

Nicholas Jankowski <nrjank>
Project Member
Tue 10 Mar 2020 08:59:33 PM UTC, comment #5: 

I think someone will have to implement this interpolation from the reference. FWIW I did some searching and found an older cached copy of the griddata documentation. The reference for the "v4" method is Sandwell, David T., "Biharmonic Spline Interpolation of GEOS-3 and SEASAT Altimeter Data", Geophysical Research Letters, 2, 139-142,1987.

As an aside, tracing the historical record of Octave's griddata is a fun exercise. It was originally added to Octave in 2007, a few months before version 3.0.0 (https://hg.savannah.gnu.org/hgweb/octave/rev/9fddcc586065). It was pulled from the Octave Forge geometry functions. And the original version of that function was written between 1999-2001, the initial version is here (http://hg.code.sf.net/p/octave/geometry/file/3c2d81569bf2/griddata.m).

Mike Miller <mtmiller>
Project Administrator
Tue 10 Mar 2020 08:42:48 PM UTC, comment #4: 

according to the matlab help, v4 uses "Biharmonic spline interpolation (MATLAB® 4 griddata method) supporting 2-D interpolation only. Unlike the other methods, this interpolation is not based on a triangulation."

Any chance griddata existed back in Octave 0.8 and was compatible? otherwise anyone suggest a biharmonic algorithm most likely to be v4 compatible? the help doesn't have a library reference for the algorithm. I haven't tried Archive.org to see if it was listed in old doc versions.

Nicholas Jankowski <nrjank>
Project Member
Tue 10 Mar 2020 07:48:17 PM UTC, comment #3: 

minor update: in another griddata bug (bug #57323), the pending patch updates the help text to indicate that Octave only supports the options linear and nearest, and a 'not yet implemented' error appears if methods v4, cubic, or natural are called.

Nicholas Jankowski <nrjank>
Project Member
Thu 24 Nov 2016 12:42:21 AM UTC, comment #2: 

This behavior is still present in Octave 4.3.0+

Juan Pablo Carbajal <juanpi>
Project Member
Thu 16 Feb 2012 10:53:58 PM UTC, comment #1: 

griddata has been expanded to accept vectors of any orientation in this changeset (http://hg.savannah.gnu.org/hgweb/octave/rev/5736d93b22d0).  The changes will show up in release 3.8 or the reporter can build from Mercurial sources.

The bug report is being re-titled to reflect the other named issue, that griddata does not accept the 'v4' argument.  See the documentation for griddata (http://www.mathworks.com/help/techdoc/ref/griddata.html)

Rik <rik5>
Project Administrator
Sat 11 Jun 2011 03:45:56 PM UTC, original submission:  

gridadata does not perform automatic transposition of arrays 

[Xi,Yi,Zi] = griddata(rand(1,10),rand(1,10),rand(1,10)',rand(1,100)',rand(1,100)); 
% this works under Matlab but crashes nuder Octave

Also griddata does not take as input the 'invdist' option

Arnaud Delorme <arnodelorme>


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

Attach Files:

Attached Files
file #48795:  griddata_v4_patch.diff added by nrjank (7KiB - application/octet-stream - griddata with v4 method)
file #48796:  griddata.m added by nrjank (11KiB - text/plain - griddata with v4 method)
file #48676:  griddata_v2.m added by hg200 (10KiB - text/x-objcsrc - Comment proposals)
file #48655:  griddata.m added by nrjank (10KiB - text/plain - modified griddata with 'v4' biharmonic interpolation method added)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by hg200 (Updated the item)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by juanpi (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by arnodelorme (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 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-04-13 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2020-04-09 nrjank Attached File- Added griddata_v4_patch.diff, #48795
        Attached File- Added griddata.m, #48796
    2020-03-25 hg200 Attached File- Added griddata_v2.m, #48676
    2020-03-23 nrjank Attached File- Added griddata.m, #48655
    2020-03-10 mtmiller Carbon-CopyRemoved 80942 -
    2014-01-19 mtmiller CategoryNone Octave Function
        Release3.4.0 dev
        Operating SystemMac OS Any
    2012-02-27 rik5 StatusNone Confirmed
    2012-02-16 rik5 Summarygriddata Matlab compatibility issue griddata does not support 'v4' method
    2012-02-16 rik5 Priority5 - Normal 3 - Low
        Item GroupIncorrect Result Matlab Compatibility

    Back to the top

    Powered by Savane 3.9