virtual crystal for ml_norm 3d by gefeichen · Pull Request #713 · UCL/STIR

gefeichen · 2020-10-14T08:26:06Z

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

…ormfactors3D.cxx

1

KrisThielemans

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 #98 (see #724). Would it be possible to revert your white-space differences, or we wait till #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)

KrisThielemans · 2020-10-22T06:10:46Z

src/buildblock/ML_norm.cxx

+        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;


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.

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.

KrisThielemans · 2020-10-22T06:13:38Z

src/buildblock/ML_norm.cxx

+//    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     ****  //


remove these

KrisThielemans · 2020-10-22T06:14:11Z

src/buildblock/ML_norm.cxx

+
+        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()) {


don't change indentation/bracing scheme

KrisThielemans · 2020-10-22T06:14:36Z

src/buildblock/ML_norm.cxx

+                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()) {


don't change indentation/bracing scheme

KrisThielemans · 2020-10-22T06:35:31Z

src/buildblock/ML_norm.cxx

+                            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;


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.~~

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.

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.

KrisThielemans · 2020-10-22T06:42:02Z

src/buildblock/ML_norm.cxx

+                        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;


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

KrisThielemans · 2020-10-22T06:42:20Z

src/buildblock/ML_norm.cxx

-
-   // 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;


same comment as above

KrisThielemans · 2020-10-22T06:45:17Z

src/utilities/apply_normfactors3D.cxx

+    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()->


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

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

KrisThielemans · 2020-10-22T06:46:29Z

src/utilities/find_ML_normfactors3D.cxx

+    const int num_axial_blocks =
+            measured_data->get_proj_data_info_sptr()->get_scanner_sptr()->
+                    get_num_axial_blocks();
+    const int virtual_axial_crystals =


KrisThielemans · 2020-10-22T06:47:02Z

src/utilities/find_ML_normfactors3D.cxx

+    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 =


as in apply*

KrisThielemans · 2020-11-15T17:17:05Z

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

gefeichen · 2020-11-16T12:21:04Z

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

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

KrisThielemans · 2020-12-19T08:28:29Z

src/buildblock/ML_norm.cxx

+                        if ((ra == num_rings - 1) || (rb == num_rings - 1) || (a == num_detectors_per_ring - 1) ||
+                            (b == num_detectors_per_ring - 1)) {
+                            continue;
+                        }


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;

KrisThielemans · 2020-12-19T08:30:12Z

src/buildblock/ML_norm.cxx

+        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;


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.

KrisThielemans · 2020-12-19T08:52:21Z

src/buildblock/ML_norm.cxx

+                            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;


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.

KrisThielemans · 2020-12-19T08:56:39Z

src/buildblock/ML_norm.cxx

+                        }
+
+                        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;


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.

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 · 2021-02-04T18:47:23Z

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

gefeichen · 2021-02-24T08:04:26Z

refer to #833

gefeichen and others added 2 commits October 14, 2020 16:04

add virtual crystal support for find_ML_normfactors3D.cxx and apply_n…

76f1ba1

…ormfactors3D.cxx

Merge pull request #1 from UCL/master

01857ba

1

KrisThielemans requested changes Oct 22, 2020

View reviewed changes

KrisThielemans reviewed Dec 19, 2020

View reviewed changes

gefeichen closed this Feb 24, 2021

	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;

Conversation

gefeichen commented Oct 14, 2020

Uh oh!

KrisThielemans left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans Oct 22, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans Feb 10, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans commented Nov 15, 2020

Uh oh!

gefeichen commented Nov 16, 2020

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gefeichen Dec 21, 2020 • edited by KrisThielemans Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans commented Feb 4, 2021

Uh oh!

gefeichen commented Feb 24, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

KrisThielemans Oct 22, 2020 •

edited

Loading

KrisThielemans Feb 10, 2021 •

edited

Loading

gefeichen Dec 21, 2020 •

edited by KrisThielemans

Loading