Dear all, I have an image which is "blurred" by a kernel which depends on the position on the image. This blurring can be expressed as a sparse matrix (A) multiplication where only the neighboring pixels have non-null contribution. Ax = b where in addition Aij>=0 x >= 0 #non negativity constrain. b >= 0 # measured signal Does anyone have some hints on where to start looking at ? Thanks for your help -- Jérôme Kieffer
Hi Jérôme, Can you explain your problem more? You know A and x and want to find b? Is this an exact solution, or is Ax = b + err? SciPy’s sparse.linalg module is where you’ll find most of your answers, I think… If you want to *build* A from some description, you might find our homography example in Elegant SciPy useful: https://github.com/elegant-scipy/elegant-scipy/blob/master/markdown/ch5.mark... Juan. On 22 Nov 2017, 1:58 AM +1100, Jerome Kieffer <google@terre-adelie.org>, wrote:
Dear all,
I have an image which is "blurred" by a kernel which depends on the position on the image. This blurring can be expressed as a sparse matrix (A) multiplication where only the neighboring pixels have non-null contribution.
Ax = b
where in addition Aij>=0 x >= 0 #non negativity constrain. b >= 0 # measured signal
Does anyone have some hints on where to start looking at ? Thanks for your help
-- Jérôme Kieffer _______________________________________________ scikit-image mailing list scikit-image@python.org https://mail.python.org/mailman/listinfo/scikit-image
On Wed, 22 Nov 2017 12:40:50 +1100 Juan Nunez-Iglesias <jni.soma@gmail.com> wrote:
Hi Jérôme,
Can you explain your problem more? You know A and x and want to find b? Is this an exact solution, or is Ax = b + err? SciPy’s sparse.linalg module is where you’ll find most of your answers, I think… If you want to *build* A from some description, you might find our homography example in Elegant SciPy useful:
Hi Juan, Thanks for your help. Nice documentation on scipy.sparse I wish I had it a couple of days ago. Indeed I have spent a week in building the matrix A using ray-tracing and now I believe it is almost correct now. The idea is to consider some image sensor (in 2D) which absorb the photon (X-ray) in volume (instead of the surface). So each pixel is considered as a voxel and the sensor is a 2D array of voxels. After this raytracing step, I know how much of a photon arriving in one pixel contributes to the neighboring pixels, this is my matrix A, (in CSR format). "b" is of course the image I read from the detector, with all element positive and I expect "x" to have all element positive as well. I tried to search in scipy.sparse.linalg the method which would allow me to retrieve "x" from "b". For now, the best I found is : res = linalg.lsmr(A, b, damp=damp, x0=b) When the damping factor is null or too small, I notice many wiggles (with negative regions) near peaks. For now I try to adjust this damp factor so that the minimum of x is positive but I am not confident on the method. -- Jérôme Kieffer tel +33 476 882 445
Hi Jérôme, Sounds great! I’ll admit that I’m not confident in this area either, but my reading of the documentation suggests that this is the right approach: the damping controls the square norm of x. By keeping it small (with large damping), you force the elements to be non-negative. I hope the final picture looks good now! =) If you get a nice result, I suggest you write a blog post about it, with pictures. It sounds like a very cool use of SciPy, and would be a valuable addition to writeups about it! Juan. On 22 Nov 2017, 7:16 PM +1100, Jerome Kieffer <google@terre-adelie.org>, wrote:
On Wed, 22 Nov 2017 12:40:50 +1100 Juan Nunez-Iglesias <jni.soma@gmail.com> wrote:
Hi Jérôme,
Can you explain your problem more? You know A and x and want to find b? Is this an exact solution, or is Ax = b + err? SciPy’s sparse.linalg module is where you’ll find most of your answers, I think… If you want to *build* A from some description, you might find our homography example in Elegant SciPy useful:
Hi Juan,
Thanks for your help. Nice documentation on scipy.sparse I wish I had it a couple of days ago.
Indeed I have spent a week in building the matrix A using ray-tracing and now I believe it is almost correct now.
The idea is to consider some image sensor (in 2D) which absorb the photon (X-ray) in volume (instead of the surface). So each pixel is considered as a voxel and the sensor is a 2D array of voxels. After this raytracing step, I know how much of a photon arriving in one pixel contributes to the neighboring pixels, this is my matrix A, (in CSR format).
"b" is of course the image I read from the detector, with all element positive and I expect "x" to have all element positive as well.
I tried to search in scipy.sparse.linalg the method which would allow me to retrieve "x" from "b". For now, the best I found is :
res = linalg.lsmr(A, b, damp=damp, x0=b)
When the damping factor is null or too small, I notice many wiggles (with negative regions) near peaks. For now I try to adjust this damp factor so that the minimum of x is positive but I am not confident on the method.
-- Jérôme Kieffer tel +33 476 882 445 _______________________________________________ scikit-image mailing list scikit-image@python.org https://mail.python.org/mailman/listinfo/scikit-image
Hi Juan, On Thu, 23 Nov 2017 12:33:21 +1100 Juan Nunez-Iglesias <jni.soma@gmail.com> wrote:
Sounds great! I’ll admit that I’m not confident in this area either, but my reading of the documentation suggests that this is the right approach: the damping controls the square norm of x. By keeping it small (with large damping), you force the elements to be non-negative.
The approach looks OK, I am just a bit "unconfident" about tweeking a factor to get the image I am expecting. Not very scientific when the resulting image has to used as input for quantitative analysis.
I hope the final picture looks good now! =) If you get a nice result, I suggest you write a blog post about it, with pictures. It sounds like a very cool use of SciPy, and would be a valuable addition to writeups about it!
This will be part of the documentation of pyFAI by the end of the year: http://pyfai.readthedocs.io/en/latest/usage/tutorial/index.html To Stefan, I am sorry I forgot about your PhD work. I will check again your work, last time it was a different topic :) Cheers, Jerome
On 22 November 2017 23:08:17 Jérôme Kieffer <google@terre-adelie.org> wrote:
On Thu, 23 Nov 2017 12:33:21 +1100 Juan Nunez-Iglesias <jni.soma@gmail.com> wrote:
Sounds great! I’ll admit that I’m not confident in this area either, but my reading of the documentation suggests that this is the right approach: the damping controls the square norm of x. By keeping it small (with large damping), you force the elements to be non-negative.
The approach looks OK, I am just a bit "unconfident" about tweeking a factor to get the image I am expecting. Not very scientific when the resulting image has to used as input for quantitative analysis.
There are various ways to do this, but you have to dampen/penalize the solution: oscillatory results often achieve the same norm as smooth results, and the optimization algorithm otherwise has no way to choose which is best. Instead of applying explicit regularization, you can also terminate conjugate gradients early, for example. Best regards Stéfan
Hi Jerome The problem you describe sounds very similar to the one solved for super-resolution imaging. My PhD thesis talks quite a bit about the methods utilized to solve it: http://mentat.za.net/phd_dissertation.html The code is here, but---written by some PhD student in 2009 ;) https://github.com/stefanv/supreme/blob/master/supreme/resolve/iterative.py#... One of the tricks we used in SR imaging was to find a good smooth estimate for the resulting image, and then used the optimization to search for updates of that "prior". You can then penalize the updates not to be too large, and thereby prevent the oscillations you mention. Best regards Stéfan P.S. Juan, w.r.t. the Elegant SciPy sparse example: https://github.com/stefanv/supreme/blob/master/supreme/resolve/operators.pyx... On Wed, Nov 22, 2017, at 00:16, Jerome Kieffer wrote:
On Wed, 22 Nov 2017 12:40:50 +1100 Juan Nunez-Iglesias <jni.soma@gmail.com> wrote:
Hi Jérôme,
Can you explain your problem more? You know A and x and want to find b? Is this an exact solution, or is Ax = b + err? SciPy’s sparse.linalg module is where you’ll find most of your answers, I think… If you want to *build* A from some description, you might find our homography example in Elegant SciPy useful:
Hi Juan,
Thanks for your help. Nice documentation on scipy.sparse I wish I had it a couple of days ago.
Indeed I have spent a week in building the matrix A using ray-tracing and now I believe it is almost correct now.
The idea is to consider some image sensor (in 2D) which absorb the photon (X-ray) in volume (instead of the surface). So each pixel is considered as a voxel and the sensor is a 2D array of voxels. After this raytracing step, I know how much of a photon arriving in one pixel contributes to the neighboring pixels, this is my matrix A, (in CSR format).
"b" is of course the image I read from the detector, with all element positive and I expect "x" to have all element positive as well.
I tried to search in scipy.sparse.linalg the method which would allow me to retrieve "x" from "b". For now, the best I found is :
res = linalg.lsmr(A, b, damp=damp, x0=b)
When the damping factor is null or too small, I notice many wiggles (with negative regions) near peaks. For now I try to adjust this damp factor so that the minimum of x is positive but I am not confident on the method.
-- Jérôme Kieffer tel +33 476 882 445 _______________________________________________ scikit-image mailing list scikit-image@python.org https://mail.python.org/mailman/listinfo/scikit-image
participants (4)
-
Jerome Kieffer -
Juan Nunez-Iglesias -
Jérôme Kieffer
-
Stefan van der Walt