bugGNU Octave - Bugs: bug #61754, Bicubic interpolation in interp2()...


bug #61754: Bicubic interpolation in interp2() incompatible with Matlab

Submitter:  Christof Kaufmann <christofkaufmann>
Submitted:  Mon 03 Jan 2022 01:25:04 AM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  Christof Kaufmann 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 08 Feb 2022 09:56:51 PM UTC, comment #17: 

If you can, use '-O2' in your CFLAGS and CXXFLAGS.  That provides a good level of optimization without being excessive.

My configure script (written in tcsh so changes to be made if you use a Bourne shell derivative) is

setenv CFLAGS "-O2 -pipe"
# Copy C compiler flags to C++ and Fortran compilers
setenv FFLAGS "${CFLAGS}"
./configure --enable-dependency-tracking -C

Rik <rik5>
Group administrator
Tue 08 Feb 2022 08:25:44 AM UTC, comment #16: 

I tried on my Laptop and got:

  • Matlab: 0.179527 s
  • Octave: 0.687583 s

So there might be some work to do, but I admit that I haven't tried any compilation optimizations. So, the bottleneck could also be partially elsewhere. Don't know.

Christof Kaufmann <christofkaufmann>
Tue 08 Feb 2022 05:30:10 AM UTC, comment #15: 

@Christof: Good analysis.  A loop count of 4 is minimal and shouldn't be a problem.  I was lazy and googled "cubic interpolation FIR filter" and saw some academic papers so I thought there might be a shortcut.  I tried the test code in comment #13 in the online Matlab interface and it only took 0.23 seconds (~4x faster), but I don't know how powerful the machine they use is.

I'll close the report as fixed.

Rik <rik5>
Group administrator
Mon 07 Feb 2022 07:50:36 AM UTC, comment #14: 

Thanks for the corrections.

The for loop has got only 4 iterations with big operations.  I tried to vectorize it, but it was not easy to do.  I would have the shifts in the third axis, but when selecting rows (or columns) of Z the shape is not preserved.  Reshaping would require transposing in one of the cases, if I did not overlook something.  My conclusion was, that this is unlikely to be faster, but complicates the code quite a lot.

Another optimization possibility could be the padding.  Padding is not always necessary.  It is only required if there are points near the boundary.  So this is not so much about performance, but more of the memory footprint, since it might avoid one copy of `Z`.  Thinking further, the full copy could always be avoided, when dealing with the near-boundary-points separately.  Then only the three first or last rows or columns had to be copied and one row or column padded, which would be cheaper.  The special case of copying the last row or column could also be caught, instead of using the algorithm with double padding there.  However, these are a lot of special cases for the code and the improvement might be rather low.

Regarding your suggestion with conv2() or filter2():  No, only if the interpolation is a pure translation (and optional crop), it is a real convolution.  In this case `h` would contain the same value, since it is the distance to the point before.  Usually it is not a convolution, or, it may be considered as a convolution with a different kernel everywhere.  So it is just similar to a convolution and conv2() or filter2() cannot be applied.

Christof Kaufmann <christofkaufmann>
Mon 07 Feb 2022 12:15:23 AM UTC, comment #13: 

I made a few small changes to follow Octave coding conventions here (http://hg.savannah.gnu.org/hgweb/octave/rev/95ee9b9605cc).

I did a minimal check on performance with

z = rand (1e3, 1e3);
tic; zi = interp2 (z, 2, "cubic"); toc

Result was 0.8 seconds.  Seems fine for me, but I don't know what size datasets people usually use.  I think I did notice a for loop which might be slow.  Since there is a kernel being applied would conv2() or filter2() be faster alternatives?

Rik <rik5>
Group administrator
Sun 06 Feb 2022 04:22:45 PM UTC, comment #12: 

Feel free to open a new report if there isn't one already.

Markus Mützel <mmuetzel>
Group administrator
Sun 06 Feb 2022 03:52:21 PM UTC, comment #11: 

Thank you!

OK, I will try to remember to summarize the file changes in the
next patches.

I search for bugs containing "spline" and could not find appropriate bug reports. I might have missed them though.

Christof Kaufmann <christofkaufmann>
Sun 06 Feb 2022 02:00:41 PM UTC, comment #10: 

Sorry. I lost track over which change was tracked where, and what was already applied or not...

I pushed your patch with some minor changes here:

It's already a very big help that those are "real" hg patches. Thank you for that.
It's a convention for GNU projects to summarize (very briefly) in the commit message what has changed in which file. It would be even more helpful if you could add that, too. (Not a requirement. It'd just help us out a bit when reviewing/pushing. 😊)

Regarding the "xtest"s: Afaict, those are checking for the results that Matlab produces for these test cases. I left them in as is in this change. But it would be nice if there was a bug report that would track those. When there is one, the "xtest"s should be replaced by "test" with that bug id.

Markus Mützel <mmuetzel>
Group administrator
Wed 02 Feb 2022 08:49:03 PM UTC, comment #9: 

I don't think it got reviewed yet. I am not sure about the `xtest`s I have added for the spline method. Maybe they should be removed again an added in a separate patch, but they are there for consistency.

Besides the patch in comment #7, you might also want to have a look at the PDF from comment #2. It tries to explain the approach and might help at reviewing the code.

Christof Kaufmann <christofkaufmann>
Wed 02 Feb 2022 08:19:37 PM UTC, comment #8: 

Did the latest patch get reviewed?  I know there was a series of other patches from Christof dealing with _pchip_deriv_ and this bug report may have been forgotten.

