Vertical streak artifacts

[rijobro]

I've noticed vertical streak artifacts on my MLEM reconstructions (also those in STIR using the NP projectors). I'm using the example Zenodo data. Am I doing something wrong or is this a bug?

I've created a small jupyter notebook (attached, saved as .txt) that demonstrates what I mean. Very similar to @casperdcl's demo, but with no cylindrical truncation or PSF. I've used norm and randoms, but no scatter or attenuation.

The important chunk of code is below, and the orthogonal views after 50 iterations is after that.

iters = 50
prompts = m['psino'].astype(np.float32)
additive = rands
multiplicative = norm

def make_non_zero(A):
    """Make an array non-zero"""
    A[A <= 0] = np.min(A[A > 0]) * 0.001
    return A

def fwd(im):
    """Forward project"""
    out = nipet.frwd_prj(im, mMRpars)
    if multiplicative is not None:
        out *= multiplicative
    if additive is not None:
        out += additive
    out = make_non_zero(out)
    return out

def bck(sino):
    """Back project"""
    sino_to_proj = np.copy(sino)
    if multiplicative is not None:
        sino_to_proj *= multiplicative
    out = nipet.back_prj(sino_to_proj, mMRpars)
    return out

bck_1s = bck(np.ones_like(prompts))
bck_1s = make_non_zero(bck_1s)
inv_bck_1s = 1 / bck_1s

recon_mlem = [np.ones_like(bck_1s)]
for k in trange(iters, desc="MLEM"):
    fprj = fwd(recon_mlem[-1])
    recon_mlem.append(recon_mlem[-1] * inv_bck_1s * bck(prompts / fprj))

Lastly, there's also some horizontal artifacts, that can be seen on this zoomed image of the chin:

[rijobro]

Hi, a little bump on this. Has anyone had a chance to have a look?
I wonder if this potential bug exists in the forward and back projector, but doesn't exist in the GPU OSEM implementation?

[pjmark]

Yes, I come across them particularly when using NAC recons. These artefacts seem to be caused by the highly defined edges of the bed and some other factors. When the mu-maps are well aligned this problem should not appear straight to the beholder. I noticed it even on Siemens reconstructions to some extent as well. Cheers, P

…

