Phase unwrapping: Branch cut algorithm #680

josteinbf · 2013-08-08T19:16:41Z

My other PR on phase unwrapping, #644, is stalled (copyright issues), but I've had another implementation lying around for a while which is perfectly safe to include.

This PR implements the branch cut algorithm. Anyone wanting to review this in detail is encouraged to have a look at Goldstein's classic paper on phase unwrapping[1] before diving in; I've implemented the same algorithm with some slight modifications, mostly to allow treating masked arrays properly.

Functionality, tests, docs and examples should be about complete. I lack a test case for masked array input---ideas for good test cases welcome!

[1] R. M. Goldstein, H. A. Zebker, C. L. Werner, "Satellite radar
interferometry: Two-dimensional phase unwrapping", Radio Science 23
(1988) 4, pp 713--720.

josteinbf · 2013-09-07T05:34:13Z

@ahojnnes @JDWarner Any chance of any of you reviewing this?

ahojnnes · 2013-09-14T08:35:10Z

@josteinbf Sorry for the delay. I am not sure if we really need a separate queue implementation. Is it performance critical (if not, we could use PyList)? Otherwise the queue implementation should be moved to skimage/_shared.

Apart from that the implementation looks very good!

josteinbf · 2013-09-14T11:33:24Z

@ahojnnes Thanks for looking into this, and thanks for your complements on the implementation. Performance-wise I guess it should be good enough, but I'm worried about readability. How does readability look to someone who did not code it? Do you feel you need to struggle to understand what's going on?

With regards to the queue, it is used rather intensively, but I have not really researched the options here much---I just figured any basic BSD-licensed C implementation would do. The average number of items in the queue is probably proportional to the number of pixels in the image (the proportionality constant is likely highly dependent on the input, but it could in bad cases get close to one), and items are enqueued and dequeued at a high frequency (all pixels enter the queue at least twice). By PyList, do you mean a python list, just accessed from C? I have no idea how efficient that is, never used something like that before. In summary, I'd very much appreciate advise on the queue implementation :) As this implementation only deals with 2D images the memory overhead of the queue is not critical, but using more memory for the queue could hurt performance through more fetches/stores (the number of flops per queue item is small).

ahojnnes · 2013-09-15T17:58:53Z

skimage/exposure/unwrap.py

