PWD based operators#706
Conversation
cako
left a comment
There was a problem hiding this comment.
@FB-I thanks a lot for this PR! This is really awesome! I found the PR very easy to read and very much in line with the rest of the operators in the library. Overall great stuff.
I have left a few comments and code suggestions, mostly bringing it up to speed with the latest PyLops features (a.k.a, reshaped, automatic shape, etc.). Otherwise just a few nits here and there.
To expand a bit more on the topic, PyLops since version 2 automatically reshapes arrays based on its inputs. So you don't need to pass x.ravel() to the operator or reshape it at the output. This is done with a bit of decorator magic (see reshaped) as well as proper super class initialization. Basically if you decorate your _matvec and _rmatvec with @Reshaped and then ensure that you call super with dims and dimsd arguments, you are all good to go. See the Sum operator for an example.
I haven't had time to look at the results yet, but I have noticed that the smoothing looks a but weird with the PWD slopes (but works with the structural-tensor slopes). Could point to some incorrect implementation but I haven't had time to dig:

Co-authored-by: Carlos da Costa <cako@users.noreply.github.com>
Co-authored-by: Carlos da Costa <cako@users.noreply.github.com>
Co-authored-by: Carlos da Costa <cako@users.noreply.github.com>
Co-authored-by: Carlos da Costa <cako@users.noreply.github.com>
Co-authored-by: Carlos da Costa <cako@users.noreply.github.com>
Co-authored-by: Carlos da Costa <cako@users.noreply.github.com>
|
@FB-I indeed, great work! I agree with Carlos that the code is well written and well organized and that these new features are going to be very useful! I have so far only looked at the example in details and cleaned up a bit (trying to make it as consistent as possible to the rest of our example). I also add a small 3D example at the end... I am not aware of any real 3D implementation of PWD, usually what one does is cascaded 2D versions (yes, one could at least solve the estimation process for all slices together and impose some smoothness in the extra dimension, but I don't think that even Madagascar does that 😉 ) and this is the approach I took here. I am now going through the other files, I will push some minor re-organization of functions in files tomorrow, then I think this is very close to be ready to be merged 🚀 |
…ble, it promotes types to float64. Also, using `float` as a function also promotes to float64. If we want to keep all computations in float32, we need to cast to that type.
|
Looking pretty good to me! I'll let you handle the final changes that @mrava87 suggested and let him approve it. On my end I made a few commits fixing some typing, improving conciseness and adjusting some Numba code. I also change the |
| inputfile = "../testdata/sigmoid.npz" | ||
|
|
||
| sigmoid = np.load(inputfile)["sigmoid"].T | ||
| sigmoid = 1e3 * np.load(inputfile)["sigmoid"].T |
There was a problem hiding this comment.
Oh, I normally do this to sigmoid to bring it to order of magnitude ~1. Totally optional.
|
I think that Madagascar is actually doing this "one could at least solve the estimation process for all slices together and impose some smoothness in the extra dimension, but I don't think that even Madagascar does that", all slices are updated in one linear solve and they use a 3d smoother for regularization to link the 3rd dim. But I can look into that if it is needed and later update the 3d implementation . Also in @mrava87 new example numpy.concat() does not seem to work unless changed in numpy.concatenate() |
Sounds good! No need to do it now 😄
|
|
@FB-I as you will see I did a bit of reshuffling of codes (mostly the internal ones) to try to keep as much consistency as possible with the rest of PyLops codebase. Other than that, I am also happy now with how things look. If you are ok, I am ready t merge? 🚀 |
|
@mrava87 Sure, thank you. It is ok to merge from my side. |
Proposing addition to PyLops of PW based Structural Smoothing (inherently based on Spraying) and PWD local slope estimation. Both operators work in 2D. Simple example provided.