Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Virtual xtal interface. #776

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
Open

Conversation

NikEfth
Copy link
Collaborator

@NikEfth NikEfth commented Dec 22, 2020

Interface for the virtual xtals.

Update CListModeDataROOT.cxx, Scanner.h, and 4 more files...

Update CListModeDataROOT.cxx, Scanner.h, and 4 more files...
@NikEfth
Copy link
Collaborator Author

NikEfth commented Dec 22, 2020

This PR addresses the issues raised in #750

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most of this is fine. However, you're changing the interpretation of the "standard" function like get_num_axials_crystal_per_block to mean the "non-virtual" ones, while it was "all of them". Although I agree that the naming is confusing, I'd rather avoid changing this as it probably has much wider consequences. For instance, what does get_num_rings mean? It was "number of all rings, non-virtual+virtual". If we'd change that, lots of things will have to be adapted.

src/IO/InterfileHeader.cxx Outdated Show resolved Hide resolved
@@ -160,7 +160,7 @@ Scanner::get_num_axial_blocks() const
int
Scanner::get_num_transaxial_blocks() const
{
return num_detectors_per_ring/num_transaxial_crystals_per_block;
return num_detectors_per_ring/(num_transaxial_crystals_per_block+num_virtual_transaxial_crystals_per_block);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cannot change this.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I still disagree on this one. Currently num_transaxial_crystals_per_block includes both. We shouldn't be adding.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then in come down to what we include in the total number of detectors in the template.
If the total number of detectors include the virtual or not.
My initial notion is that they should not. BUT if you don't you have problems with the get_bin() function (wrong inversions).

So, until we restructure the class, the total number of rings and detectors per ring should include the virtual.

ok = false;
}

if (scanner_sptr->get_num_detectors_per_ring() != root_file_sptr->get_num_dets_per_ring())
if ((scanner_sptr->get_num_detectors_per_ring() - scanner_sptr->get_num_virtual_transaxial_crystals()) != root_file_sptr->get_num_dets_per_ring())
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cannot change this

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess we could change it after all, if we say that the root_file_sptr->get_num_dets_per_ring() returns the number of physical detectors per ring, which is probably logical.

Needs serious documentation though

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But I would still find it very confusing if get_num_dets_per_ring() means different things for Scanner and the ROOT stuff

// default:
// return 0;
// }
return num_virtual_axial_crystals_per_block * get_num_axial_blocks();
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is probably going to be incorrect. The mCT doesn't use a virtual crystal for the last block of course (apparently Tahereh did for her scanner).

Would have to be derived from num_rings somehow. Do we need this function?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I know.
I don't know what to propose here. Say that you add one that the end. I don't think that it actually changes much. Where could this create a problem? Normalisation? Multiple bed positions? I guess the axial middle of the scanner will include that ring and that might be a problem.

What we could do is check with the num_rings + gaps - 1,

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@tniknejad Is this something you could change?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm very confused. you introduced this function. We have get_num_virtual_axial_crystals_per_block which is well-defined. why do we need this?