+
+
+def find_branch_cuts(residues):
+    '''Connect residues with branch cuts such that the length of the cuts


One-liner as a short description.

ahojnnes · 2013-09-15T18:17:03Z

@josteinbf Let's keep the current queue implementation, but move it to skimage/_shared. Maybe we can also wrap it into a nice class, e.g. http://docs.cython.org/src/tutorial/clibraries.html ?

stefanv · 2013-09-15T21:23:53Z

Here's one description of the allocation pattern: http://www.laurentluce.com/posts/python-list-implementation/

ahojnnes · 2013-09-16T12:34:55Z

skimage/exposure/_unwrap_naive.pyx

+from libc.math cimport M_PI
+
+
+def unwrap_naive_1d(double[::1] image, double[::1] unwrapped_image):


double[::1] => double[:]

Can you please explain the difference? I thought double[::1] would assume a C-ordered array (not patched together by memory chunks scattered in different places, but one continuous allocation), and allow cython to deal with this more efficiently?

Sorry, you are right, of course. Disregard my suggestion ;-)

josteinbf · 2013-09-18T21:46:43Z

@ahojnnes Thanks for your comments. This week is too busy for me to respond properly to this, but I hope to get to it next week.

josteinbf · 2013-09-27T18:58:54Z

@ahojnnes The queue OO wrapper is done now and in shared. It does indeed look a lot better, good suggestion. One question: I have assumed that cython will deallocate the queue for me when the queue variable goes out of scope, can you confirm that this is correct? If I'm wrong about that, the code will leak memory.

josteinbf · 2013-09-29T07:48:04Z

@stefanv Thanks for the link on allocation patterns. The contents convince me that PyList is not an option for this code, since pop(0) and insert(0) are common operations in the code, both of which are O(n). With the separate queue implementation, these are both O(1).

josteinbf · 2013-09-29T09:23:57Z

@ahojnnes I think I have addressed all of your comments now. Where I have not answered to your comment, I have just implemented your suggestion. There are a few questions for you that I hope you can answer above.

josteinbf · 2013-11-24T09:08:06Z

I'm in the process of rebasing this, which will require a fair few changes, since there is a substantial overlap with #644, which was just merged. Some decisions need to be made; I plan to do the following:

The public API will be a single function with a method kwarg choosing between the Goldstein implementation (this PR) and the more sophisticated method implemented by Phase unwrapping #644: unwrap_phase(..., method=None).
find_residues will be a public function, as it is useful for diagnosing problems when unwrapping.
The example will demonstrate the reliability method (the default), so it will not change much from Phase unwrapping #644.

Comments welcome.

ahojnnes · 2013-11-24T09:40:59Z

@josteinbf This sounds like a perfect plan to me! find_residues should also definitely be public.

Maybe you can combine both examples into one script comparing the methods side-by-side?

josteinbf · 2013-11-24T10:24:19Z

I'll have a look at what the differences between the methods really are. However, I fear that since we have only relatively unrealistic example data, we will not be able to demonstrate how the methods are different. If we want to do that, we would probably have to include an extra image. Would that be possible? It would of course increase the size of the distribution.

ahojnnes · 2013-11-24T10:59:51Z

What kind of images do you have and what size are they?

The C code comes from the C-algorithms library: http://c-algorithms.sourceforge.net/ We provide the cython wrapper class developed as part of a Cython tutorial: http://docs.cython.org/src/tutorial/clibraries.html

josteinbf · 2013-11-24T16:05:48Z

@ahojnnes Rebased.

Is there an easy way I can view a test coverage report from this fancy coverall thing? With the modified API, I'm certain there are some untested lines lying around relating to the method kwarg.

josteinbf · 2013-11-24T16:07:23Z

@ahojnnes Regarding comparisons in the example: I think we should do that in a separate PR (if we should do it at all---the papers cited will certainly have information on this). This one's big enough already.

josteinbf · 2013-11-24T16:27:49Z

A single doctest is failing, and I'm having trouble seeing what the problem is:

======================================================================
FAIL: Doctest: skimage.exposure.unwrap.find_phase_residues
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/usr/lib/python2.7/doctest.py", line 2201, in runTest
    raise self.failureException(self.format_failure(new.getvalue()))
AssertionError: Failed doctest test for skimage.exposure.unwrap.find_phase_residues
  File "/home/jostein/Programming/scikit-image/skimage/exposure/unwrap.py", line 308, in find_phase_residues

----------------------------------------------------------------------
File "/home/jostein/Programming/scikit-image/skimage/exposure/unwrap.py", line 349, in skimage.exposure.unwrap.find_phase_residues
Failed example:
    find_phase_residues(image)   # "--" indicates a masked element
Expected:
    masked_array(data =
     [[0 0 0 --]
     [0 1 0 --]
     [0 0 0 --]
     [-- -- -- --]],
                 mask =
     [[False False False  True]
     [False False False  True]
     [False False False  True]
     [ True  True  True  True]],
           fill_value = 999999)
Got:
    masked_array(data =
     [[0 0 0 --]
     [0 1 0 --]
     [0 0 0 --]
     [-- -- -- --]],
                 mask =
     [[False False False  True]
     [False False False  True]
     [False False False  True]
     [ True  True  True  True]],
           fill_value = 999999)
    <BLANKLINE>


----------------------------------------------------------------------

I can only see a single difference, the <BLANKLINE> at the end. Is that really the problem? If so, how to fix it? I've already tried adding extra blank lines in the docstring.

josteinbf · 2013-11-24T16:33:15Z

A bit of googling turned up the solution to the blank line issue. I expect tests to pass now.

coveralls · 2013-11-24T16:49:59Z

Coverage decreased (-0.79%) when pulling 5ab2ad4 on josteinbf:phase-unwrap-goldstein into caf183a on scikit-image:master.

josteinbf · 2013-11-24T18:46:03Z

c86aae2 should have failing tests. These tests may mean that the branch cuts implementation has a bug related to handling masked regions---I don't know yet. Until I know, please don't merge.

stefanv · 2014-11-10T12:43:58Z

This PR's birthday is almost here!

aaron-parsons · 2018-10-30T13:39:50Z

Hi.

Myself, @darrenbatey and @pierrethibault are currently looking at phase unwrapping, and possibly extending this functionality in scikit-image.

What's the state of this PR? If we were to update it so that the files are in the right place, what's left to do to get the PR through?

Aaron

stefanv · 2018-10-30T19:25:13Z

@aaron-parsons That's excellent news; a thorough review at this point would be most welcome and helpful, especially given that the PR has been dormant since 2014. Thank you!

ahojnnes reviewed Sep 15, 2013
View reviewed changes

ahojnnes reviewed Sep 16, 2013
View reviewed changes

josteinbf mentioned this pull request Nov 22, 2013

Hausdorff distance #742

Closed

josteinbf added 12 commits November 24, 2013 16:45

unwrap: Rename unwrap_naive extension.

b22aeca

unwrap: Rename branch cuts extension.

4744637

Add a C queue implementation to _shared.

88d0cd9

The C code comes from the C-algorithms library: http://c-algorithms.sourceforge.net/ We provide the cython wrapper class developed as part of a Cython tutorial: http://docs.cython.org/src/tutorial/clibraries.html

unwrap: Use an OO wrapper for the Calg queue.

acea43c

Add skimage._shared.queue to bento.info.

a0b3daa

unwrap: Shorten first sentence of docstrings.

cf0fba5

unwrap: Prefer Py_ssize_t for indices.

faea4fc

unwrap: Fix indentation.

5fed5a9

unwrap: Turn off cython boundschecking.

07e52ad

unwrap: Improve comments in cython code.

f2df5b9

unwrap: Clarify difference between pixel and intersection coordinates.

daf57f6

unwrap: Do not start cython function names with underscores.

7da12bb

unwrap: Fix doctest.

6deab11

unwrap: Fix doctest.

5ab2ad4

unwrap: Add test for masked arrays using branch cuts.

c86aae2

soupault added the type: new feature label Nov 7, 2018

Base automatically changed from master to main February 18, 2021 18:22

mkcor mentioned this pull request Feb 22, 2021

2021's calendar of community management #5169

Closed

grlee77 added 🙏 Feature request and removed type: new feature labels Feb 22, 2022

mkcor mentioned this pull request Feb 28, 2022

2022's calendar of community management #6165

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Phase unwrapping: Branch cut algorithm #680

Phase unwrapping: Branch cut algorithm #680

josteinbf commented Aug 8, 2013

josteinbf commented Sep 7, 2013

ahojnnes commented Sep 14, 2013

josteinbf commented Sep 14, 2013

ahojnnes Sep 15, 2013

ahojnnes commented Sep 15, 2013

stefanv commented Sep 15, 2013

ahojnnes Sep 16, 2013

josteinbf Sep 29, 2013

ahojnnes Sep 29, 2013

josteinbf commented Sep 18, 2013

josteinbf commented Sep 27, 2013

josteinbf commented Sep 29, 2013

josteinbf commented Sep 29, 2013

josteinbf commented Nov 24, 2013

ahojnnes commented Nov 24, 2013

josteinbf commented Nov 24, 2013

ahojnnes commented Nov 24, 2013

josteinbf commented Nov 24, 2013

josteinbf commented Nov 24, 2013

josteinbf commented Nov 24, 2013

josteinbf commented Nov 24, 2013

coveralls commented Nov 24, 2013

josteinbf commented Nov 24, 2013

stefanv commented Nov 10, 2014

aaron-parsons commented Oct 30, 2018

stefanv commented Oct 30, 2018



		def find_branch_cuts(residues):
		'''Connect residues with branch cuts such that the length of the cuts

		from libc.math cimport M_PI


		def unwrap_naive_1d(double[::1] image, double[::1] unwrapped_image):

