Help with deconvolving photon reads distributed across neighboring pixels I'm building a radiation detector that collects photons in CCD pixels and we can relate the energy of the photon from the intensity of the pixel. To test the detector, I took the following spectrum of Iron 55, which is expected to have a single peak at 6 keV:



The tail at lower values of the energy is interesting because it's not supposed to be there. Taking a closer look at pixels in this range, I see that they are clustered next to each other, with what seems to be energies from a single photon distributed into multiple neighboring pixels.



I am intrigued by the possibility of getting a cleaner spectrum by collecting these pixels together. I could try manually going through and clustering at each pixel, but I thought a more elegant way would be to deconvolve the image using a simple point spread function. In the result of this deconvolution I need 2 things to be true:


the energy of the system should be conserved, so the sum of the pixels (the $L_1$ norm) should be conserved.
there should be no negative pixel intensities, since those are unphysical. 


I did that, using a Gaussian PSF with a half-pixel standard deviation. However, I find that I either need to violate condition 1 of my requirements, or condition 2 of my requirements. Just using the gaussian PSF with $L_1$ norm 1, which conserves the energy, gives me a plot like this:



The white pixels are negative, unphysical values. At this point I'm stumped, but also obsessed with this problem, because it seems like it should have a clean solution. Any advice?
There is some work on resolving point sources from their smoothed signal called super-resolution, along with some elegant methods to solve it using convex-programming. The idea is that if you know that your signal consists of spikes (single points of energy), you can resolve them from the smoothed signal, even when the smoothing causes interference between spikes.

I used this method in my PhD to find the amplitude of 1D spikes from low-pass filtered signals, using the research published here (not mine). I didn't look at 2D signals like yours, but I recall that some work has been done on it, so you might be able to find some code.

For your method, I'm guessing the negative pixels occur because the smearing isn't exactly modelled by your PSF

What about something simple:


Assume that the smear occurs within a 3 x 3 radius.
Convolve the image with a 3 x 3 filter of all ones. This will add all the values in that area.
Find the points of local maxima in a 3 x 3 radius in the original image.
Replace the value of those points with those of the convolved image and zero the rest.