CT正投影

所用知识

Beer-Lambert定律；
线性衰减系数：线性衰减系数 µ 表示光子穿过单位路径与物质发生作用的概率。
量纲是长度的倒数（单位 $c{m}^{-1}$ ）。

编写程序，实现特殊几何形状。如圆、椭圆、三角形、正方形等的正投影

夹心圆的正投影

示例图片如下：
$\theta =0$：

$\theta >0$

其中，大圆半径为4，小圆半径为1；大圆（不包括小圆部分）质量衰减系数为0.5；小圆部分质量衰减系数为1；大圆以外部分质量衰减系数为0；
大圆方程为：$$x^2+y^2=16$$
小圆方程为：$$(x{)}^2+(y-2{)}^2=1$$
对每条光线通过物体发生的反应叠加求每个角度下灰度值图像。

数学模型

分段考虑，有两种情况：

1. 光线经过的路径全部为大圆(0.5);
2. 光线经过的路径一部分大圆，一部分小圆（0.5+0.5）。
   若为图所示情况，则数学模型中分为三段：
   $$f(x)=\{\begin{array}{rlr}2\sqrt{16-x^2}\ast 0.5& & x\in (-4,-1)\\ 2\sqrt{16-x^2}\ast 0.5+2\sqrt{1-x^2}\ast 0.5& & x\in (-1,1)\\ 2\sqrt{16-x^2}\ast 0.5& & x\in (1,4)\end{array}$$

探测器旋转

若探测器旋转，则将坐标轴变换。探测器逆时针旋转，就相当于坐标轴顺时针旋转相应的度数。
关于$\theta$函数如下所示：
$$f(x)=\{\begin{array}{rlr}2\sqrt{16-x^2}\ast 0.5& & x\in (-4,\sqrt{2sin\theta }-1)\\ 2\sqrt{16-x^2}\ast 0.5+2\sqrt{1-(\sqrt{2sin\theta }-1{)}^2}\ast 0.5& & x\in (\sqrt{2sin\theta }-1,\sqrt{2sin\theta }+1)\\ 2\sqrt{16-x^2}\ast 0.5& & x\in (\sqrt{2sin\theta }+1,4)\end{array}$$

代码如下：

import math
import struct
import numpy as np
np.set_printoptions(threshold=np.inf)
import matplotlib.pyplot as plt #plt用于显示图片
#import matplotlib.image as mpimg #mpimg用于读取图片
#p为旋转角度
x = np.linspace(-4,4,2000) 

p = np.linspace(0,2*(math.pi),int((math.pi)*1000))

#a = [[0 for i in range(2000)] for i in range(3141)]
a = np.arange(3141*2000, dtype=np.float64).reshape((3141,2000))
b = np.arange(3141*2000, dtype=np.int32).reshape((3141,2000))#np.int32 表示整数；reshape为二维数组

#print (math.sin(p[3100]))
#matplotlib中没有函数可以将RGB图转换为灰度图，所以我们定义一个
def banfill(x):
    grayvalue=10*(math.sqrt(16-np.power(x,2)))
    return grayvalue

def manfill(x,p):
	b=math.sin(p)
	grayvalue1=10*(np.sqrt(16-np.power(x,2))+np.sqrt(1-np.power(2*b-x,2)))
	#print(grayvalue)
	return grayvalue1


def grayresult(x,p):
    a=math.sin(p)
    if (x>-4 and x<(2*a-1)):
        grayv1=banfill(x)
        return grayv1
    elif(x>(2*a-1) and x<(2*a+1)):
        grayv2=manfill(x,p)
        return grayv2
    elif (x>(2*a+1) and x<4):
        grayv=banfill(x)
        return grayv
for i in range(3141):
    for j in range(2000):
        a[i][j]=grayresult(x[j], p[i])
        
        b[i][j]=np.nan_to_num(a[i][j])##将数组中的值化为整数，方便渲染图像

#print(b)	   #检验b是否正确
        
plt.imshow(b,cmap='gray')

plt.show()

程序结果

正方形ct正投影

思路

其中，将x轴分为三部分：
最左边和最右边部分求弦长原理相同，都是使用相似三角形方法；
中间部分无论x值为多少，弦长均相等。

离散图像求正投影

旋转图像（适当扩充），求一条射线经过像素灰度值的和。

import cv2
import matplotlib.pyplot as plt
import numpy as np
import math
from PIL import Image

from time import perf_counter

def sum1(p,image3):
    (h, w) = image3.shape[0:2]
    (cX, cY) = (w/2, h/2)
    M= cv2.getRotationMatrix2D((cX, cY), p, 1.0)
    image1 = cv2.warpAffine(image3, M, (w, h))
    #print(image1.shape)
    m=[0]*724
    for i in range(image1.shape[0]):
        for j in range(image1.shape[1]):
            m[i]=m[i]+image1[i][j]
        m[i]=m[i]/301
    return m

def findmax(x):
    max1=0
    for i in range(len(x)):
        if (x[i]>max1):
            max1=x[i]
    return max1
if __name__ == '__main__':
    t0 = perf_counter()
    img = cv2.imread('lena512gray.png', cv2.IMREAD_GRAYSCALE)
    a=round((724-img.shape[0])/2)
    b=round((722-img.shape[0])/2)
    image = cv2.copyMakeBorder(img,a,a,b,b, cv2.BORDER_CONSTANT,value=[0,0,0])
    theta = np.linspace(0,360,3600)
    grayresult = np.arange(len(theta)*724, dtype=np.float64).reshape((len(theta),724))
    for i in range(len(theta)):
            grayresult[i]=sum1(theta[i],image)
            print(grayresult[i])
    print(grayresult)
    plt.imshow(grayresult,cmap='gray')
    plt.show()
    cv2.imwrite('result.png',grayresult);

输入图像

输出结果