Non-identical results if simplices is provided or not

[kanglcn]

I have a point cloud interferogram intf (with the coordinates gix) to be unwrapped

intf.shape # (292623,) 
gix.shape # (292623, 2)

The edges and simplices is calculated as:

from scipy.spatial import Delaunay
tri = Delaunay(gix)
a_edges = np.sort(tri.simplices[:,[0,1]])
b_edges = np.sort(tri.simplices[:,[0,2]])
c_edges = np.sort(tri.simplices[:,[1,2]])
edges = np.unique(np.concatenate([a_edges,b_edges,c_edges],axis=0),axis=0)
simplices = tri.simplices

And I unwrap the intf by

unw_edgelist = kamui.unwrap_arbitrary(np.angle(intf),edges)

and

unw_mcf = kamui.unwrap_arbitrary(np.angle(intf),edges,tri.simplices)

But the unwrapped results are very different. It seems unw_mcf is not correct.

unw_edgelist:

unw_mcf:

I want to know if I did anything wrong?

BTW, I also find unw_mcf is 10 times faster than unw_edgelist which is difference as @yoyololicon said in #2 .

[yoyolicoris]

Hi @kanglcn, thanks for reporting this!

Hmm, your code looks about right.
Here are a few directions I'll try:

1. Manually check the final cost value after the ILP converged (attached on the res object at https://github.com/yoyololicon/kamui/blob/7fdca9529f07de64b5d2050609c5199d88055a37/kamui/core.py#L103 and https://github.com/yoyololicon/kamui/blob/7fdca9529f07de64b5d2050609c5199d88055a37/kamui/core.py#L148.) I should (and will) add the result object into return values. Although MCF and edgelist methods use the same cost function, their matrix equations differ and could affect the underlying solver. In my HRTF paper, we found that 40% of the time, they have (slightly) different solutions, but not as drastically as you show above.
2. If they have very close cost value, then I'll suspect the cause is on the integration function: https://github.com/yoyololicon/kamui/blob/7fdca9529f07de64b5d2050609c5199d88055a37/kamui/core.py#L15. I ran into a similar issue when using scipy.sparse.csgraph.breadth_first_order and fixed with depth_first_order.

Since you mentioned, MCF is 10 times faster, I suspect (1) more.
I'll look into this when I have time.

[kanglcn]

Hi @yoyololicon , thanks for your kind reply. I follows your views and found that the final cost value for the two approaches are same which implies calculate_m and calculate_k are correct. The integrate function looks also good to me. I have no idea what is wrong there but I found the order of the input point cloud affects the unwrapping result greatly which is abnormal.

I also have a question regarding to the speed. The main optimization is done with scipy.optimize.linprog in this package. Have you ever compare the speed with other solver? I.e., the PuLP for general ILP or the google or-tools for MCF. My test shows a commercial InSAR phase unwrapping program (which implements MCF algorithm) runs ~600 times faster than kamui with the edgelist approach which is very discouraging to me.

[yoyolicoris]

The integrate function looks also good to me. I have no idea what is wrong there but I found the order of the input point cloud affects the unwrapping result greatly which is abnormal.

Ah, then I guess it's a problem of Delaunay.
Try set qhull_options to "QJ" like the following:
https://github.com/yoyololicon/hrtf-ilp/blob/0447f7e09e5cc2ffcb900d625bd9315f670d3213/toa.py#L132

I also have a question regarding to the speed. The main optimization is done with scipy.optimize.linprog in this package. Have you ever compare the speed with other solver? I.e., the PuLP for general ILP or the google or-tools for MCF. My test shows a commercial InSAR phase unwrapping program (which implements MCF algorithm) runs ~600 times faster than kamui with the edgelist approach which is very discouraging to me.

Yeah, that's true. If a dedicated MCF algorithm is used, it's usually faster. Even networks.network_simplex can perform quite well, a few times faster than ILP, as far as I can remember. The problem is that MCF only works on 2D data, while ILP is more general. If efficiency is important, I can add a dedicated MCF-based method for 2D data, or you can reference my old implementation here.