-
Notifications
You must be signed in to change notification settings - Fork 0
/
radon.py
128 lines (109 loc) · 5.61 KB
/
radon.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
import numpy as np
iterationCount = 10
def createBackProjectionFilter(len):
if len%2==0: raise NameError("Incorrect len")
out = np.zeros(len)
middle=int(len/2)
out[middle]=1
for i in range(1,middle,2):
val=-4.35/(np.power(np.pi,2) * np.power(i,2))
out[middle-i]=val
out[middle+i] = val
return out
def radonTransform(input, stepSize=1, stepsArray=None, detectorsNumber=100, detectorsWidth=140, output=None, normalize=True):
if stepsArray is None: stepsArray = np.arange(0,180,stepSize)
xSize = int(180/stepSize+1)
if output is None: output = np.zeros((detectorsNumber,xSize))
circleRadius = np.sqrt(np.power(len(input)/2,2)+np.power(len(input[0])/2,2) )
center = (len(input[0])/2,len(input)/2)
output = radonCircleLoop(input,output, stepsArray, stepSize, center, circleRadius,detectorsNumber,detectorsWidth)
if normalize: output /= max(output.flatten()) #Normalizacja
return output
def inverseRadonTransform(input, stepSize=1, stepsArray=None, detectorsWidth=140, output=None, outputWidth=None, outputHeight=None, normalize=True):
if stepsArray is None: stepsArray = np.arange(0,180, stepSize)
if output is None:
if outputHeight is None: outputHeight = len(input)
if outputWidth is None: outputWidth = outputHeight
output = np.zeros((outputHeight,outputWidth))
circleRadius = np.sqrt(np.power(outputWidth/2,2)+np.power(outputHeight/2,2) )
center = (outputWidth/2,outputHeight/2)
output = radonCircleLoop(input,output, stepsArray, stepSize, center,circleRadius,len(input),detectorsWidth,inverse=True)
if normalize:
for i in range(len(output)):
output[i] /= max(output[i].flatten())
return output
def radonCircleLoop(input, output, stepsArray, step, center, circleRadius, detectorsNumber, detectorsWidth, inverse=False):
detectorDistance = (circleRadius * 2 * detectorsWidth / 180) / detectorsNumber
iterations = [0]
imagesIter = []
if inverse == True:
count = len(stepsArray)
nextIter = np.floor(count/iterationCount)
for i in range(1,iterationCount-1):
iterations.append(i*nextIter)
iterations.append(count-1)
for index,stepAngle in enumerate(stepsArray):
centralEmiterPos = (center[0] + circleRadius * np.sin(np.radians(stepAngle)),center[1] + np.cos(np.radians(stepAngle)) * circleRadius)
centralReceiverPos = (center[0] - circleRadius * np.sin(np.radians(stepAngle)),center[1] - np.cos(np.radians(stepAngle)) * circleRadius)
for currentDetector in range(0,detectorsNumber):
distanceFromMainDetector = (currentDetector - (detectorsNumber / 2)) * detectorDistance
cos = np.cos(np.radians(stepAngle))
sin = np.sin(np.radians(stepAngle))
emiterPos = centralEmiterPos[0] + distanceFromMainDetector * cos, centralEmiterPos[1] - distanceFromMainDetector * sin
receiverPos = centralReceiverPos[0] + distanceFromMainDetector * cos, centralReceiverPos[1] - distanceFromMainDetector * sin
if not inverse:
points = BresenhamAlgorithm(input, emiterPos, receiverPos)
if len(points) > 0: output[currentDetector][int(stepAngle/step)] = sum(points) # Normalizacja
else:
color = input[currentDetector, int(stepAngle/step)]
output = BresenhamAlgorithm(input, emiterPos, receiverPos, output, returnOrDraw=False, lineColor=color)
if inverse == True:
if index in iterations:
imagesIter.append(output.copy())
if inverse == True:
return imagesIter
return output
def filterSinogram(input):
lines = len(input[0])
s = 9
while s >= len(input): s -= 2
mask = createBackProjectionFilter(s)
for num in range(0,lines):
input[:, num] = filter1D(input[:,num], mask=mask)
return input
def filter1D(input, mask=None):
output = np.zeros(len(input))
if mask is None: mask = createBackProjectionFilter(9)
maskSize, inputSize = len(mask), len(input)
for X in range(0,inputSize):
for maskX in range(-int(maskSize / 2), int(maskSize / 2) + 1):
cX = X + maskX
if cX < 0: cX += inputSize
if cX >= inputSize: cX -= inputSize
output[X] += input[cX] * mask[maskX + int(maskSize / 2)]
return output
# returnOrDraw : Return list of values if True ; Draw line if False
def BresenhamAlgorithm(input, A, B, output=None, moreThanZeroValues=True, returnOrDraw=True, lineColor=0.5):
if returnOrDraw and output is None: output = []
if not returnOrDraw:
if output is None: raise NameError("output must be given")
outSizeX, outSizeY = len(output[0]), len(output)
inputSizeX, inputSizeY, = len(input[0]), len(input)
X, Y = int(A[0]), int(A[1])
X2, Y2 = int(B[0]), int(B[1])
dx, dy = abs(X-X2), abs(Y-Y2)
xAdd, yAdd = 1 if X<X2 else -1, 1 if Y<Y2 else -1
def bresenhamLoop(X,Y,output):
if returnOrDraw and X >= 0 and Y >= 0 and X < inputSizeX and Y < inputSizeY:
color = input[inputSizeY - 1 - int(Y)][int(X)]
if not moreThanZeroValues or color > 0: output.append(color)
if not returnOrDraw and X>=0 and Y>=0 and Y<outSizeY and X < outSizeX:
output[outSizeY - 1 - int(Y)][int(X)] = max(0, output[outSizeY - 1 - int(Y)][int(X)]+lineColor)
return X+xAdd,Y+yAdd,output
if dx >= dy :
yAdd = float(abs(Y-Y2))/abs(X-X2)*yAdd
while X != X2: X,Y,output = bresenhamLoop(X,Y,output)
else:
xAdd = float(abs(X-X2))/abs(Y-Y2)*xAdd
while Y != Y2: X,Y,output = bresenhamLoop(X,Y,output)
return output