________________________________ From: Richard Brown <notifications@github.com> Sent: 29 June 2020 17:49 To: NiftyPET/NIPET <NIPET@noreply.github.com> Cc: Subscribed <subscribed@noreply.github.com> Subject: Re: [NiftyPET/NIPET] Vertical streak artifacts in MLEM (#20) Hi, a little bump on this. Has anyone had a chance to have a look? I wonder if this potential bug exists in the forward and back projector, but doesn't exist in the GPU OSEM implementation? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#20 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABPEVHTLGO4LOZ4MZYBVWN3RZDA2ZANCNFSM4NYXAXUA>.

[casperdcl]

@rijobro / @pjmark I can confirm that I can reproduce these images. However if I add attenuation correction (using the object & hardware attenuation as in the example notebook) as well as scatter correction then I get this (basically the example notebook MLEM without PSF):

[KrisThielemans]

that's interesting. I guess you could try "NAC" with only hardware mu-maps (i.e. without patient) to check if that's it?

Weird thing though is that I don't think we have this with STIR projectors.

[casperdcl]

hmm... I've tried a few permutations as @rijobro describes with no scatter nor PSF (click on images for full size):

Hardware AC	Object AC	Both

Given that the artefacts disappear only if both mu maps are added together before being projected (but are present if only one map is used), it means this is probably not a NiftyPET bug. Really looks like full AC is required...

Are STIR projectors Siddon-based (infinitely narrow LoRs)? NiftyPET uses Siddon so really should be using ~2.5mm FWHM PSF to approximate a more realistic projector - maybe if STIR does this the NAC issue gets hidden?

EDIT: no, looks like NAC+PSF still doesn't work:

[rijobro]

I wonder if there's anything we can test with simulated data. For example, what would we expect if we projected the brainweb data and then NAC reconstructed that?

[pjmark]

Those would be pretty good. I have done lots of them for simulation. But whenever there is discrepancy between attenuation and emission data, you get those artefacts P

…

________________________________ From: Richard Brown <notifications@github.com> Sent: 02 December 2020 19:24 To: NiftyPET/NIPET <NIPET@noreply.github.com> Cc: Markiewicz, Pawel <p.markiewicz@ucl.ac.uk>; Mention <mention@noreply.github.com> Subject: Re: [NiftyPET/NIPET] Vertical streak artifacts in MLEM (#20) I wonder if there's anything we can test with simulated data. For example, what would we expect if we projected the brainweb data and then NAC reconstructed that? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#20 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABPEVHQVILAQ2QSQ7DQWHZTSS2H5HANCNFSM4NYXAXUA>.

[casperdcl]

I wonder if there's anything we can test with simulated data. For example, what would we expect if we projected the brainweb data and then NAC reconstructed that?

could use nipet.simulate_sino(phantom, mu_map, mMRpars, simulate_d=True, mu_input=True) and then poisson() followed by NAC recon...

[KrisThielemans]

Are STIR projectors Siddon-based (infinitely narrow LoRs)? NiftyPET uses Siddon so really should be using ~2.5mm FWHM PSF to approximate a more realistic projector - maybe if STIR does this the NAC issue gets hidden?

STIR's ray-tracing projector does ray-tracing (!). It happens to do it via Siddon algorithm. However, we rarely use it with just a single ray (although that's still the default). Using 5-10 rays improves discretisation artefacts dramatically.

I guess we could do a STIR NAC with only 1 ray if you'd like to compare.

Normally I'd expect that if you simulate and reconstruct with the same projector, everything works nicely...

[casperdcl]

FYI artefacts also visible in NAC OSEM:

I guess we could do a STIR NAC with only 1 ray if you'd like to compare.

yes would be a nice to see done on the same data

[rijobro]

@casperdcl my original reconstruction and the one in your notebook use projectors where the data is copied between the GPU and CPU. But there's also an OSEM implementation that works purely on the GPU. Do you think you'd be able to run an NAC recon using that functionality? I'm sure that the result would be the same if it's down to ray tracing, but I just wanted to confirm.

[casperdcl]

@rijobro my last comment (just above yours) is the NAC OSEM you refer to...

[rijobro]

ah. well good to tick that one off the list at least!

[rijobro]

Just reconstructed the same data using STIR/SIRF. Both NACs with norm and randoms, top row is 1 ray, bottom row is 10 rays.

[KrisThielemans]

interesting... I don't think I've seen this with other scanners, so should be a feature of the gaps that make it worse.

@rijobro what about trying with 100 LORs. shouldn't take too much longer (weirdly enough).

I have no clue though why it would disappear with full AC only. There's obviously going to be similar discretisation artefacts in the AC factors and it's marginally conceivable that they cancel out a little bit, but I'd have expected it then to mostly disappear in the "object AC" recon.

Still, there might be a good thing here: Richard might have wrapped NiftyPET's projectors correctly into STIR!

[pjmark]

This is a tough one. But there are clues here and there. For example, when the bed is removed they never appear. Look at this 24hr scan of Ge68 cylindrical phantom reconstructed without AC:

It looks pure and nice.

[KrisThielemans]

interesting. As this cylinder with head-coil? Likely not.

I can see that the bed/coil introduces high-frequency AC artefacts (and the coil symmetric patterns). I can also see that if you don't do pre-smoothing of your mu-map that you have resolution mismatch with the AC, and therefore have remaining artefacts when you do "hardware AC only", but it's a mystery to me why they then would disappear when adding in the object AC.

@pjmark would you be able to make the cylinder data available? Ideally we'd reconstruct it with STIR (with both its own projector and the NiftyPET-one) as well.

[rijobro]

1, 10, 100 rays:

[KrisThielemans]

ok. nice. so it doesn't go away, but reduces dramatically.

[ashgillman]

One thing I can note in the 1-ray reconstructions is that the artefacts appear to be rotationally symmetrical around the gantry centre, and that the repetition angle seems to be either 90 or 180 degrees.

To me they look a bit like hyperbolas

I don't know the significance of the shape.. But the symmetries about x and y also point towards VOI/LOI errors with the ray tracing. And they're somehow exacerbated by poor AC (or rather somehow compensated with good AC maybe?)

But I think you already worked that out...

[pjmark]

interesting. As this cylinder with head-coil? Likely not.

I can see that the bed/coil introduces high-frequency AC artefacts (and the coil symmetric patterns). I can also see that if you don't do pre-smoothing of your mu-map that you have resolution mismatch with the AC, and therefore have remaining artefacts when you do "hardware AC only", but it's a mystery to me why they then would disappear when adding in the object AC.

@pjmark would you be able to make the cylinder data available? Ideally we'd reconstruct it with STIR (with both its own projector and the NiftyPET-one) as well.

Yes, but the data size is >100GB so difficult to share. Any ideas?

[KrisThielemans]

the data size is >100GB so difficult to share. Any ideas?

It looks like OneDrive can do 100GB per file these days. Nevertheless probably best to split. see e.g. https://www.baeldung.com/linux/split-files