Help me: Offsetscan

[ProjX0]

Thank you for sharing the excellent toolkit.

I am currently interested in offset scan and am trying it with the data I have. However, I seem to have made a mistake, and the reconstruction isn’t working well.

Could you possibly help me if I share the geometry information and data with you?

Thank you very much.

[kylechampley]

Sure. Please also send a picture of one of the projections.

[hws203]

I guess that your right_side offscan mechanism has some issue. I tested your "d23_offscan.py" which shows some strange reconed slice and like attached one. I just change the +100 offset by -100 as like centerCols = 0.5*(numCols-1)-100 and also change numCols by 512 too. Please, check it by changing the offset value by -100 and the numCols by 512.

[ProjX0]

For reference, the images should look like the one above after reconstruction. I am currently trying to achieve this using LEAP, but it appears that I am making some mistakes, as the images are not coming out as intended.

[kylechampley]

Thank you @hws203 for pointing out that bug. The cause for this was quite minor and was easy to fix. Currently the bug fix is in my development branch champley_dev. This branch will be merged into the main branch and a new LEAP release will come out later this week. You can either use this new code now or there is a workaround if you want to just work with the current version. The work-around is just add the following command after you set the reconstruction volume parameters: leapct.set_diameterFOV(leapct.get_voxelWidth()*leapct.get_numX()) The issue you noticed was caused by LEAP calculating the field of view incorrectly. Thus the above line forces LEAP to not restrict the field of view, so you can see the whole thing.

OK, now for @ProjX0 specific issue. Could you post your whole LEAP script? Above you have only shown a portion of it. One question I have is where the "shift_distance" value came from. I did not see this in the geometry file you posted. Also, why do you shift centerCol and centerRow by this same amount? I also don't understand the parameter names in that file like p0.x, p0.y, phi, etc. Can you explain what these mean?

[ProjX0]

In our half-beam geometry, the centerline of the source is positioned 38.139 mm from the center of the detector.

p0.x, p0.y: coordinates of the center of the image on the detector. Phi (φ): rotation angle of the detector relative to the x-axis. Eta (η): rotation angle of the detector relative to the y-axis. Theta (θ): rotation angle of the gantry.

The parameters in the geometry.txt file are currently fixed for testing with SOD = 410 mm and SDD = 660 mm, and the other variables have not been used.

Here is the whole LEAP script:

[kylechampley]

Thanks for the code. I was able to make some improvements. I fixed some more small bugs in the code, so please use the code in the champley_dev branch until I merge in into the main branch this weekend. The main two things I changes in the processing were:

1) Change the direction of the rotation. Instead of an angular range of 360.0, just put -360.0.

2) Crop off the border of the projection images which contain bad pixels. With the binned data, I did it like this: g = leapct.crop_projections([2, leapct.get_numRows()-3], [2, leapct.get_numCols()-3], g)

With those changes, here is what I got.

With those changes, it is still not a great result, but I noticed in the geometry file you posted that from projection-to-projection there seemed to be a lot of variation in the geometry. You'll need to use LEAP's "modular-beam" geometry and specify all those parameters. I did not know how to interpret all those parameters and information from some of the measured projections seemed to be missing from this file, so I couldn't try this myself.

[kylechampley]

Oh, one more thing. Strictly speaking, this data cannot be properly reconstructed with the "offset scan" strategy because many of the projections are truncated on both sides. Offset scan implies that the projections are only truncated on one side. Yes, this is true for most of your projections, but there are enough projections that violate this, so you'll still get some artifacts. Regardless, I think the next thing you need to do is use the modular-beam geometry as I discussed in my last post.

[ProjX0]

I first downloaded the code from the champley_dev branch, deleted the existing version, and reinstalled it.

I then made the changes as you suggested:

1. Change the direction of the rotation. Instead of an angular range of 360.0, just put -360.0. leapct.set_conebeam(numAngles-start_index, numRows, numCols, pixelSize, pixelSize, centerRow, centerCol, leapct.setAngleArray(numAngles-start_index, -360.0), 410, 660)

However, I am still encountering problems. Is there something I missed?

For reference, the whole LEAP script is as follows:

If this works, I will thoroughly examine the projection data and geometry information, and then attempt it using the modular-beam geometry. Thank you.