Phase unwrapping: Branch cut algorithm #680

Are you sure you want to change the base?

Phase unwrapping: Branch cut algorithm #680

Conversation

josteinbf commented Aug 8, 2013

josteinbf commented Sep 7, 2013

ahojnnes commented Sep 14, 2013

josteinbf commented Sep 14, 2013

ahojnnes Sep 15, 2013

Choose a reason for hiding this comment

ahojnnes commented Sep 15, 2013

stefanv commented Sep 15, 2013

ahojnnes Sep 16, 2013

Choose a reason for hiding this comment

josteinbf Sep 29, 2013

Choose a reason for hiding this comment

ahojnnes Sep 29, 2013

Choose a reason for hiding this comment

josteinbf commented Sep 18, 2013

josteinbf commented Sep 27, 2013

josteinbf commented Sep 29, 2013

josteinbf commented Sep 29, 2013

josteinbf commented Nov 24, 2013

ahojnnes commented Nov 24, 2013

josteinbf commented Nov 24, 2013

ahojnnes commented Nov 24, 2013

josteinbf commented Nov 24, 2013

josteinbf commented Nov 24, 2013

josteinbf commented Nov 24, 2013

josteinbf commented Nov 24, 2013

coveralls commented Nov 24, 2013

josteinbf commented Nov 24, 2013

stefanv commented Nov 10, 2014

aaron-parsons commented Oct 30, 2018

stefanv commented Oct 30, 2018