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 crystal for ml_norm 3d #713

Closed
wants to merge 2 commits into from
Closed

Conversation

gefeichen
Copy link
Contributor

new_max_delta in ml_norm maybe not correct.
when reviewing, ignore whitespace and empty lines should be chosen.

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.

thanks. some detailed comments. 2 general ones.

  • I'd rather not introduce a lot of white-space differences in this PR. I'd rather do that by fixing fix and enforce convention on white-spaces, code-indentation and other coding style #98 (see configuration for white-space #724). Would it be possible to revert your white-space differences, or we wait till configuration for white-space #724 is merged (could take a bit)
  • I'd rather replace make_fan_data code with the one from make_fan_data_remove_gapsetc and remove the explicit*remove_gapsfunctions. That'd mean we can only do this if the*remove_gaps` functions give exactly the same result if there are no virtual crystals. We have no tests for anything like this unfortunately so I guess you have to do it by hand or even create a small test first! (create text files with norm factors, create some projdata (e.g. just 1), find norm factors, check)

const int new_fan_size = fan_size - num_transaxial_blocks_in_fansize*virtual_transaxial_crystals;
const int new_half_fan_size = new_fan_size/2;
const int num_axial_blocks_in_max_delta = max_delta/(num_axial_crystals_per_block);
const int new_max_delta = max_delta - (num_axial_blocks_in_max_delta)*virtual_axial_crystals-1;
Copy link
Collaborator

Choose a reason for hiding this comment

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

this formula fits the previous one, but it seems wrong. If virtual_axial_crystals==0, it results in new_max_delta = max_delta - 1, which would of course be really problematic if max_delta==0.

I believe we should not have the -1.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Possibly this came about as @tniknejad added a virtual ring at the end (while the mCT doesn't). Might need a bit more thought on how to handle that.

Comment on lines +1072 to +1078
// std::cerr<<new_num_rings<<std::endl;
// std::cerr<<new_num_detectors_per_ring<<std::endl;
// std::cerr<<new_max_delta<<std::endl;
// std::cerr<<new_fan_size<<std::endl;


// **** End **** //
Copy link
Collaborator

Choose a reason for hiding this comment

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

remove these


shared_ptr<SegmentBySinogram<float> > segment_ptr;
Bin bin;
for (bin.segment_num() = proj_data.get_min_segment_num(); bin.segment_num() <= proj_data.get_max_segment_num(); ++ bin.segment_num()) {
Copy link
Collaborator

Choose a reason for hiding this comment

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

don't change indentation/bracing scheme

for (bin.view_num() = 0; bin.view_num() < num_detectors_per_ring / 2; bin.view_num()++)
for (bin.tangential_pos_num() = -half_fan_size;
bin.tangential_pos_num() <= half_fan_size;
++bin.tangential_pos_num()) {
Copy link
Collaborator

Choose a reason for hiding this comment

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

don't change indentation/bracing scheme

continue;
}

int offset = (a%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(a%num_transaxial_crystals_per_block%new_num_transaxial_crystals_per_block):0;
Copy link
Collaborator

@KrisThielemans KrisThielemans Oct 22, 2020

Choose a reason for hiding this comment

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

can you clarify the rationale for the offset? I'd rather have it zero, I believe this is how it's done in #617. (There is of course a question if the virtual crystal is put "before" or "after" the block in transaxial direction, but I'd rather resolve that with #181. Otherwise we need yet another variable. In any case, that's now what you seem to have in mind here).

Note that in your offset formula you have a%n%n which I guess is a%n.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The position of virtual crystal should be considered and I have another question: if the number of virtual crystal is 0 or 1, it's not necessary to set this offset.

12(transaxial crystal)+3(virtual crystal) = 15

if index is 10
10%15=10 < 12 offset = 0
index is 14
14%15 = 14> 12 offset = 14%15%12 = 14%12 = 2
This is set for multiple virtual crystals.

Copy link
Collaborator

Choose a reason for hiding this comment

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

from my comment above, index=14 is a virtual crystal in your example, so should not be stored. (the bin-value will be zero, and should not be entered in the data without gaps).

Continuing your example, the first non-virtual crystal in the next block would have index=15 (as zero-based). it should be assigned to new_index=12. Using the formula below without offset

int new_a = a - a/num_transaxial_crystals_per_block*virtual_transaxial_crystals;

we'd get, new_a = 15 - (15/15)*3 = 12, which is correct.

Similarly, if a=16, new_a=16-3=13. Going to the next block a=30, new_a=30-(30/15)*3=3--6=24.

So I believe no offset needed in all cases.

if ((rb+1)% num_axial_crystals_per_block == 0) continue;
if ((a+1)% num_transaxial_crystals_per_block == 0) continue;
if ((b+1)% num_transaxial_crystals_per_block == 0) continue;
if ((ra == num_rings - 1) || (rb == num_rings - 1) ) continue;
Copy link
Collaborator

Choose a reason for hiding this comment

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

same comment as above. I don't think this should be here


// const int num_axial_crystals_per_block = num_axial_detectors/num_axial_blocks; // number of axial detector per block in fan_data w/o gaps
// const int num_transaxial_crystals_per_block = num_transaxial_detectors/num_transaxial_blocks; // number of transaxial detector per block in fan_data w/o gaps
int offset = (a%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(a%num_transaxial_crystals_per_block%new_num_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.

same comment as above

Comment on lines +90 to +94
const int num_rings =
measured_data->get_proj_data_info_sptr()->get_scanner_sptr()->
get_num_rings() -(num_axial_blocks-1)*virtual_axial_crystals;
const int num_detectors_per_ring =
measured_data->get_proj_data_info_sptr()->get_scanner_sptr()->
Copy link
Collaborator

Choose a reason for hiding this comment

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

can we call these variables new_* or (better) non_virtual_*? even better would be to add member functions that do this calculation in Scanner

Copy link
Collaborator

@KrisThielemans KrisThielemans Feb 10, 2021

Choose a reason for hiding this comment

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

Let's call these variables num_physical_num_detectors_per_ring etc but leave corresponding adding member functions to #776

const int num_axial_blocks =
measured_data->get_proj_data_info_sptr()->get_scanner_sptr()->
get_num_axial_blocks();
const int virtual_axial_crystals =
Copy link
Collaborator

Choose a reason for hiding this comment

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

as above

Comment on lines +120 to +123
const int num_rings =
measured_data->get_proj_data_info_sptr()->get_scanner_sptr()->
get_num_rings() -(num_axial_blocks-1)*virtual_axial_crystals;
const int 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.

as in apply*

@KrisThielemans
Copy link
Collaborator

@gefeichen any chance you can have another look at this?

@gefeichen
Copy link
Contributor Author

@gefeichen any chance you can have another look at this?

@gefeichen any chance you can have another look at this?

sorry for the later reply. I will finish it by December

Comment on lines +1099 to +1102
if ((ra == num_rings - 1) || (rb == num_rings - 1) || (a == num_detectors_per_ring - 1) ||
(b == num_detectors_per_ring - 1)) {
continue;
}
Copy link
Collaborator

Choose a reason for hiding this comment

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

ok. I agree that for the mCT this condition is dangerous.

I believe that what this condition attempts to do is to detect if it's a virtual crystal. If so, it shouldn't set the fan-data at all. However, the condition is flawed. What about the following.

Add something like this (but adding 2 more args)

bool is_virtual_crystal(int a, int ra)
{
   const int a_in_block = a % num_transaxial_crystals_in_block;
   if (a_in_block >= num_non_virtual_transaxial_crystals_in_block)
     return true;
   const int ra_in_block = ra % num_axial_crystals_in_block;
   return (ra_in_block >= num_non_virtual_axial_crystals_in_block)
}

Then the condition should become

  if (is_virtual_crystal(a,ra) || is_virtual_crystal_b,rb))
    continue;

const int new_fan_size = fan_size - num_transaxial_blocks_in_fansize*virtual_transaxial_crystals;
const int new_half_fan_size = new_fan_size/2;
const int num_axial_blocks_in_max_delta = max_delta/(num_axial_crystals_per_block);
const int new_max_delta = max_delta - (num_axial_blocks_in_max_delta)*virtual_axial_crystals-1;
Copy link
Collaborator

Choose a reason for hiding this comment

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

Possibly this came about as @tniknejad added a virtual ring at the end (while the mCT doesn't). Might need a bit more thought on how to handle that.

continue;
}

int offset = (a%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(a%num_transaxial_crystals_per_block%new_num_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.

from my comment above, index=14 is a virtual crystal in your example, so should not be stored. (the bin-value will be zero, and should not be entered in the data without gaps).

Continuing your example, the first non-virtual crystal in the next block would have index=15 (as zero-based). it should be assigned to new_index=12. Using the formula below without offset

int new_a = a - a/num_transaxial_crystals_per_block*virtual_transaxial_crystals;

we'd get, new_a = 15 - (15/15)*3 = 12, which is correct.

Similarly, if a=16, new_a=16-3=13. Going to the next block a=30, new_a=30-(30/15)*3=3--6=24.

So I believe no offset needed in all cases.

}

int offset = (a%num_transaxial_crystals_per_block)>new_num_transaxial_crystals_per_block?(a%num_transaxial_crystals_per_block%new_num_transaxial_crystals_per_block):0;
int new_a = a - a/num_transaxial_crystals_per_block*virtual_transaxial_crystals - offset;
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
int new_a = a - a/num_transaxial_crystals_per_block*virtual_transaxial_crystals - offset;
int new_a = a - (a/num_transaxial_crystals_per_block)*virtual_transaxial_crystals;

I'd feel safer with brackets.

Copy link
Contributor Author

@gefeichen gefeichen Dec 21, 2020

Choose a reason for hiding this comment

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

bool is_virtual_crystal(int a, int ra)
{
   const int a_in_block = a % num_transaxial_crystals_in_block;
   if (virtual_crystal_transaxial_after_block){
   if (a_in_block >= num_non_virtual_transaxial_crystals_in_block)
     return true;
   }
   else{
   if (a_in_block <= num_of_virtual_transaxial_crystals_in_block)
     return true;
   }
   const int ra_in_block = ra % num_axial_crystals_in_block;
   if (virtual_crystal_axial_after_block){
      return (ra_in_block >= num_non_virtual_axial_crystals_in_block)
   else
       return (ra_in_block <= num_of_virtual_axial_crystals_in_block)
   }
}

it is easy to add virtual_crystal_transaxial_before_or_after_block. And it can deal with outermost virtual rings in the axial.
Maybe it's not necessary to do such things now.
I think I can do some tests for mCT with and without virtual crystals.

@KrisThielemans
Copy link
Collaborator

@gefeichen @tniknejad shall we discuss this in a tcon? Please send me and @NikEfth an email with your availability.

@gefeichen gefeichen closed this Feb 24, 2021
@gefeichen
Copy link
Contributor Author

refer to #833

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

Successfully merging this pull request may close these issues.

2 participants