Skip to content

Commit

Permalink
Add more examples.
Browse files Browse the repository at this point in the history
  • Loading branch information
@z0gSh1u
z0gSh1u committed Jul 7, 2024
1 parent b54b543 commit 6b8843e
Show file tree
Hide file tree
Showing 17 changed files with 1,141 additions and 29 deletions.
12 changes: 7 additions & 5 deletions crip/physics.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -255,14 +255,10 @@ def computeMu(atten: Atten, spec: Spectrum, energyConversion: EnergyConversion)
return np.sum(spec.omega * eff * mus) / np.sum(spec.omega * eff)


def computeAttenedSpectrum(spec: Spectrum, attens: Or[Atten, List[Atten]], Ls: Or[float, List[float]]) -> Spectrum:
def computeAttenedSpectrum(spec: Spectrum, attens: List[Atten], Ls: List[float]) -> Spectrum:
''' Calculate the spectrum after attenuation through `attens` with thickness `Ls` [mm]
using Beer-Lambert law.
'''
if isType(attens, Atten):
attens = [attens]
if isType(Ls, float):
Ls = [Ls]
cripAssert(len(attens) == len(Ls), '`attens` should be paired with `Ls`.')

N = len(attens)
Expand All @@ -275,6 +271,12 @@ def computeAttenedSpectrum(spec: Spectrum, attens: Or[Atten, List[Atten]], Ls: O
return Spectrum(omega, spec.unit)


def normalizeSpectrum(spec: Spectrum):
''' Normalize a Spectrum.
'''
return Spectrum(spec.omega / spec.sumOmega, spec.unit)


def forwardProjectWithSpectrum(lengths: List[TwoD], materials: List[Atten], spec: Spectrum, energyConversion: str):
''' Perform forward projection using `spec`. `lengths` is a list of corresponding length [mm] images
(projection or sinogram) of `materials`.
Expand Down
24 changes: 6 additions & 18 deletions crip/spec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,25 +18,13 @@ def decompMaterial(src: Atten, base1: Atten, base2: Atten, mode='coeff') -> NDAr
''' Decompose material `src`'s attenuation onto two materials `base1` and `base2`.
Return decomposing coefficients in `DiagEnergyRange` when `mode='coeff'`, or proportion when `mode='prop'`.
'''
cripAssert(mode in ['coeff', 'prop'], f'Invalid mode: {mode}.')
cripAssert(False, 'The implementation of `decompMaterial` requires more investigation.')

srcMu = src.mu
baseMu1 = base1.mu
baseMu2 = base2.mu
M = np.array([baseMu1, baseMu2], dtype=DefaultFloatDType).T

if mode == 'prop':
M /= srcMu
srcMu = np.ones_like(srcMu)

res = np.linalg.pinv(M) @ (srcMu.T)

return res


def deDecompCoeff(lowSpec: Spectrum, highSpec: Spectrum, base1: Atten, len1: Or[NDArray, Iterable], base2: Atten,
def deDecompProjCoeff(lowSpec: Spectrum, highSpec: Spectrum, base1: Atten, len1: Or[NDArray, Iterable], base2: Atten,
len2: Or[NDArray, Iterable]):
''' Compute the decomposing coefficient (Order 2 with bias term) of two spectra onto two material bases.
''' Compute the decomposing coefficient (Order 2 with bias term) of two spectra onto two material bases
used in the projection domain.
'''

def computePostlog(spec, attens, Ls):
Expand All @@ -63,7 +51,7 @@ def computePostlog(spec, attens, Ls):
def deDecompProj(lowProj: TwoOrThreeD, highProj: TwoOrThreeD, coeff1: NDArray,
coeff2: NDArray) -> Tuple[TwoOrThreeD, TwoOrThreeD]:
''' Perform dual-energy decompose in projection domain point-by-point using coeffs.
Coefficients can be generated using @see `deDecompCoeff`.
Coefficients can be generated using @see `deDecompProjCoeff`.
'''
cripAssert(isOfSameShape(lowProj, highProj), 'Two projection sets should have same shape.')
cripAssert(
Expand Down Expand Up @@ -194,7 +182,7 @@ def compose3(b1, b2, b3, v1, v2, v3):


def vmi2Mat(b1, b2, b1Mat: Atten, b2Mat: Atten, E: int):
''' Virtual Monoenergetic Imaging using two-material decomposition.
''' Virtual Monoenergetic Imaging using two-material decomposition at energy `E` [keV].
'''
return b1 * b1Mat.mu[E] + b2 * b2Mat.mu[E]

Expand Down
97 changes: 97 additions & 0 deletions example/.ipynb_checkpoints/01-LAC-of-Al-and-Cu-checkpoint.ipynb
View file

Large diffs are not rendered by default.

169 changes: 169 additions & 0 deletions example/.ipynb_checkpoints/02-Compute-LAC-checkpoint.ipynb
View file

Large diffs are not rendered by default.

105 changes: 105 additions & 0 deletions example/.ipynb_checkpoints/03-DECT-checkpoint.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "e280d524-c2ae-4546-ba08-569b23d3e6b3",
"metadata": {},
"outputs": [],
"source": [
"# This notebook is about crip and Dual-Energy CT."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "7f5ba35c-cf14-4f66-a6b7-7d7702eaf2cf",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from crip.physics import Atten, Spectrum\n",
"from crip.spec import decompMaterial, deDecompProjCoeff"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "3225da4b-d248-4e2e-9449-c91111b28daf",
"metadata": {},
"outputs": [],
"source": [
"spec80 = Spectrum.fromFile('03-Spectrum-80.txt', 'keV')\n",
"spec140 = Spectrum.fromFile('03-Spectrum-140.txt', 'keV')\n",
"PMMA = Atten.fromBuiltIn('PMMA')\n",
"Al = Atten.fromBuiltIn('Al')"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "20481870-1b22-47ff-a8b3-3f6f18edc443",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"c:\\Users\\13637\\anaconda3\\envs\\main\\lib\\site-packages\\crip\\spec.py:34: RuntimeWarning: divide by zero encountered in log\n",
" return -np.log(attenedSpec.omega / flat)\n"
]
},
{
"ename": "CripException",
"evalue": "Inputs should be 1D sequence.",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mCripException\u001b[0m Traceback (most recent call last)",
"Cell \u001b[1;32mIn[9], line 4\u001b[0m\n\u001b[0;32m 2\u001b[0m L1 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m20\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m2\u001b[39m) \u001b[38;5;241m*\u001b[39m \u001b[38;5;241m10\u001b[39m \u001b[38;5;66;03m# lengths to compute on\u001b[39;00m\n\u001b[0;32m 3\u001b[0m L2 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m5\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m) \u001b[38;5;241m*\u001b[39m \u001b[38;5;241m10\u001b[39m\n\u001b[1;32m----> 4\u001b[0m c1, c2 \u001b[38;5;241m=\u001b[39m \u001b[43mdeDecompProjCoeff\u001b[49m\u001b[43m(\u001b[49m\u001b[43mspec80\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mspec140\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mPMMA\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mL1\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mAl\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mL2\u001b[49m\u001b[43m)\u001b[49m\n",
"File \u001b[1;32mc:\\Users\\13637\\anaconda3\\envs\\main\\lib\\site-packages\\crip\\spec.py:45\u001b[0m, in \u001b[0;36mdeDecompProjCoeff\u001b[1;34m(lowSpec, highSpec, base1, len1, base2, len2)\u001b[0m\n\u001b[0;32m 42\u001b[0m postlogLow\u001b[38;5;241m.\u001b[39mappend(computePostlog(lowSpec, [base1, base2], [i, j]))\n\u001b[0;32m 43\u001b[0m postlogHigh\u001b[38;5;241m.\u001b[39mappend(computePostlog(highSpec, [base1, base2], [i, j]))\n\u001b[1;32m---> 45\u001b[0m beta, gamma \u001b[38;5;241m=\u001b[39m \u001b[43mfitPolyV2D2\u001b[49m\u001b[43m(\u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpostlogLow\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpostlogHigh\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mlenCombo\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241m.\u001b[39mT\n\u001b[0;32m 47\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m beta, gamma\n",
"File \u001b[1;32mc:\\Users\\13637\\anaconda3\\envs\\main\\lib\\site-packages\\crip\\shared.py:235\u001b[0m, in \u001b[0;36mfitPolyV2D2\u001b[1;34m(x1, x2, y, bias)\u001b[0m\n\u001b[0;32m 229\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfitPolyV2D2\u001b[39m(x1: NDArray, x2: NDArray, y: NDArray, bias: \u001b[38;5;28mbool\u001b[39m \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m NDArray:\n\u001b[0;32m 230\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m''' Fit a degree-2 polynomial using pseudo-inverse to a pair of variables `x1, x2`.\u001b[39;00m\n\u001b[0;32m 231\u001b[0m \u001b[38;5;124;03m Output 6 coefficients `c[0~5]`, minimizing the error of `y` and\u001b[39;00m\n\u001b[0;32m 232\u001b[0m \u001b[38;5;124;03m `c[0] * x1**2 + c[1] * x2**2 + c[2] * x1 * x2 + c[3] * x1 + c[4] * x2 + c[5]`.\u001b[39;00m\n\u001b[0;32m 233\u001b[0m \u001b[38;5;124;03m If `bias` is False, `c[5]` will be 0.\u001b[39;00m\n\u001b[0;32m 234\u001b[0m \u001b[38;5;124;03m '''\u001b[39;00m\n\u001b[1;32m--> 235\u001b[0m \u001b[43mcripAssert\u001b[49m\u001b[43m(\u001b[49m\u001b[43mis1D\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx1\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01mand\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mis1D\u001b[49m\u001b[43m(\u001b[49m\u001b[43mx2\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;129;43;01mand\u001b[39;49;00m\u001b[43m \u001b[49m\u001b[43mis1D\u001b[49m\u001b[43m(\u001b[49m\u001b[43my\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mInputs should be 1D sequence.\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m 236\u001b[0m cripAssert(isOfSameShape(x1, x2) \u001b[38;5;129;01mand\u001b[39;00m isOfSameShape(x1, y), \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m`x1`, `x2` and `y` should have same shape.\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m 238\u001b[0m x1sq \u001b[38;5;241m=\u001b[39m x1\u001b[38;5;241m.\u001b[39mT\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\n",
"File \u001b[1;32mc:\\Users\\13637\\anaconda3\\envs\\main\\lib\\site-packages\\crip\\utils.py:40\u001b[0m, in \u001b[0;36mcripAssert\u001b[1;34m(cond, hint)\u001b[0m\n\u001b[0;32m 37\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m''' The only assert method for crip.\u001b[39;00m\n\u001b[0;32m 38\u001b[0m \u001b[38;5;124;03m'''\u001b[39;00m\n\u001b[0;32m 39\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m cond:\n\u001b[1;32m---> 40\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m CripException(hint)\n",
"\u001b[1;31mCripException\u001b[0m: Inputs should be 1D sequence."
]
}
],
"source": [
"# Get the decomposing coefficients in the projection domain onto two basis materials.\n",
"L1 = np.arange(0, 20 + 1, 2) * 10 # lengths to compute on\n",
"L2 = np.arange(0, 5 + 1, 1) * 10\n",
"c1, c2 = deDecompProjCoeff(spec80, spec140, PMMA, L1, Al, L2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2fe9c9bc-61bd-433d-a398-69bfbadea993",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
97 changes: 97 additions & 0 deletions example/01-LAC-of-Al-and-Cu.ipynb
View file

Large diffs are not rendered by default.

169 changes: 169 additions & 0 deletions example/02-Compute-LAC.ipynb
View file

Large diffs are not rendered by default.

132 changes: 132 additions & 0 deletions example/02-Spectrum.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
Energy(keV) Omega
10 0
11 0
12 45.289
13 613.43
14 4601.5
15 20979
16 68137
17 1.6875e+05
18 3.4126e+05
19 5.9946e+05
20 9.3977e+05
21 1.3477e+06
22 1.8016e+06
23 2.279e+06
24 2.7575e+06
25 3.1929e+06
26 3.6283e+06
27 4.0158e+06
28 4.3626e+06
29 4.6606e+06
30 4.8893e+06
31 5.0867e+06
32 5.251e+06
33 5.3719e+06
34 5.4585e+06
35 5.5074e+06
36 5.5308e+06
37 5.5326e+06
38 5.5145e+06
39 5.4751e+06
40 5.4248e+06
41 5.3624e+06
42 5.2848e+06
43 5.2017e+06
44 5.1133e+06
45 5.0173e+06
46 4.9213e+06
47 4.8186e+06
48 4.7132e+06
49 4.6087e+06
50 4.5009e+06
51 4.3939e+06
52 4.2851e+06
53 4.1813e+06
54 4.0759e+06
55 3.9712e+06
56 3.8667e+06
57 3.7647e+06
58 1.9414e+07
59 3.5678e+06
60 3.4702e+06
61 3.3747e+06
62 3.2823e+06
63 3.1906e+06
64 3.1032e+06
65 3.0146e+06
66 2.93e+06
67 1.1982e+07
68 2.7642e+06
69 5.0691e+06
70 2.4438e+06
71 2.3776e+06
72 2.314e+06
73 2.2517e+06
74 2.1899e+06
75 2.1296e+06
76 2.0681e+06
77 2.0102e+06
78 1.9527e+06
79 1.8953e+06
80 1.839e+06
81 1.7836e+06
82 1.7312e+06
83 1.6778e+06
84 1.6252e+06
85 1.5733e+06
86 1.5231e+06
87 1.4741e+06
88 1.4253e+06
89 1.3772e+06
90 1.3299e+06
91 1.2832e+06
92 1.238e+06
93 1.1924e+06
94 1.1481e+06
95 1.1047e+06
96 1.0618e+06
97 1.0188e+06
98 9.7683e+05
99 9.3532e+05
100 8.948e+05
101 8.5467e+05
102 8.1398e+05
103 7.7442e+05
104 7.3478e+05
105 6.955e+05
106 6.5686e+05
107 6.1755e+05
108 5.7906e+05
109 5.3978e+05
110 5.014e+05
111 4.6263e+05
112 4.2349e+05
113 3.8472e+05
114 3.4154e+05
115 3.0002e+05
116 2.5884e+05
117 2.1824e+05
118 1.6157e+05
119 79522
120 0
121 0
122 0
123 0
124 0
125 0
126 0
127 0
128 0
129 0
130 0
131 0
132 0
133 0
134 0
135 0
136 0
137 0
138 0
139 0
140 0
Loading

0 comments on commit 6b8843e

Please sign in to comment.