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#L58 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#L34 On Wed, Nov 22, 2017, at 00:16, Jerome Kieffer wrote: > On Wed, 22 Nov 2017 12:40:50 +1100 > Juan Nunez-Iglesias <jni.s...@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 _______________________________________________ scikit-image mailing list scikit-image@python.org https://mail.python.org/mailman/listinfo/scikit-image
