Repository navigation
Fix per-angle geometry indexing in iterative algorithms (MATLAB) and blocksize-1 angle blocks (Python) - #801
Merged
Merged
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
When using non-circular trajectories (e.g., per-angle
DSO,DSD,offDetector,rotDetector, such as on Loop-X), several iterative algorithms apply the geometry or weights (W) of a different angle than the projection currently being updated. Because these values are constant across all angles in circular scans, the bug was previously invisible. However, with per-angle geometry, this mismatch causes the reconstruction to diverge or turn into noise.In addition, this PR fixes a bug in the Python blocksize-1 algorithms where three angles were being projected per update instead of one.
Python Changes
File:
tigre/algorithms/iterative_recon_alg.py(specificallyset_angle_index)blocksize=1,order_subsetsreturns each block as a 1D(3,)row[alpha, theta, psi]with a scalar index. Since commita7c1f8a("fix TypeError: only 0-dimensional arrays can be converted to Python scalars #741") removed the blocksize-1 reshaping inart_data_minimizing, this row is passed directly toAx/Atb.check_georeads this 1D array as a list of three separate rotation angles. As a result, every SART-type update back-projectsW·(p_alpha − A_0°x)twice at angle 0°.set_valso buildsVfrom the same three back-projections.(k,3)and every index block is 1D(k,). Forblocksize > 1, this does nothing as the shapes are already correct.sart,sart_tv,asd_pocs,awasd_pocs,piccs,pcsd, andaw_pcsd.sirt,fista,cgls,lsqr,lsmr, andmlem) yield bit-identical results.DSOstays(1,)since commit26d8e2e. (Checked with NumPy's "ndim>0 to scalar" deprecation turned into an error).MATLAB Changes
1. Corrected Angle Indexing in Loops
In non-OS loops,
jjis the position in the shuffled order, whileindex_angles(jj)is the actual angle. In OS loops,jjis the block number, andorig_index{jj}holds the block's angles. All of these algorithms default toOrderStrategy='random', making correct indexing critical.PICCS.m,SART_TV.m,B_ASD_POCS_beta.m: Updatedgeo.DSD=DSD(jj)/geo.DSO=DSO(jj)to useDSD(index_angles(:,jj))/DSO(index_angles(:,jj)). This now correctly matches the indexing behavior foroffOrigin,offDetector, androtDetector.SART_TV.m: The update now usesW(:,:,index_angles(:,jj))instead ofW(:,:,jj). Also,resL2previously comparedproj(:,:,index_angles)(shuffled) withAx(res,geo,angles)(input order), resulting in a wrong residual that triggered the "error increased" early exit. It now correctly usesprojin input order.OS_AwASD_POCS.m: ChangedDSD(jj)/DSO(jj)toDSD(:,orig_index{jj})/DSO(:,orig_index{jj}). This also fixesAwASD_POCS(which wraps this file).OS_PCSD.m,OS_AwPCSD.m:offOrigin,offDetector,rotDetector,DSD,DSO, andWwere incorrectly accessed by block number (jj); they now properly useorig_index{jj}. This also fixes thePCSDandAwPCSDwrappers.2. Geometry Restoration Before Full-Data
AxCallsSeveral fields were previously left unrestored before error/stop checks,
est_proj,resL2, andBregmansteps, causing these global calls to mistakenly use a single block's geometry.OS_SART.m,OSEM.m,SART_TV.m,B_ASD_POCS_beta.m: Added missinggeo.DSO = DSO.OS_PCSD.m,OS_AwPCSD.m: Addedgeo.DSD,geo.DSO, andgeo.rotDetectorto both restore blocks.OS_SART/OSEM: RequestingresL2used to throw an error (geo.DSD Should be 1x1 or 1xsize(angles,2)). BecauseSART(a wrapper) always requests it, it was computing an incorrect residual and could stop early. This is now resolved.3. Fixed Existing Crashes/Typos
OS_PCSD.m,OS_AwPCSD.m:f0was being used indp_vec=(f-f0)before being initialized. Setf0=fat the start of each iteration (matchingASD_POCS). Also updatedresL2(ii)toresL2(iter).OS_AwPCSD.m: Correctednitertomaxiter.B_ASD_POCS_beta.m: The defaultepsiloncall toim3Dnormwas missing the'L2'argument.