Skip to content
Snippets Groups Projects
Commit 34259d8a authored by David Johansen's avatar David Johansen
Browse files

fix dtype bugs and forgetting to use array method on CIL objects

parent 6f286e34
No related branches found
No related tags found
No related merge requests found
{
// Use IntelliSense to learn about possible attributes.
// Hover to view descriptions of existing attributes.
// For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
"version": "0.2.0",
"configurations": [
{
"name": "Python Debugger: Remote Attach",
"type": "debugpy",
"request": "attach",
"connect": {
"host": "10.66.20.9",
"port": 44123
},
"pathMappings": [
{
"localRoot": "${workspaceFolder}",
"remoteRoot": "."
}
]
}
]
}
\ No newline at end of file
import zarr
import numpy as np
import os
import matplotlib.pyplot as plt
from cil.io import ZEISSDataReader, NikonDataReader
from cil.framework import AcquisitionData
......@@ -25,25 +27,30 @@ def write_projection_chunks(config_path, zarr_path, chunk_size, proj_shape):
out_path: Directory to write zarr file.
"""
nchunks = (proj_shape[0]-1)//chunk_size + 1
compressors = zarr.codecs.Blosc(cname='lz4')
zarr_array = zarr.open(
# compressor = zarr.codecs.Blosc(cname='lz4')
compressor = zarr.codecs.Blosc(cname='zstd')
zarr_array = zarr.open_array(
store=zarr_path,
mode='w',
dtype=np.uint16,
dtype=np.float32,
shape=proj_shape,
chunks=(chunk_size, *proj_shape[1:]),
compressors=compressors)
compressor=compressor
)
for i in range(nchunks):
start = i*chunk_size
end = min((i+1)*chunk_size, proj_shape[0])
reader = create_reader(config_path, roi={'angle': (start, end, 1)})
chunk = reader.read()
chunk = reader.read().as_array()
print(f'Write {i}: chunk.min(),chunk.mean()={chunk.min()},{chunk.mean()}')
expected_shape = (end-start, *proj_shape[1:])
if chunk.shape != expected_shape:
raise ValueError(f'Chunk shape mismatch: expected {expected_shape}, got {chunk.shape}')
zarr_array[start:end,:,:] = chunk
print(f'Finished writing chunk {i}.')
def read_projection_chunk(zarr_path, chunk_size, chunk_index):
zarr_array = zarr.open_array(store=zarr_path, mode='r')
......@@ -51,6 +58,7 @@ def read_projection_chunk(zarr_path, chunk_size, chunk_index):
start = chunk_index*chunk_size
end = min((chunk_index+1)*chunk_size, proj_shape[0])
chunk = zarr_array[start:end,:,:]
print(f'Finished reading projection {chunk_index}.')
return chunk
def recon_projection_chunk(ag, ig, zarr_path, chunk_size, chunk_index):
......@@ -62,20 +70,25 @@ def recon_projection_chunk(ag, ig, zarr_path, chunk_size, chunk_index):
end = min((chunk_index+1)*chunk_size, ag.shape[0])
ag_chunk = ag.copy().set_angles(ag.angles[start:end])
chunk = read_projection_chunk(zarr_path, chunk_size, chunk_index)
chunk = chunk.astype(np.float32)
print(f'Read {chunk_index}: chunk.min(),chunk.mean()={chunk.min()},{chunk.mean()}')
data_chunk = AcquisitionData(array=chunk, deep_copy=False, geometry=ag_chunk)
data_chunk = TransmissionAbsorptionConverter()(data_chunk)
data_chunk.reorder(order='tigre')
recon_chunk = FDK(data_chunk, image_geometry=ig).run(verbose=1)
print(f'Finished reconstruction chunk {chunk_index}.')
return recon_chunk
def chunk_volumes():
pass
def write_numpy(arr, file_name):
np.save(f"{file_name}.npy", arr)
os.sync()
def main():
path = '/dtu/3d-imaging-center/projects/2024_QIM_platform/analysis/data/valnut/valnut_2014-03-21_643_28/tomo-A/valnut_tomo-A.txrm'
chunk_size = 600
zarr_path = 'ct_chunk.zarr'
path = '/dtu/3d-imaging-center/projects/2021_DANFIX_Casper/raw_data_3DIM/Casper_top_3_2 [2021-03-17 16.54.39]/Casper_top_3_2_recon.xtekct'
chunk_size = 300
zarr_path = 'out/ct_chunk.zarr'
ag = create_reader(path).get_geometry()
ig = ag.get_ImageGeometry()
......@@ -88,5 +101,16 @@ def main():
nchunks = (proj_shape[0]-1)//chunk_size + 1
recon_accumulated = np.zeros(ig.shape, dtype=np.float32)
for chunk_index in range(nchunks):
recon_accumulated += recon_projection_chunk(ag, ig, zarr_path, chunk_size=chunk_size, chunk_index=chunk_index)
recon_chunk = recon_projection_chunk(ag, ig, zarr_path, chunk_size=chunk_size, chunk_index=chunk_index)
recon_accumulated += recon_chunk.as_array()
plt.figure(figsize=(8,8))
plt.imshow(recon_chunk.as_array()[ig.shape[0]//2,:,:])
plt.colorbar()
plt.savefig(f'out/center_recon_{chunk_index}.png')
# zarr_store = zarr.open('ct_recon.zarr', mode='w', shape=ig.shape, dtype=np.float32)
write_numpy(recon_accumulated.as_array(), 'recon')
if __name__ == '__main__':
main()
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment