jacg / petalo

PET scanner image reconstruction
1 stars 4 forks source link

Sensor simulation #22

Closed andLaing closed 1 year ago

andLaing commented 2 years ago

Implementation of reconstruction of LORs from MC reconstruction

andLaing commented 2 years ago

More or less all the names will have to be improved. So far a very basic reconstruction is implemented. Still need to implement the DOI correction and KMeans clustering among other things.

jacg commented 2 years ago

Do you want to talk about this urgently, or can it wait a few days? (I'm in the middle of an optimization which is ploughing through much more code than I initially hoped, and I'd rather get that finished before I get too lost.)

andLaing commented 2 years ago

Nothing urgent. Just figured there was enough that it could go up in the draft PR

andLaing commented 2 years ago

Seem to have developed a magical mystery test fail: sum_of_weights_equals_length_through_box

On master too in my local.

jacg commented 2 years ago

Yes, this is irksome. It's (almost certainly) a spurious floating point inaccuracy failure found by proptest. It remembers (the hashes needed to generate) the counterexamples it finds in proptest-regressions/*.txt. In your case it's probably weights.txt. Then it uses those to check the couterexamples before generating any other cases. Similar to what hypothesis does. So just remove the last line from proptest-regressions/weights.txt and it should work again ... until it finds another spurious counterexample. The joy of floating point arithmetic. Dealing with this is on my TODO list, but it's not in any danger of reaching the top any time soon.

andLaing commented 2 years ago

The commit history of this branch is a complete shambles at the moment. If there's no complaints I'm going to rebase and do some tidying.

jacg commented 2 years ago

No problem, go ahead.

jacg commented 2 years ago

Is there anything I can do to help this PR along?

I'm a bit worried that, as I clean up, refactor and reorganize, master will diverge more and more from this branch, making merging more difficult.

andLaing commented 2 years ago

I've not been working on it so much the last week but I think it's in a state that I can do some comparisons with the other code. Maybe it's worthwhile to try and get the work so far integrated and then add things like the KMEANS separator option in a later PR.

jacg commented 2 years ago

It would be good to get it merged. Let me know if I can help in any way.

andLaing commented 1 year ago

In it's current state (only the dot product separation for clustering, only zrms correction for DOI, ...) I've compared with the existing Julia code and I get the same results within error for y, z, and dt. For x there's maybe a slight difference at the edges which I'm looking into.

I think the code can be reviewed now.

jacg commented 1 year ago

Have you any idea of how the run-time and memory usage compare to the Julia version?

jacg commented 1 year ago

What datasets would you suggest I use to take this for a spin?

jacg commented 1 year ago

I would like to merge as soon as

I think it's important to have the whole reconstruction chain available in this repository, and this is the only missing piece. If it does the job, let's merge. Further tweaks and improvements can be done later.

andLaing commented 1 year ago

Have you any idea of how the run-time and memory usage compare to the Julia version?

There isn't an identical Julia executable. Julia does the intermediate analysis to be able to compare with MC truth and get calibrations whereas here we just use the calibration to go straight to LORs.

I'd have to write another script (it's been on the todo for a while) in Julia to do the same for performance comparisons.

jacg commented 1 year ago

Hmm, OK, I'd rather not merge something that has horrible performance characteristics, so we'll need to do at least a little bit of measurement to make sure that it's acceptable. Further optimization can be done later.

andLaing commented 1 year ago

What datasets would you suggest I use to take this for a spin?

I was using the "sanity" run that you gave me at the beginning of the year (think) when I was looking into the biases in the reconstruction. I think it's still on one of the servers, I'll have a look and if not I can copy my local files.

I would like to merge as soon as

The HDF5 module names have been fixed
You are happy that the results aren't wrong.
We are convinced that it's not significantly slower than the Julia version.

I think it's important to have the whole reconstruction chain available in this repository, and this is the only missing piece. If it does the job, let's merge. Further tweaks and improvements can be done later.

So I'll try and get to the name changes today.

The results look good, modulo the limitations on the options for DOI reconstruction and time parameter (here only minimum, in julia you can use an average). I'll pass you the plots later. There was one difference I wanted to look at more closely but I suspect it's not a show stopper.

As I mentioned above, there isn't exact equivalence because the Julia looks through the vertex and primary tables, saves reco and truth parameters and allows for a calibration etc then there's a separate script that makes the LORs from that 'reduced' output. At the moment the RUST version is much faster but as it does less it's not the fairest of comparisons.

jacg commented 1 year ago

How much work would it take to port the missing part to Rust? Or is it something that requires an interactive environment for human interaction?

jacg commented 1 year ago

the limitations on the options for DOI reconstruction and time parameter (here only minimum, in julia you can use an average).

Any difficulties here, or just a question of getting around to it?

andLaing commented 1 year ago

It depends which "missing part" you mean. There's the other cluster algorithm (KMeans, I have another branch where I've been doing some tests with the linfa_clustering::KMeans implementation), the average time and the DOI options. Most is just getting round to doing it with the possible exception of allowing the use of average time which might need a change in approach.

I'd prefer to do these additions in separate PRs since this one is already a bit of a monster and contains the things most likely to cause rebase conflicts (moving existing code to different files etc).

jacg commented 1 year ago

Anything else would definitely go in a separate PR. I just wanted to have some idea of how much is still missing.

It depends which "missing part" you mean.

In the first instance, the part which calculates the calibrations. I want to be able to do a complete reconstruction entirely using code in this repo. From this perspective, alternative algorithms are less urgent than missing links in the chain.

andLaing commented 1 year ago

In the first instance, the part which calculates the calibrations. I want to be able to do a complete reconstruction entirely using code in this repo. From this perspective, alternative algorithms are less urgent than missing links in the chain.

Calculating the calibrations requires not only the reduced data with mc truth saved but energy filtering, fitting distributions etc. I hadn't even put this on my TODO.

jacg commented 1 year ago

Does this part require human interaction, or can it be automated completely?

andLaing commented 1 year ago

I've tended to supervise but it's more or less automated.

jacg commented 1 year ago

I've added a few stylistic changes. Further such fiddling can wait until I update the toolchain to the recently-released Rust 1.64.

Is there anything else you want to do here before we merge?

andLaing commented 1 year ago

I think that's all that I thought to do for this first PR. The outputs are statistically identical to what I get from Julia with the exception of x where there seems to be a bit of a difference at the edges. I suspect that it's a slight difference in the application of one of the selection criteria but it doesn't seem like a big enough difference to worry about at this level (see plots attached). It can be checked again in the subsequent PRs. rust-julia_dt rust-julia_q1 rust-julia_x1 rust-julia_y1 rust-julia_z1