HRIR interpolation point selection

[hagenw]

I'm not sure if the simple usage of nearest neighbor points for performing the interpolation is the best approach for impulse responses. Here is an easy example:

x0 = [45 46 46 50; 0 -2 2 0];
xs = [47; 0];
x0_int = findnearestneighbour(x0,xs,3);
figure; plot(x0(1,:),x0(2,:),'xb')
axis square; axis([43 52 -3 3]); hold on
plot(xs(1),xs(2),'xr')
plot(x0_int(1,:),x0_int(2,:),'or','MarkerSize',10)
xlabel('Azimuth / deg')
ylabel('Elevation / deg')

As you can see we want to get the point shown as red cross. The red circles are the three nearest neighbors that are selected for the interpolation. Afterwards it is checked if the nearest and second nearest neighbor have the same azimuth or elevation. If so only these points are used. In this example all three points will be used for interpolation.
The problem is that in this case it would in my opinion make more sense to use just the two points that have both an elevation of 0°. The problem is I don't know if we are able to write an algorithm that does this sort of clever selection that works with every grid and point spacing.
Has anyone some experience on this?

[fietew]

Maybe this is interesting:
https://en.wikipedia.org/wiki/Natural_neighbor
the German version explains it a little bit better:
https://de.wikipedia.org/wiki/Voronoi-Interpolation

[fietew]

Voronoi regions for a sphere:
https://github.com/tylerjereddy/py_sphere_Voronoi
https://www.jasondavies.com/maps/voronoi/

[hagenw]

Cool, especially the second link.

Now, we only need a person who is going to implement this.

[trettberg]

Pushed some code to a branch named vbap_style_interpolation.

Here is the example from above.
(The upper point actually gets a weight of zero. Just the markersize has an offset.)

azi = [0; 1; 1; 5;]*pi/180;
elev = [0; 1 ; -1; 0;]*pi/180; 
[x,y,z] = sph2cart(azi,elev,ones(4,1));
x0 = [x y z];
[x,y,z] = sph2cart(2*pi/180,0,1);
xs = [x y z];

[indeces,weights] = findconvexcone(x0,xs);
x0_int = x0(indeces,:);

scatter3(x0(:,1),x0(:,2),x0(:,3),[],[],'bo');
markersize = 70*weights+30*ones(3,1)
hold on
scatter3(xs(:,1),xs(:,2),xs(:,3),[],[],'rx')
scatter3(x0_int(:,1),x0_int(:,2),x0_int(:,3), markersize ,[],'ro');
xlabel('y'); zlabel('z');
view([90,0]); grid off;

[fietew]

Is there literature describing the interpolation scheme in detail?

[hagenw]

That looks nice. I guess it would be good to specify several interpolation test cases to be able to compare different approaches.

[trettberg]

(This is in fact VBAP if all points lie on a sphere.)

[VeraE]

I just played around with Till's idea, which works nicely. I still wonder which points we really want to use for interpolation, though. See the picture below (just a modification of Till's example): In both cases the algorithm uses the three points marked in red but the upper point has a weight of zero. Does this selection make sense in the right case as the upper point is so close to the desired position?

Another question is which three points (forming a triangle containing the desired point) are selected by the algorithm if the desired point lies on an intersection of several triangle edges. For example if the desired point in the pictures below would lie exactly between the upper and the lower point. Till is just checking on that.

[fietew]

Very interesting:
https://www.youtube.com/watch?v=gxNa9BD5CnQ
https://www.youtube.com/watch?v=zaGd5tXkCnE&index=35&list=PLGVZCDnMOq0qfJkoiFj-hN7lSHgQzXtqQ

[trettberg]

Sorry, I made a mistake in these examples:
I compute the triangulation via the 3D convex hull. For this to work, the origin must be contained in the convex hull. That is not the case here.
Done properly, the upper, lower and rightmost point are selected in this example.

Still it may be a worthwhile improvement over the current version: The three nearest neighbors are selected such that the point lies within the triangle.

[hagenw]

I created a pull request for this.
As a next step it would be very cool, if someone could define a few test cases in a function called test_interpolation_point_selection.m in the validation folder. Otherwise it is not so easy to say if this is finished now or if we have to change something.

[hagenw]

I have run the test function and everything works fine.
The good thing is that we also should be able to implement the interpolation function in a general way as findconvexcone() returns always three loudspeakers, but sets the weights to zeros if not all are necessary for interpolation.
I will first do a new release of the Toolbox, after that I will merge #102 into master.

[hagenw]

In findconvexcone() you state that it may fail when

The convex hull of x0 does not contain the origin.

which I'm not completely understand. Can this happen for typical loudspeaker setups? If so, could you please add an example for such a case to the test function as well.
The other restriction that the x0 needs to be convex should be fine as we assume it for WFS as well.

[trettberg]

As an example, let's assume a grid which covers only a spherical cap around the northpole (of the unit sphere).
Then origin is outside the convex hull.

Two kinds of errors may occur:

a) requested point: northpole.
A triple of points is returned. It will work for interpolation but the triangle is not necessarily part of a spherical Delaunay triangulation.
(This is the case the examples discussed above.)

b) requested point: southpole.
This is "outside" the grid. No triple can be found. An error is raised.

I added testcases for these.

The workaround: add dummy points
Even if the origin is contained in the hull, this will be desirable in practice for a grid with "large holes".

[hagenw]

Thanks for all test cases, I added now findconvexcone() by merging #102 and think this is a nice addition as we can also use it for VBAP if we plan to add this at some point.

For the actual HRTF interpolation we are still free to choose how to do it and should also carefully have a look at all the test cases where findconvexcone() fails.
@fietew: Matlab and Octave both contain voronoi() and voronoin(), so it should also be possible to try such kind of interpolation, but I guess it will have the same problems for the failing test cases.

The update of the interpolation method could be done together with #102 or done in a new pull request afterwards as #102 incorporates a more special case of interpolation that works only for some HRTF data sets that have no noise in the phase at low frequencies.

[fietew]

@fietew: Matlab and Octave both contain voronoi() and voronoin(), so it should also be possible to try such kind of interpolation, but I guess it will have the same problems for the failing test cases.

As far as I know, this functions compute the voronoi regions for cartesian coordinates. We are looking for the voronoi interpolation on the surface of a sphere, which requires to use a different metric. I however agree that we can postpone this implementation and create a new pull request, when it's ready.

[hagenw]

Is this issue covered by #130 and we can close it?

Our implementations do not use Voronoi, would this still be better?
If so, I would propose to create a new issue requesting Voronoi-based point selection.

[VeraE]

I'm not sure which one is better. In some cases Voronoi interpolation makes intuitively more sense to me, e.g. for larger gaps in a dataset. I don't see myself trying out Voronoi interpolation in the near future to compare it, though. Looks more a research question to me, though I might not be aware of all relevant literature.
If we close this now (which would be alright with me), we should maybe remove the voronoi option in the code for #130 because it looks like it won't be implemented in the near future.