Rik <rik5>
Group administrator
Sun 30 Jan 2022 11:57:32 PM UTC, comment #7: 

Rebased on current default and improved performance.

(file #52777)

Christof Kaufmann <christofkaufmann>
Sat 22 Jan 2022 10:26:11 AM UTC, comment #6: 

Rebased on new patch of bug #61863.

(file #52707)

Christof Kaufmann <christofkaufmann>
Thu 20 Jan 2022 11:03:23 AM UTC, comment #5: 

Sorry, I forgot to change the documentation about "cubic". Updated patch attached (still rebased onto the patch of bug #61863).

(file #52699)

Christof Kaufmann <christofkaufmann>
Thu 20 Jan 2022 02:46:56 AM UTC, comment #4: 

OK, I submitted a new bug report (bug #61863) for the complex matrix issue. While doing this, I found out, that "pchip" needed some extra care. I also added a patch there.

Then I rebased the bicubic interpolation fix on the patch of bug #61863 with some minor improvements compared to the submission in comment #2. The new patch is attached.

(file #52694)

Christof Kaufmann <christofkaufmann>
Wed 19 Jan 2022 04:54:11 PM UTC, comment #3: 

Thank you very much for your work.

Instead of modified source files, patch files are the preferred form of contribution. That way, it would be easier to see what you changed and it would reduce the risk of accidentally reverting changes in the development sources (by replacing an existing file with some other file).

Additionally, it would be preferable if bug fixes and new features come in distinct patch files. (You can still use the same bug report for submitting those when they are closely related.)

See also: https://wiki.octave.org/Contribute

Markus Mützel <mmuetzel>
Group administrator
Wed 19 Jan 2022 04:43:51 PM UTC, comment #2: 

I implemented bicubic interpolation in interp2.m, see attached file. To make the implementation easier to understand, I also attached a PDF (LaTeX source) to explain the approach. I also added tests, not only for the bicubic interpolation, but for the sake of consistency also for the existing methods. I found a few bugs:

  • Interpolating complex matrixes gave an NA + 0i instead of NA + NAi. I could fix this.
  • Spline interpolation does only support meshgrid format, but not index pairs. Matlab supports that. There is already a note in _splinen_.m. I added a failing test using xtest.
  • Spline interpolation returns a row, when interpolating a single column, though with the correct values. I added a failing test using xtest.

I have not opened new bug reports nor checked for existing ones regarding the two spline interpolation issues.

Where to go from here?

  • I have not checked interp1.m, interp3.m, interpn.m. The cubic interpolation might also be Matlab incompatible.
  • Bilinear interpolation could also be implemented by convolution easily. This might be faster, maybe only for specific situations. When switching from interp2 to the convolutional approach in imresize.m to support anti-aliasing (submitted in bug #61760, not yet pushed), the speed-up was quite big for some tests.

Please let me know if the implementation needs improvements.

(file #52682, file #52683)

Christof Kaufmann <christofkaufmann>
Mon 03 Jan 2022 01:47:23 PM UTC, comment #1: 

Patches to fix compatibility are always welcome.

Markus Mützel <mmuetzel>
Group administrator
Mon 03 Jan 2022 01:25:04 AM UTC, original submission:  

Using "cubic" for the method argument of interp2 is just an alias for "pchip" (piecewise cubic Hermite interpolating polynomial). However in Matlab a cubic convolution is used, as stated in their documentation of interp2. This method is also described on Wikipedia's Bicubic Interpolation page.

A tiny example to show the difference:

> interp2(magic(4), 1, 1.5, 'cubic')
ans = 8.1875

Matlab (R2021a):

> interp2(magic(4), 1, 1.5, 'cubic')
ans = 8.6250

The topic has already been discussed in bug #60322 (regarding imresize) and there is even a patch for imresize available to reach Matlab compatibility. Maybe this patch can be adapted for interp2 as well? What is your opinion on this?

Christof Kaufmann <christofkaufmann>


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

Attach Files:

Attached Files
file #52777:  interp2_fix_bicubic4.patch added by christofkaufmann (13KiB - text/x-patch - bicubic interpolation - v4 (rebased and performance improved))
file #52707:  interp2_fix_bicubic3.diff added by christofkaufmann (14KiB - text/x-patch - bicubic interpolation - v3 (rebased))
file #52699:  interp2_fix_bicubic2.diff added by christofkaufmann (14KiB - text/x-patch - bicubic interpolation - v3)
file #52694:  interp2_fix_bicubic1.diff added by christofkaufmann (13KiB - text/x-patch - bicubic interpolation - v2)
file #52682:  interp2.m added by christofkaufmann (26KiB - text/x-objcsrc - bicubic interpolation - v1)
file #52683:  bicubic.pdf added by christofkaufmann (104KiB - application/pdf - bicubic interpolation - v1)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by christofkaufmann (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 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-02-08 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-02-06 mmuetzel StatusPatch Submitted Ready For Test
    2022-02-02 rik5 StatusNone Patch Submitted
    2022-01-30 christofkaufmann Attached File- Added interp2_fix_bicubic4.patch, #52777
    2022-01-22 christofkaufmann Attached File- Added interp2_fix_bicubic3.diff, #52707
    2022-01-20 christofkaufmann Attached File- Added interp2_fix_bicubic2.diff, #52699
    2022-01-20 christofkaufmann Attached File- Added interp2_fix_bicubic1.diff, #52694
    2022-01-19 christofkaufmann Attached File- Added interp2.m, #52682
        Attached File- Added bicubic.pdf, #52683

    Back to the top

    Powered by Savane 3.13-8484.
    Corresponding source code