Bug in interpolate surface

[Roxana-P]

Problem
The routine splinepy.BSpline.interpolate_surface does not work for an unequal number of control points in the two parametric directions.

Example
Please find an example below which illustrates how an easy surface can be interpolated for an equal number of control points, but not for an unequal number.
The code to produce the example is documented below the images.

Equal number of control points (3,3)

Unequal number of control points (3,4)

Unequal number of control points (4,3)

import numpy as np
import splinepy


def interpolate_function(f, samples_per_direction):
    """Interpolate the function f with a given number of samples 
    in each direction on the interval [-1, 1] x [-1, 1]

    Parameters
    ----------
    f : function
        Function to interpolate
    samples_per_direction : list
        Number of samples per direction

    Returns
    -------
    np.ndarray, splinepy.BSpline
        Target points and interpolating BSpline
    """

    # Create dataset
    x = [np.linspace(-1, 1, samples_per_direction[0]), 
        np.linspace(-1, 1, samples_per_direction[1])]
    xx, yy = np.meshgrid(x[0], x[1])

    target_points = np.vstack(
        (xx.flatten(), yy.flatten(), h(xx, yy).flatten())
    ).T

    interpolated_surface = splinepy.BSpline.interpolate_surface(
        target_points,
        size_u=samples_per_direction[0],
        size_v=samples_per_direction[1],
        degree_u=2,
        degree_v=2,
    )
    return target_points, interpolated_surface
    
def h(x, y):
    return x + y

if __name__ == "__main__":
    try:
        import gustaf as gus

        gustaf_available = True
    except ImportError:
        gustaf_available = False


    for sample_sizes in [[3, 3], [3, 4], [4, 3]]:
        target_points, interpolated_surface = interpolate_function(h, sample_sizes)
        
        if gustaf_available:
            bspline = gus.BSpline(**interpolated_surface.todict())
            gus.show.show_vedo([bspline])

Solution
In cpp/splinepy/fitting/fitting.cpp the points are ordered according to their x-coordinate in line 85. There is an index mistake here, which affects the computation only for size_u != size_v.
I am happy to provide a PR with a solution.

[jzwar]

Fixed in PR #100