[ProjX0]

I apologize, but I entered the start_index incorrectly. With start_index = 47, the measured projection matches the information in geometry.txt. Can these parameters be applied to LEAP's "modular-beam" geometry?

[kylechampley]

Oh, I forgot. I had to modify the detector shift amounts. I used the following: shift_distance_col = -47.138 # mm shift_distance_row = -38.139 # mm

I am also concerned that you are still using an older version of LEAP. To check the LEAP version run the following command: leapct.about() It should say version 1.11

[ProjX0]

The versions were different. Everything is working fine now that I matched the version to 1.11. Thank you very much. By the way, could you please explain how you estimated shift_distance_col and shift_distance_row?

[kylechampley]

Good question. If your data were NOT truncated, I would have just used this algorithm: find_centerCol

Unfortunately for offset scan data this algorithm does not work properly yet. For centerRow, I just looked at the projections and saw something that looked like a rotation stage at the bottom of the images. It looked like the projections were tangent to this flat surface, so I just flipped the sign of the value you gave me which approximately put this parameter in the right spot. Good image quality is achievable if this parameter is close to the true value, but it does not have to be perfect.

For centerCol, I estimated this parameter by performing several single-slice reconstructions with different values. Then I just choose the one that looks best. For this strategy, I just added a new algorithm in LEAP (again, in the champley_dev branch) to do this automatically. Here is how to run it: from preprocessing_algorithms import * centerCol = leapct.get_centerCol() centerCols = np.array(range(-5,5+1))+centerCol print(centerCols) f_stack = parameter_sweep(leapct, g, 'centerCol', centerCols) leapct.display(f_stack)

The preprocessing_algorithms.py file is in the utils directory.

[kylechampley]

Oh, wait, I just realized that centerCol and centerRow can be calculated from your geometry file like this: centerCol = p0.x/pitch centerRow = numRows - 1 - p0.y/pitch

Once I used these parameters, I got even better results.

[ProjX0]

thank you very much. It's working well now. However, I discovered one more unexpected result during my attempt.

In binning mode, when using leapct.set_default_volume(), the reconstruction works well. It also works fine when using leapct.set_volume(750, 750, 450, 0.2, 0.2).

However, when not using binning mode (using full data), the images are different when using leapct.set_default_volume() and leapct.set_volume(750, 750, 450, 0.2, 0.2).

When using leapct.set_default_volume(), the image comes out well, but when using leapct.set_volume, it looks as if the FBP filter is not applied.

The whole code is as follows:

[kylechampley]

I actually got a cudaMalloc error when I tried this. Turns out this was caused by LEAP trying to allocate data with a negative number of elements. I fixed this bug and pushed the changes.

OK, so the real issue here is that you need to do the following: leapct.set_volume(750, 750, 450, 0.2, 0.2, 0.0, 0.0, -47.274899)

Note that the z=0 plane is defined by the centerRow parameter. Since your centerRow parameter is far from numRows/2, your field of view is not centered around zero, it is centered around z=-47.274899, so you need to shift your volume to be centered in the field of view. One way to automate this is the following leapct.set_default_volume() leapct.set_volume(750, 750, 450, 0.2, 0.2, 0.0, 0.0, leapct.get_offsetZ())

And one more thing. To get rid of that bright ring around your reconstruction, add the following line to your script: leapct.set_truncatedScan(True)

[kylechampley]

Thanks @ProjX0 and @hws203 for this discussion. It has helped me find and fix several bugs in LEAP, making this toolbox better. If you find LEAP useful, please consider giving it a "Star" by clicking the star icon on the LEAP main page. Stars increase the visibility of LEAP, which brings in more users, which increases the chance of people reporting bugs such as yourselves, which helps improve LEAP for myself and everyone else who uses it.

[ProjX0]

Of course. I clicked "Star" because it is such a useful toolkit.

You mentioned two methods regarding the offset Z as follows:

1. leapct.set_volume(750, 750, 450, 0.2, 0.2, 0.0, 0.0, -47.274899)
2. leapct.set_default_volume() leapct.set_volume(750, 750, 450, 0.2, 0.2, 0.0, 0.0, leapct.get_offsetZ())

I have entered these(1 or 2) into the code, but the output still appears blurry, as shown in the image below.

I changed only the offset Z setting in the whole code that I shared just before, as you suggested.

P.S. The method to solve the truncation artifact, leapct.set_truncatedScan(True), works well. Thank you.

[kylechampley]

Are you using the latest code? I just fixed this bug about 48 hours ago. You can just pull from the main branch because I merged all changes to main. Also I recommend cropping the data like this:

if bin_mode == True: g = leapct.crop_projections([3,numRows-4], [3,numCols-4], g) else: g = leapct.crop_projections([6,numRows-7], [6,numCols-7], g)

[hws203]

I have checked your v1.11 version about the offscan issue. And I could see that there is still a FOV limit issue at the case of left-side offscan with the modification of "centerCols = 0.5*(numCols-1)-100" .

As you recommend if I add the leapct.set_diameterFOV(leapct.get_voxelWidth()*leapct.get_numX()), then there is no such issue.

But I hope that it will work as like righ-side offscan example.

I added star mark at your main page too.

Best regards.

[ProjX0]

@kylechampley Yes, I am using the latest code. Just in case, I downloaded the main code again and reinstalled it. I also tried cropping the data as you suggested, but the issue is still occurring. (It only occurs when bin_mode is false)

The whole code is as follows:

[ProjX0]

Another issue is when using iterative reconstruction techniques in offset scans, the resulting image shows artifacts in the center as shown below. Is there a way to resolve this artifact?

[kylechampley]

Thanks for your patience. I found and fixed another bug that hopefully solves your issue. Apparently, I forgot to account for the extra memory used by texture memory and I was running out of memory. Now I account for this and the data is chunking into smaller pieces. The change just got pushed to the champley_dev branch; it should say version 1.12. Please check this out and let me know if the issue is resolved. By the way, what GPU are you using?

Now let's talk about iterative reconstruction. That ring appears there because you are not converged; you need more iterations. I can tell because your result looks blurry. If you start with an FBP reconstruction this ring artifact won't be so strong because it'll be closer to convergence.

Also notice that the ring is the boundary of the edge of your detector as projected into the volume. Inside the ring you have 360 degrees of projections and outside you have less than 360 degrees. One way to mitigate these artifacts quicker with iterative reconstruction is to use an algorithm that accounts for this variable range in angular coverage. SART, SIRT, and ASDPOCS due this naturally. If you want to use LS, WLS, RLS, or RWLS make sure you add this argument: preconditioner='SQS'

Finally, I would like to remark about other issues with iterative reconstruction. Iterative reconstruction can outperform FBP only when you have an accurate model. If your geometry specification is off by too much or you have a lot of beam hardening, scatter, etc., iterative reconstruction might actually give a poorer result due to overfitting. I would focus on getting all your geometry specifications into LEAP before doing iterative reconstruction. You will also have issues with iterative reconstruction because of truncation artifacts. To avoid these, make sure your reconstruction is reconstructing the whole object. This includes turning off the windowing like this: leapct.set_diameterFOV(2.0leapct.get_voxelWidth()leapct.get_numX())

Notice I added a 2.0 in the above calculation which will remove all the windowing.

Good luck!

[hws203]

I have checked the version 1.12 which included in charmley-dev branch, and I can see the FOV-limit issue of off-scan has been cleared now. My GPU is RTX-3090TI with 24GB memory.

[ProjX0]

@kylechampley First issue I am using an RTX 4090 with a single GPU. In version 1.12, I tested the FBP reconstruction with "bin_mode = false", and the results are shown below.

The image blurring still occurred with the specific FOV settings, but when the FOV was set as on the right, the image blurring did not occur.

Second issue I tested the recommendations you gave previously for reducing ring artifacts and achieving better iterative reconstruction results. The main modifications and additions to the code are as follows.

bin_mode = True ... leapct.set_diameterFOV(2.0leapct.get_voxelWidth()leapct.get_numX())

f = leapct.FBP(g) filters = filterSequence(2e2) filters.append(TV(leapct, delta=0.01 / 100.0, p=1.0)) filters.append(histogramSparsity(leapct, mus=[0.0, 0.02], weight=0.00001)) leapct.RWLS(g, f, iteration_number, filters, preconditioner='SQS')

I used the RWLS reconstruction method and compared the number of iterations at 5, 10, and 50.

As you can see from the results, the ring artifact inside became more pronounced as the number of iterations increased. And when using other reconstruction methods (SART, ASD-POCS), the ring artifacts are even more pronounced compared to the RWLS method.

I think this seems to be a very normal tendency as you mentioned above. So, rather than increasing the number of iterations related to reconstruction, it might be necessary to use another process to reduce the ring artifacts.

[hws203]

@kylechampley. I found the almost same issue which named as "ramp-filter missing issue" by me. As you can see below image, after doFBP, the recon slice shows very blurred one as like there is no ramp-filer adapted during back-propagation.

And I attached my test script too. For your side test, I shared my image files(https://drive.google.com/file/d/1KvSFzIEzhvaSmuOSkTGAsPNymWqwTm0I/view?usp=sharing). Please check this case. From Jeff.

[hws203]

[kylechampley]

I just cannot reproduce those blurry results that you guys are getting. I've tried both of your data sets on two different computers (Windows with 2 1080Ti GPUs and Linux with 2 4090 GPUs). All I can say is run the leapct.about() command and make sure it says version 1.12. I know you said you checked this, but I am at a loss as to the cause of those blurry results.

@hws203 I don't understand why you included the modular-beam stuff in that script. I skipped all of that. The reconstruction didn't look great, but it definitely didn't look blurry like what you posted. There is likely something wrong with the geometry specification. Anyway, here is what I got.

@ProjX0 I would not include the histogram sparsity regularizer. This filter is only for severely ill-posed reconstructions like limited angle or few-view. In addition, the values you give it do not agree with the LAC values in your object, so it will likely give very strange results with the parameters you specified. If I were you, I would focus on getting all your geometry parameters correct before I started working on iterative reconstruction. From the looks of your geometry file, your system has a fair amount of variability from projection to projection.

Oh, and be careful when you call them "ring artifacts". I agree that the artifact is a ring, but "ring artifacts" is a specific term in CT that is caused by detector pixel-to-pixel gain variations that cause several ring artifacts in reconstructions.

[ProjX0]

Thank you so much for the detailed and great feedback.

As you pointed out in https://github.com/LLNL/LEAP/issues/35#issuecomment-2115606743, when centerCol and centerRow are not the center, where should I make changes in the following modularbeam code?

T_phi = 2.0np.pi/float(numAngles) for n in range(numAngles): phi = nT_phi-0.5*np.pi

sourcePositions[n,0] = sod*np.cos(phi)
sourcePositions[n,1] = sod*np.sin(phi)

moduleCenters[n,0] = (sod-sdd)*np.cos(phi)
moduleCenters[n,1] = (sod-sdd)*np.sin(phi)

rowVectors[n,2] = 1.0

colVectors[n,0] = -np.sin(phi)
colVectors[n,1] = np.cos(phi)

leapct.set_modularbeam(numAngles, numRows, numCols, pixelSize, pixelSize, sourcePositions, moduleCenters, rowVectors, colVectors)

[kylechampley]

I recommend starting with a cone-beam specification, converting it to modular-beam, and then tweaking the geometry specification to match your geometry. Like this:

leapct.set_conebeam(...) leapct.convert_to_modularbeam() ...

See demo 25 for an example.

You can also utilize scipy for making your rotation matrices. For example: from scipy.spatial.transform import Rotation as R A = R.from_euler("xyz", [1.0, 2.0, 3.0], degrees=True).as_matrix()

[hws203]

I found the key function for making this blurry result, that is the custom volume setting function which is 'leapct.set_volume()'. If I use leapct.set_default_volume() function, then there is no such blurry result. I can not understand why this happens with leapct.set_volume() function. Please explain this reason.

Below is the reconed slice with leapct.set_default_volume(). Below is my new script.

[ProjX0]

@hws203 In my cases, I've encountered an issue with my geometry and data where the image blurs when using leapct.set_default_volume(). Despite asking @kylechampley several times, we haven't been able to identify the problem. I'm still at a loss about what's causing this.

[hws203]

@ProjX0 I guess that this blurry issue is also related to the allocated memory size of GPU. But I can not identify it clearly. But at my case study, there might be some clue for this issue.

[kylechampley]

The problem with this bug is that I cannot reproduce it. I have tried on two different computers and on both computers I get the expected result, no matter how I set the volume.

There is nothing magical about the set_default_volume function. For cone-beam it sets the voxel size like this: voxelWidth = sod / sdd pixelWidth voxelHeight = sod / sdd pixelHeight and then it sets the number of voxels and offsetZ so that the volume fills the field of view and is centered in the field of view.

The likely issue is that there is some failure in the CUDA functions that do the filtering and it fails silently (i.e., no error message is printed to the screen). The failure is likely a malloc error which would explain why it seems to depend on the size of the volume.

I will work on catching all possible errors in the filtering functions of LEAP and push a new version in about 12 hours from now. I'll let you know when it is ready.

[kylechampley]

OK, I pushed changes to the champley_dev branch. Please download, compile, and try using it. Be sure to run the command leapct.about() and make sure it says it was compiled today.

This new version catches all CUDA errors and prints the error message to the screen if something went wrong. Please watch your screen for any error messages. I also added a logging utility to LEAP. I will be adding more status logs in the future, but currently I only included some. To see log messages, run this command: leapct.set_log_debug().

If you encounter an error, please send me all LEAP logging and error statements printed to the screen. You can send them directly to my email address if you wish.

[hws203]

I checked your logged version from champley-branch, but unfortunately, there is no error message. I will check it more Today.

[hws203]

I got some clue after changing the leapct.set_volume() to leapct.set_default_volume(), there are some logging about slice generation at GPU. So, I guess that at the blurry result case, your code does not reach this line code. Please check this point and add more log if you need. And let me try to check it again. And if I change the volume size by 2404x2404x2196 at leapct.set_volume(), then there is no blurry issue too.

[hws203]

@kylechampley leap.find_centerCol() function is also affected by this volume size, this find_centerCol() function also works correctly only at the case of default_volume size. Please check this issue too.

[kylechampley]

The volume specification (e.g., volume size) has absolutely no impact on the output of leapct.find_center().

But the output of leapct.find_center() will affect leapct.set_default_volume(). And it definitely should. The placement of your detector affects the size of the field of view that you can reconstruction. Thus, different values of centerCol will make the volume specified by leapct.set_default_volume() have different values.

This is all intentional.

[ProjX0]

Thank you once again for your kind and prompt response.

I conducted tests on three cases. In cases (1) and (3), when I set the volume and FOV as per the code inserted in the figure, there was no image blurring. However, in case (2), when the settings were applied, image blurring occurred, and no log messages were generated. (For reference, the right side of the figure shows the log messages output for each case.)

I have sent the code I used via email.

[hws203]

I found that the cudaFFT library makes this issue. So if I build Leap v1.12 without "__INCLUDE_CUFFT" option, then there is no such blurry issue. I guess that cudaFFT usage is not matched to their intent sometimes.

[kylechampley]

OK, I finally fixed the bug! Again, please pull from the champley_dev branch.

@hws203 I do not recommend turning off CUFFT. You will lose some features in LEAP and things will go slower.

The reason why the ramp filter was not being applied was that I was launching a CUDA kernel with a block size that was too big. I always test LEAP on machines with multiple GPUs and with multiple GPUs the job gets divided across all the GPUs and then it never happened that I used a block size that was too big. To reproduce the effect you guys saw, I had to tell LEAP to use just one GPU.

Also, you guys sometimes did not encounter this bug because LEAP was dividing the job into smaller pieces when LEAP determined that you'd run out of GPU memory if it wasn't split up.

Regardless, I am confident that this is fixed.

[hws203]

Ramp-filter with CudaFFT is now working well at your new champly-dev version.

[ProjX0]

I also did not experience any image blurring with my setup. Thank you.

However, while performing modularbeam using my geometry file, I encountered another issue during testing. I discovered an issue where the code does not stop at leapct.FBP(g, f) (resulting in infinite loading). Upon checking the Task Manager, I noticed that the CUDA operations were not active.

I have sent you detailed information about this problem via email. You can check the most recent email (title: Issue: modularbeam) I sent for detailed information about this problem.

[kylechampley]

Yes, I saw your email. I solved that problem and a few others in the script you sent me, but the reconstruction still looks a bit blurry. I'll respond to you by email in a day or two with my corrections.

I am closing this issue because I believe the offset scan issues have been resolved in LEAP.