src/buildblock/Scanner.cxx Outdated Show resolved Hide resolved
@@ -746,12 +790,12 @@ check_consistency() const
{
const int dets_per_ring =
get_num_transaxial_blocks() *
get_num_transaxial_crystals_per_block();
(get_num_transaxial_crystals_per_block() + get_num_virtual_transaxial_crystals_per_block());
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't be changed

@@ -782,7 +826,7 @@ check_consistency() const
{
const int dets_axial =
get_num_axial_blocks() *
get_num_axial_crystals_per_block();
(get_num_axial_crystals_per_block() + get_num_virtual_axial_crystals_per_block());
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't be changed I believe

@@ -817,12 +861,12 @@ check_consistency() const
this->get_name().c_str());
else
{
if ( get_num_detectors_per_ring() % get_num_transaxial_crystals_per_singles_unit() != 0)
if ( get_num_detectors_per_ring() % (get_num_transaxial_crystals_per_singles_unit() + get_num_virtual_transaxial_crystals_per_block()) != 0)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't be changed

{
warning("Scanner %s: inconsistent transaxial singles unit info:\n"
"\tnum_detectors_per_ring %d should be a multiple of num_transaxial_crystals_per_singles_unit %d",
"\tnum_detectors_per_ring %d should be a multiple of num_transaxial_crystals_per_singles_unit + virtual %d",
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't be changed

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But this is what you have in mCT.
This is an important comment: If you don't include the virtual crystals in the total number of detectors per ring the unlisting fails. The sinogram is broken. I could send an image to show you.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It has something to do with the inversion in the tangestial positions when the tang_pos is larger than the number of views. Which currently is calculated from the total number of detectors per ring.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you don't include the virtual crystals in the total number of detectors per ring the unlisting fails. The sinogram is broken.

of course. Siemens does all calculations/sinograms in terms of "total number". Nowhere in their files there are any "virtual" ones (even not in the norm, presumably because they do gap-filling)

this->get_name().c_str(),
get_num_detectors_per_ring(), get_num_transaxial_crystals_per_singles_unit());
get_num_detectors_per_ring(), (get_num_transaxial_crystals_per_singles_unit() + get_num_virtual_transaxial_crystals_per_block()));
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't be changed

@NikEfth
Copy link
Collaborator Author

NikEfth commented Dec 22, 2020

The definition of the scanner overall is very circular and confusing and probably we should re-write from scratch.

However, I took the lead from the mCT example were the number of the detectors per ring includes the virtual crystals (13+1)*48. Respectively the num_rings should include the virtual rings.

The changes, I made allow us to design a scanner if a 2-fold way, 1 gap per block (in mCT every 13 crystals) or 1 gap every several blocks. The later is in better agreement with most scanners I would know, where the block separation between blocks is small of e.g. half crystal, but there is 1 or 2 crystal between the rotated (bigger) blocks (r-sectors in GATE).

For example, my scanner has 18 sectors with 4 transaxial blocks with 8 transaxial crystals.
Now both number of detectors per ring 18*(4*8+1) and 18*4*(8+1) can pass the consistency test.
And soon the CListModeDataROOT will be able to know were to push the gap.

I think that bein able to reposition the gaps is quite important.

@NikEfth
Copy link
Collaborator Author

NikEfth commented Dec 22, 2020

I would like to add that the modification does not interfere with the singles units, which can still be independent.

@KrisThielemans
Copy link
Collaborator

The definition of the scanner overall is very circular and confusing and probably we should re-write from scratch

yeah. at some point... You cleaned at least the set_params stuff up in #304. Let's not attempt a major rewrite until #304 and #577 are merged.

The changes, I made allow us to design a scanner if a 2-fold way, 1 gap per block (in mCT every 13 crystals) or 1 gap every several blocks. The later is in better agreement with most scanners I would know, where the block separation between blocks is small of e.g. half crystal, but there is 1 or 2 crystal between the rotated (bigger) blocks (r-sectors in GATE).

I certainly agree that many scanners would have bigger gaps between "buckets" than "blocks". However, I hadn't realised that you wanted to fix that with this PR. The name of the PR and previous conversation led me to think you only wanted to be able to set the number of virtual crystals (which is what most of this PR is about).

In fact, I cannot see how you could fix this without having additional members num_axial_virtual_crystals_per_bucket etc (or proably better num_additional_axial_virtual_crystals_per_bucket as it'd add up to the one between blocks).

To me adding that facility seems to need a fair amount of work, but maybe I'm wrong.

@NikEfth
Copy link
Collaborator Author

NikEfth commented Dec 22, 2020

Well, I am not really fixing it, for that we would need a num_virtual_crystal_per_bucket() function, but just letting the consistency test pass and the make the CListModeROOT do the right thing.
The key logic is that based on the number of detector per ring you set, we can figure out how many virtual crystals there are. Is the divisor of block number or block*crystal.

I understand that some tests succeed but there might be further implications down the line.

NikEfth and others added 3 commits December 22, 2020 13:58
Co-authored-by: Kris Thielemans <KrisThielemans@users.noreply.github.com>
Co-authored-by: Kris Thielemans <KrisThielemans@users.noreply.github.com>
@NikEfth
Copy link
Collaborator Author

NikEfth commented Dec 23, 2020

hmm, after all, I introduced the gaps in buckets.
My main question is whether you prefer the num_of_detectors_per_ring in the template to be the physical detectors only, or both gaps and detectors
I think that the second is easier. Otherwise, we have to go everywhere and replace the get_numXX() functions with get_num_logicalXX()
While if the total number in the template includes the virtual crystals we have to introduce functions get_num_phisicalXX() in fewer cases

Actually, I think that in the template we should set the number of physical detectors and the then get functions to do the calculations. However, this means that we cannot have set functions. There have to become explicitly set_physical, or very very very well documented. What do you think?

int
Scanner::get_num_rings() const
{  return get_num_physical_rings() + (get_num_virtual_axial_crystals_per_bucket() - 1) *  get_num_axial_buckets() +
            (get_num_virtual_axial_crystals_per_block()-1) * get_num_axial_blocks();}
int
Scanner::get_num_detectors_per_ring() const
{
  return num_detectors_per_ring + (get_num_virtual_transaxial_crystals_per_bucket() - 1) * get_num_transaxial_buckets() +
          (get_num_virtual_transaxial_crystals_per_block()-1) * get_num_transaxial_blocks();}
int
Scanner::get_num_physical_rings() const
{  return num_rings; 
}
int
Scanner::get_num_physical_detectors_per_ring() const
{
  return num_detectors_per_ring;
}

@KrisThielemans
Copy link
Collaborator

I think that in the template we should set the number of physical detectors and the then get functions to do the calculations. However, this means that we cannot have set functions. There have to become explicitly set_physical, or very very very well documented. What do you think?

While reading/setting things in terms of physical detectors would make a lot of sense, changing the meaning of get_num_rings() to suddenly be get_num_physical_rings() would break quite a lot of stuff. Given my pace of accepting PRs, I'd strongly recommend not attempting that... Of course, making it clearer in the documentation what is what would help.

Of course, having members get_num_physical_rings(), get_num_non_physical_rings() (with the sum get_num_rings() ) is an excellent idea. Providing an alias for get_num_rings() that makes it clear what it means is of course fine as well. We could then consider deprecating the old names.

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've given this another look and added some comments. I didn't check consistency of everything. I also didn't think about the "singles units" stuff.

I guess I'm fine with changing the internal meaning of some of the private variables like Scanner::num_rings if you prefer, but I don't think we should change the meaning of get_num_rings() etc (which you don't of course, so that's good). It could be confusing for a developer of Scanner, but it can be addressed by documentation.

I really think though that we cannot change the meaning of num_rings etc in the Interfile header. It's too late for that.

inline int get_num_rings() const;
//! get the number of detectors per ring
//! get the number of detectors per ring. This function returns the total
//! number of detectors + gasp in the ring.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
//! number of detectors + gasp in the ring.
//! number of (physical and virtual) detectors in the ring.

@@ -188,9 +192,10 @@ class Scanner
//! \name Functions returning geometrical info
//@{

//! get number of rings
//! get number of rings. This function returns the total number of rings + gaps
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
//! get number of rings. This function returns the total number of rings + gaps
//! get number of rings. This function returns the total number of (physical and virtual) rings

@@ -302,9 +311,11 @@ class Scanner
// zlong, 08-04-2004, add set_methods
//! set scanner type
inline void set_type(const Type & new_type);
//! set number of rings
//! set number of rings. This function will set the total number of rings
//! reseting the virtual crystals to the axial direction to 0.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
//! reseting the virtual crystals to the axial direction to 0.
//! resetting the virtual crystals in the axial direction to 0.

Comment on lines +317 to +318
//! set the namber of detectors per ring. This function will set the total number
//! of detectors per ring, reseting the number of virtual crystals to 0.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
//! set the namber of detectors per ring. This function will set the total number
//! of detectors per ring, reseting the number of virtual crystals to 0.
//! set the number of detectors per ring. This function will set the total number
//! of detectors per ring, resetting the number of virtual crystals (in transaxial direction) to 0.

}

int
Scanner::get_num_transaxial_blocks() const
{
return num_detectors_per_ring/num_transaxial_crystals_per_block;
return num_transaxial_crystals_per_block == 0 ? num_detectors_per_ring :
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

num_transaxial_crystals_per_block should never be zero. It doesn't make sense in STIR. You might have only 1 block, but no crystals in a block means no crystals

}

int
Scanner::get_num_axial_buckets() const
{
return get_num_axial_blocks()/num_axial_blocks_per_bucket;
return num_axial_blocks_per_bucket == 0 ? get_num_axial_blocks() :
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as above

ok = false;
}

if (scanner_sptr->get_num_detectors_per_ring() != root_file_sptr->get_num_dets_per_ring())
if ((scanner_sptr->get_num_detectors_per_ring() - scanner_sptr->get_num_virtual_transaxial_crystals()) != root_file_sptr->get_num_dets_per_ring())
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But I would still find it very confusing if get_num_dets_per_ring() means different things for Scanner and the ROOT stuff

@@ -600,6 +600,12 @@ InterfilePDFSHeader::InterfilePDFSHeader()
num_transaxial_crystals_per_block = 0;
add_key("number of crystals_per_block in transaxial direction",
&num_transaxial_crystals_per_block);
num_virtual_axial_crystals_per_block = 0;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

seems we'd need some bucket stuff here as well

@KrisThielemans
Copy link
Collaborator

@NikEfth I cannot recall where we were with this, but aside from the comments above, it probably needs a revisit now that #833 has been merged.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants