Commit 79ebc609 authored by Edward Andò's avatar Edward Andò
Browse files

rough psi-scan working example

parent 125f41fe
Pipeline #68614 failed with stage
in 1 minute and 3 seconds
"""
EA and BM sooo early in the morning of 2021-06-08
Now that you've run tomopackTo3DguessExperimental-psiScan.py you should have nice 3D positions in psi-Scan.npy
Let's make sure they're OK and then try to use the sample to calibrate the attenuation curve
"""
import radioSphere
import numpy
import tifffile
import matplotlib.pyplot as plt
pos3D = numpy.load("psi-Scan.npy")
radiusMM = 1
rootDir = "./data/2021-05-19-EPFL/"
SDD = 403.784
SOD = 36.3206
zoomPosCentre = -22.90
IO_for_radio = tifffile.imread(f"{rootDir}/01-holder-AVG16/Proj/00025.tif")
radio = tifffile.imread( f"{rootDir}/02-sample-AVG16/Proj/00025.tif") / IO_for_radio
logRadio = -numpy.log(radio)
projMM = radioSphere.projectSphereMM(pos3D, numpy.ones(pos3D.shape[0]), sourceDetectorDistMM=SDD, pixelSizeMM=0.139*4, detectorResolution=[524, 428])
plt.subplot(1, 2, 1)
plt.imshow(logRadio/logRadio.max())
plt.title("Does this measurement...")
plt.subplot(1, 2, 2)
plt.imshow(projMM/projMM.max())
plt.title("look like this identification? (if not Ctrl+C)")
plt.show()
#tifffile.imsave("psi-Scan.tif", im.astype('<f4'))
"""
Now try to use sample for calibration
"""
plt.plot(logRadio.ravel(), projMM.ravel(), 'x')
plt.xlabel("-log(I/I0)")
plt.ylabel("mm")
plt.show()
......@@ -19,7 +19,7 @@ import glob
radioCrop = 0
nSphere = 25
nSphere = 26
radiusMM = 1
......
......@@ -403,38 +403,37 @@ def psiSeriesScanTo3DPositions(radio, # not (necessarily) MM
#tifffile.imsave(cachePsiFile, psiXseries.astype('<f4'))
#print("done.")
L_x = 20 # TODO: SCALING IN X DIRECTION SHOULD BE A FUNCTION OF THE CONE ANGLE
L_yz = 2 # TODO: THIS SHOULD BE A FUNCTION OF THE PIXELS PER RADIUS
#L_x = 20 # TODO: SCALING IN X DIRECTION SHOULD BE A FUNCTION OF THE CONE ANGLE
#L_yz = 2 # TODO: THIS SHOULD BE A FUNCTION OF THE PIXELS PER RADIUS
struct = psiXseries[(psiXseries.shape[0])//2 - L_x:(psiXseries.shape[0])//2 + L_x + 1,
(psiXseries.shape[1])//2 - L_yz:(psiXseries.shape[1])//2 + L_yz + 1,
(psiXseries.shape[2])//2 - L_yz:(psiXseries.shape[2])//2 + L_yz + 1]
#struct = psiXseries[(psiXseries.shape[0])//2 - L_x:(psiXseries.shape[0])//2 + L_x + 1,
#(psiXseries.shape[1])//2 - L_yz:(psiXseries.shape[1])//2 + L_yz + 1,
#(psiXseries.shape[2])//2 - L_yz:(psiXseries.shape[2])//2 + L_yz + 1]
fXconvolvedSeries = scipy.ndimage.convolve(fXseries,struct/struct.sum())
#if useCache and not loadedCache:
#tifffile.imsave(f'{cacheFile[:-4]}_struct.tif', struct.astype('<f4'))
#tifffile.imsave(f'{cacheFile[:-4]}_fXconvolvedSeries.tif', fXconvolvedSeries.astype('<f4'))
#fXconvolvedSeries = scipy.ndimage.convolve(fXseries,struct/struct.sum())
##if useCache and not loadedCache:
##tifffile.imsave(f'{cacheFile[:-4]}_struct.tif', struct.astype('<f4'))
##tifffile.imsave(f'{cacheFile[:-4]}_fXconvolvedSeries.tif', fXconvolvedSeries.astype('<f4'))
binaryPeaks = fXconvolvedSeries > massThreshold
#binaryPeaks = fXconvolvedSeries > massThreshold
zoomLevel = sourceDetectorDistMM/((CORxPositions[0] + CORxPositions[-1])/2)
CORxDelta = numpy.abs(CORxPositions[0]-CORxPositions[1])
# Look in a volume of +/- half a radius in all directions for the highest value (+/- 1 radius keeps overlapping and causing issues, half a radius doesn't overlap particles, but still contains one clean peak)
fXconvolvedMaximumFiltered = scipy.ndimage.maximum_filter(fXconvolvedSeries,
size=(numpy.int(numpy.floor(radiusMM/CORxDelta)),
fXseriesMaximumFiltered = scipy.ndimage.maximum_filter(fXseries,
size=(3*numpy.int(numpy.floor(radiusMM/CORxDelta)),
numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)),
numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel))
)
)
allPeaks = fXconvolvedSeries == fXconvolvedMaximumFiltered
masses = allPeaks*fXconvolvedSeries
numpy.int(numpy.floor(radiusMM/pixelSizeMM*zoomLevel)))
)
allPeaks = fXseries == fXseriesMaximumFiltered
masses = allPeaks*fXseries
if verbose:
tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype('<f4'))
tifffile.imsave(cacheFile[:-4] + '_peaks.tif', allPeaks.astype('<f4'))
tifffile.imsave(cacheFile[:-4] + '_masses.tif', masses.astype('<f4'))
tifffile.imsave(cacheFile[:-4] + '_fXseries.tif', fXseries.astype('<f4'))
tifffile.imsave(cacheFile[:-4] + '_fXconvolvedMaximumFiltered.tif', fXconvolvedMaximumFiltered.astype('<f4'))
tifffile.imsave(cacheFile[:-4] + '_fXseriesMaximumFiltered.tif', fXseriesMaximumFiltered.astype('<f4'))
if scanFixedNumber:
# get the indices of all of the peaks, from highest to lowest
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment