Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Expand Gabor Kernel Generation to nD #2739

Open
wants to merge 53 commits into
base: main
Choose a base branch
from
Open

Conversation

kne42
Copy link
Contributor

@kne42 kne42 commented Jul 28, 2017

Description

Expands the Gabor filter to work on images of any dimension.

Checklist

[For detailed information on these and other aspects see scikit-image contribution guidelines]

References

For reviewers

(Don't remove the checklist below.)

  • Check that the PR title is short, concise, and will make sense 1 year
    later.
  • Check that new functions are imported in corresponding __init__.py.
  • Check that new features, API changes, and deprecations are mentioned in
    doc/release/release_dev.rst.

@pep8speaks
Copy link

pep8speaks commented Jul 28, 2017

Hello @kne42! Thanks for updating the PR.

Line 181:80: E501 line too long (89 > 79 characters)
Line 182:80: E501 line too long (89 > 79 characters)
Line 183:80: E501 line too long (89 > 79 characters)
Line 184:80: E501 line too long (89 > 79 characters)
Line 185:80: E501 line too long (89 > 79 characters)

Line 59:47: E127 continuation line over-indented for visual indent
Line 60:47: E127 continuation line over-indented for visual indent
Line 61:47: E127 continuation line over-indented for visual indent
Line 62:47: E127 continuation line over-indented for visual indent

Comment last updated on May 30, 2018 at 21:40 Hours UTC

@kne42 kne42 changed the title [WIP][Enhancement] N-Dimensional Gabor Filter ENH: N-Dimensional Gabor Filter Jul 28, 2017
@kne42 kne42 changed the title ENH: N-Dimensional Gabor Filter (4/4) N-Dimensional Gabor Filter Aug 6, 2017
@kne42 kne42 force-pushed the gabor-nD branch 2 times, most recently from c62318b to 5fbfe7a Compare April 19, 2018 19:06
@kne42 kne42 changed the title (4/4) N-Dimensional Gabor Filter n-Dimensional Gabor Filter Apr 19, 2018
@kne42 kne42 force-pushed the gabor-nD branch 2 times, most recently from f4dac92 to 923bcb2 Compare April 21, 2018 22:55
@kne42 kne42 changed the title n-Dimensional Gabor Filter n-Dimensional Gabor Filter (WIP) Apr 25, 2018

# handle deprecation
message = ('Using deprecated, 2D-only interface to'
'gabor_kernel. This interface will be'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You need spaces at the end of each sub-string, since concatenation within parentheses won't add spaces. =)

message = ('Using deprecated, 2D-only interface to'
'gabor_kernel. This interface will be'
'removed in scikit-image 0.16. Use'
'gabor_kernel(frequency, rotation=theta,'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You say rotation=theta here but the signature does not have a rotation keyword. fwiw I think theta as a kwarg can be maintained, silently adding support for multidimensional theta.

if sigma_y is None:
sigma_y = _sigma_prefactor(bandwidth) / frequency
def __new__(cls, frequency, thetas=None, bandwidth=1, sigma=None,
sigma_y=None, n_stds=3, offset=None, ndim=2, **kwargs):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure that the complexity of a class is necessary. Why did you pick this design instead of extending the existing function?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jni I chose this design because it encourages reusing the existing gabor kernel. As of right now, the filters.gabor function creates a new gabor kernel every time you use it. This becomes very expensive if I, say, want to apply the same gabor filter to 1000 images.

While there is already a method of eliminating this expense by creating a gabor kernel and applying it to the images yourself, I would posit that most users would not consider this as they would immediately gravitate towards the gabor function, which is not obviously expensive from the documentation when applied many times. Furthermore, if one wanted to reproduce what the gabor function did, they would have to look at the source code and figure out exactly how we were able to apply the filter to an image.

Although I generally prefer functional design to object-oriented design, I think that this choice is justified for the following reasons:

  • NumPy is already object-oriented. If we're representing the kernel as an np.ndarray anyways, it would be most useful to extend this interface to apply it as well.
  • This tags along the filter function to the kernel instead of having to use one function to make the kernel and another to apply it (this would be the functional alternative). This is also a lot more intuitive.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kne42 I think the simplest option is to add to filters.gabor a kernel parameter that is the (bare) NumPy array, and avoid recomputing the kernel if this argument is given.

The NumPy community has moved away from inheriting from NumPy arrays because it raises a whole bunch of problems. An example: suppose someone comes up with a fancy new "Gabor++" method that works by creating a Gabor kernel, modifying in some weird way, then applying it. Any call to np.copy(kernel) will strip away the class and return a bare array. Then you've lost your way of applying the kernel.

The other issue is the "surprise" factor: this is a pattern that does not exist anywhere in skimage, or NumPy/SciPy for that matter. So now users have to wrap their heads around a whole new API for Gabor kernels. Or they run g = filters.gabor_kernel(freq) and then type type(g) and get a strange class that they've never seen.

@jni
Copy link
Member

jni commented Apr 26, 2018

@kne42 this is looking nice! Please ping when it's ready for a full review.

Speaking of pinging, did you receive my email? I used the address on your GH profile.

@kne42
Copy link
Contributor Author

kne42 commented Apr 30, 2018 via email

@kne42
Copy link
Contributor Author

kne42 commented May 16, 2018

@jni @emmanuelle This PR is now ready for its penultimate review as I start cleaning up documentation/implementation, working out kinks in the functions, and expanding unit test coverage.

Some points of discussion:

  • What to do with all of these helper functions. Which, if any, do you think should be made public? If we make them public, should we introduce this refactor in a separate PR for cleanliness or just make the changes in this one?
  • Deprecation handling & backwards-compatibility have resulted in some less-than-ideal design choices. How do you think we should address this?
  • Thoughts on general design/clarity?

@kne42 kne42 changed the title n-Dimensional Gabor Filter (WIP) Expand Gabor Kernel Generation to nD May 16, 2018
@codecov-io
Copy link

codecov-io commented May 21, 2018

Codecov Report

Merging #2739 into master will increase coverage by 0.05%.
The diff coverage is 96.25%.

Impacted file tree graph

@@            Coverage Diff            @@
##           master   #2739      +/-   ##
=========================================
+ Coverage   86.05%   86.1%   +0.05%     
=========================================
  Files         338     338              
  Lines       27364   27504     +140     
=========================================
+ Hits        23547   23683     +136     
- Misses       3817    3821       +4
Impacted Files Coverage Δ
skimage/filters/tests/test_gabor.py 98.09% <100%> (+2.01%) ⬆️
skimage/filters/_gabor.py 94.73% <94.28%> (-5.27%) ⬇️
skimage/draw/_random_shapes.py 97.89% <0%> (+2.1%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ef78faf...abea636. Read the comment docs.


norm = np.linalg.norm(u)

if not np.isclose(norm, 1):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this check necessary?



def _compute_projection_matrix(axis, indices=None):
"""Generates a matrix that projects an axis onto the 0th coordinate axis.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the 0th coordinate axis?


References
----------
.. [1] Bart M. ter Haar Romeny. Front-End Vision and Multi-Scale Image
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a collection of articles; could you add the name (and perhaps page number) of the relevant one?

>>> Y = np.asarray([0.5, 0.5])

>>> M = _compute_rotation_matrix(X, Y)
>>> Z = np.matmul(M, X)
Copy link
Member

@stefanv stefanv May 28, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Z = M @ X

Parameters
----------
image : non-complex array
Linear space to map the filter to.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How is the image a linear space here?

Standard deviation for Gaussian kernel. The standard deviations of the
Gaussian filter are given for each axis as a sequence, or as a single
number, in which case it is equal for all axes.
ndim : int, optional
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm wondering here: should we allow center and sigma to be scalars? If not, we can derive ndim from center and sigma.

offset : float, optional
Phase offset of harmonic function in radians.
axes : int or sequence of int, optional
Ordering of axes that defines the plane of rotation. Ordering not
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure what this parameter does?

if sigma_y is None:
sigma_y = _sigma_prefactor(bandwidth) / frequency
# handle deprecation
message = ('Using deprecated, 2D-only interface to gabor_kernel. '
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there no way to support the old style calling convention still?

m = np.asarray(np.meshgrid(*[range(-c, c + 1) for c in spatial_size],
indexing='ij'))

rotm = np.matmul(m.T, rot.T).T
Copy link
Member

@stefanv stefanv May 28, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rot @ m?

R = np.eye(ndim) # Initial rotation matrix = Identity matrix

# Loop to create matrices of stages
for step in np.round(2 ** np.arange(np.log2(ndim))).astype(int):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Clarify what happens here: why the log2, the 2**, etc.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest linking to

https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions

and pointing out that this function will compute the correct matrices for rotation along individual axes in the same way that is described for 3D in that page.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add the same links to the general generate_rotation_matrix function

if n + step >= len(w):
break

r2 = x[w[n]] * x[w[n]] + x[w[n + step]] * x[w[n + step]]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

r2 = x[w[n]] ** 2 + x[w[n + step]] ** 2

Ah, but then we take the sqrt, so how about np.hypot?

r2 = x[w[n]] * x[w[n]] + x[w[n + step]] * x[w[n + step]]
if r2 > 0:
r = np.sqrt(r2)
pcos = x[w[n]] / r # Calculation of coefficients
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, x[w[n]] and x[w[n + step]]. Since these indices are re-used, should we name them?

dst : (N, ) array
Vector of desired direction.
use_homogeneous_coords : bool, optional
If the input vectors should be treated as homoegeneous coordinates.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

homogeneous

r : float
Radial coordinate.
thetas : (N, ) array
Quasipolar angles.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some explanation of where these angles lie would be helpful (given the "quasi" nature you explained to me).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants