Add voronoi interpolation to HRIR processing

[feliximmohr]

Initial voronoi interpolation implementation as discussed in #151.

[hagenw]

Cool, we are finally on the way to have a proper point selection for 2D and 3D cases.

@feliximmohr could you have a look at the test_interpolation_point_selection() test function that we have for the findconvexcone() method and see if you could integrate the new Voronoi method to compare to the old one? If you are aware of more interesting test cases for the Voronoi method, feel also free to add them there.

[feliximmohr]

test_interpolation_point_selection()

• integrate findvoronoi() including reference data for testing

• add legends to plots

findvoronoi()

• add compatibility for 2D and partial grids as well as collinear test cases

• fix a problem in calculation of surface area leading to negative weights in some cases

• minor fixes

[hagenw]

Great work, I was able to run test_interpolation_point_selection() under Matlab and Octave and it looks very good.

I still have a few questions:

1. For the 3D cases the Voronoi algorithm selects most of the time 4 points, even if it is not very intuitive (see the first plots by the test function). I guess this is this desired behaviour?
2. The weights for the test case "2D grid, arc with gap smaller than 180 deg" are different between findconvexcone and findvoronoi. The one from findconvexcone look a little bit more plausible to me.
3. In general we don't care for the actual weight values in the tests (use NaN instead). Should we change this and add the expected values or isn't it worth the effort?
4. The point selection for the non-sphere case worked quite well. Is this always the case with Voronoi? Then we could disable the corresponding warning.

[fietew]

1. For the 3D cases the Voronoi algorithm selects most of the time 4 points, even if it is not very intuitive (see the first plots by the test function). I guess this is this desired behaviour?

It terms of Voronoi regions, it is plausible.

2. The weights for the test case "2D grid, arc with gap smaller than 180 deg" are different between findconvexcone and findvoronoi. The one from findconvexcone look a little bit more plausible to me.

findconvexcone is based on the euclidean distances between the two selected points: You may guess, that the intersection point of the arrow and the line (not plotted) between the selected points is quite in the middle between the two points. Hence, the weights are almost equal. findvoronoi is based upon great-circle-distances, which in 2D is just difference in azimuth: You can see, that one of the selected point is closer to xs with respect to the azimuth angle. Hence, it gets way more weight.
Note, that both interpolation principles become the same, if the distance between the selected points is small.

3. In general we don't care for the actual weight values in the tests (use NaN instead). Should we change this and add the expected values or isn't it worth the effort?

For some cases, I don't know, what weights exactly are to be expected.

4. The point selection for the non-sphere case worked quite well. Is this always the case with Voronoi? Then we could disable the corresponding warning.

Yes

[hagenw]

We should add an entry for Voronoi in SFS_config.m as well, see https://github.com/sfstoolbox/sfs-matlab/blob/voronoi/SFS_config.m#L384-L393
@fietew could you do this?

Is it also desirable to switch to Voronoi as a new default as it works in 2D and 3D or should we stay with the current default?

[fietew]

FYI: There are still some cases, were I am not completely satisfied with the results of the algorithm, so its still work in progress.

[hagenw]

What is the status of this?
Should we include it, but stay with the current default methods or would you like to improve this first and make it the default after merging?

[fietew]

At the moment. I don't have the time to push this forward.
We should stay with the old version and merge this after we have covered all the special cases.

[fietew]

@feliximmohr will be working on this PR, again. Some thoughts how the code could be (re)structured:

the term "old" corresponds to x0
the term "new" corresponds to [x0;xs]

Detection of dimensionality (to be checked by a function, e.g. check_dimensionality)

• 1D: xs is co-linear with or equal to one x0
• 2D: all x0 and xs are in one plane
• 2.5D: all x0 are in one plane (can be checked via rotate_to_principal_axes), but xs is not
• 3D: otherwise

Interpolation

• 1D: no interpolation needed. We take the co-linear x0.
• 2D:
  • rotate the coordinate system, such that [x0;xs] are in the xy-plane. T
  • order [x0;xs] with respect to the azimuth angle
  • linear interpolation for xs using the neighboring x0
• 2.5D:
  • rotate the coordinate system, such that x0 are in the xy-plane.
  • order x0 with respect to the azimuth angle
  • Calculate new convex hull for [x0;xs]
  • extract all x0 sharing a triangle with xs, denoted as x0_s
  • calculate the difference in azimuth between each point of x0_s and its according neighbors in x0 with respect to the azimuth angle. Multiply the arc length by pi to get the old voronoi area for each x0_s
  • calculate new voronoi region and areas for all x0_s
  • calculate weights
• 3D
  • Calculate old convex hull for x0
  • Calculate new convex hull for [x0;xs]
  • extract all x0 sharing a triangle with xs, denoted as x0_s
  • calculate old and new Voronoi region and areas for all x0_s
  • calculate weights

Remarks on computation of Voronoi regions, if all points are in one hemisphere

• if all points are in one hemisphere, the resulting convex hull does not include the center of the unit sphere
• all circumcenters of the tetrahedrons (triangles of the convex hull augmented by the center of the sphere) are currently projected into the same hemisphere (wrong)
• the orientation of each triangle (is its normal vector pointing away or towards the center of the sphere?) has to be taken into account for the projection

[hagenw]

@feliximmohr: thanks for your effort. I'm wondering what the current status of this pull request is. Are there still changes you would like to make or should I start and review the current code?

[feliximmohr]

I just pushed the last changes i wanted to make. From my point of view, the code review can be started.

[hagenw]

Cool, thanks a lot. I will have a look at it during the next days, hopefully until mid of November.

[hagenw]

Great work, I added a few comments and checking of input args, otherwise all tests run fine and I think it is ready to merge.

The only thing missing at the moment is documentaton about the method in SFS_config in the section before conf.ir.interpolationpointselection. At the moment the its default setting is 'nearestneighbour'. Is the new method better suited as a default, or should we stick with the current method?

[fietew]

There were some functions, which we initially planned to (re-)use for findvoronoi and findconvexcone. However, it turned out that these are either not needed or do not exhibit the exact same behavior. I removed the according functions. As a side effect, this PR does no longer include any changes for findconvexcone.

[fietew]

At the moment the its default setting is 'nearestneighbour'. Is the new method better suited as a default, or should we stick with the current method?

I would stick with the current (simpler) method. The Voronoi interpolation is for sure better as it covers many cases (1D, 2D, etc.) , but quiet complex (also with regard to computation time).

I think, this is ready for merge. Squash would be fine for me.