Skip to content

fix dicom2nifti for julia v1.2 #27

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

fix dicom2nifti for julia v1.2 #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 64 additions & 69 deletions examples/dicom2nifti.jl
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

# dicom2nifti.jl
# Simple DICOM to NIfTI converter

Expand All @@ -8,7 +14,7 @@
# work with other scanners. It may give incorrect information. Test thoroughly
# and use at your own risk.

using DICOM, NIfTI, ArgParse
using DICOM, NIfTI, ArgParse, LinearAlgebra

function parse_commandline()
s = ArgParseSettings()
Expand Down Expand Up @@ -39,28 +45,6 @@ function parse_commandline()
parse_args(s)
end

function walk(fn::Function, path::String)
for fname in readdir(path)
if fname == "." || fname == ".."
continue
end
fpath = joinpath(path, fname)
if isdir(lstat(fpath))
walk(fn, fpath)
else
fn(fpath)
end
end
end

function findtag(d, tag, from, what)
tag = lookup(d, tag)
if tag == false
error("$what missing from $from")
end
tag.data
end

function main()
cmd = parse_commandline()

Expand All @@ -79,79 +63,85 @@ function main()

# Load all DICOMs into an array, indexed by series number
dicoms = Dict{Int, Any}()
walk(cmd["dicomdir"]) do fpath
if !endswith(fpath, cmd["extension"])
return
end
for (root, dirs, files) in walkdir(cmd["dicomdir"])
for file in files
fpath = joinpath(root, file)
if !endswith(fpath, cmd["extension"])
return
end

local d
f = open(fpath, "r")
d = DICOM.dcm_parse(f)
local d
println(fpath)
d = DICOM.dcm_parse(fpath)

series_number_index = findfirst((x) -> x.tag == (0x0020, 0x0011), d)
if series_number_index == 0
println("No series number found for $fpath; ignoring")
return
end
series_number = d[series_number_index].data[1]
if !(tag"Series Number" in keys(d))
println("No series number found for $fpath; ignoring")
return
end
series_number = d[tag"Series Number"]

dicom_arr = get(dicoms, series_number, nothing)
if dicom_arr == nothing
dicoms[series_number] = dicom_arr = {}
dicom_arr = get(dicoms, series_number, nothing)
if dicom_arr == nothing
dicoms[series_number] = dicom_arr = []
end
push!(dicom_arr, d)
end
push!(dicom_arr, d)
end

# Sort DICOMs
for (series_number, slices) in dicoms
d = slices[1]

protocol_name = findtag(d, (0x0018, 0x1030), series_number, "Protocol Name")[1]
orientation = findtag(d, (0x0020, 0x0037), series_number,
"Image Orientation (Patient)")
pixel_spacing = findtag(d, (0x0028, 0x0030), series_number, "Pixel Spacing")
slice_thickness = findtag(d, (0x0018, 0x0050), series_number, "Slice Thickness")[1]
time_step = lookup(d, (0x0018, 0x0080)) # This is the TR, but may not be the time step
phase_encoding_direction = lookup(d, (0x0018, 0x1312))
protocol_name = d[tag"Protocol Name"]
orientation = d[tag"Image Orientation (Patient)"]
pixel_spacing = d[tag"Pixel Spacing"]
slice_spacing = slice_thickness = d[tag"Slice Thickness"]
if tag"Spacing Between Slices" in keys(d)
slice_spacing = d[tag"Spacing Between Slices"]
end
time_step = false
if tag"Repetition Time" in keys(d)
time_step = d[tag"Repetition Time"] # This is the TR, but may not be the time step
end
phase_encoding_direction = d[tag"In-plane Phase Encoding Direction"]
pixel_rows = d[tag"Rows"]
pixel_cols = d[tag"Columns"]

# Convert orentation to RAS
# Convert orientation to RAS
orientation = transform*reshape(orientation, (3, 2))

# Determine Z as cross product of X and Y
orientation = hcat(orientation, cross(orientation[:, 1], orientation[:, 2]))

# Scale by pixel size
orientation = orientation*[pixel_spacing[1] 0 0; 0 pixel_spacing[2] 0; 0 0 slice_thickness]
orientation = orientation*[pixel_spacing[1] 0 0; 0 pixel_spacing[2] 0; 0 0 slice_spacing]

# Find Z coordinates of each slice
image_positions = transform*hcat(
[float32(findtag(d, (0x0020, 0x0032), series_number, "Image Position (Patient)"))
for d in slices]...)
image_positions = transform * hcat([d[tag"Image Position (Patient)"] for d in slices]...)
z_coords = [dot(image_positions[:, i], orientation[:, 3])
for i = 1:size(image_positions, 2)]

# Sort Z coordinates
p = sortperm(z_coords)

# Concatenate slices along Z axis according to increasing Z coordinate
raw = cat(3,
[findtag(slices[i], (0x7FE0, 0x0010), series_number, "Pixel Data")[1] for i in p]...)
raw = cat([slices[i][tag"Pixel Data"] for i in p]..., dims = 3)

orientation = float32(hcat(orientation, image_positions[:, p[1]]))
orientation = hcat(orientation, image_positions[:, p[1]])

dim_info = phase_encoding_direction.data[1] == "ROW" ? (1, 2, 3) :
phase_encoding_direction.data[1] == "COL" ? (2, 1, 3) :
dim_info = phase_encoding_direction == "ROW" ? (1, 2, 3) :
phase_encoding_direction == "COL" ? (2, 1, 3) :
(0, 0, 0)

voxel_size = tuple(pixel_spacing..., slice_thickness)
voxel_size = tuple(pixel_spacing..., slice_spacing)

# Permute according to axes specified on command line
if axes != nothing
# Determine permutation of current volume to RAS
ras = Array(Int, 3)
sign = Array(Bool, 3)
ras = Array{Int,1}(undef, 3)
sign = Array{Bool,1}(undef, 3)
rg = [1:3;]
abs_orientation = abs(orientation)
abs_orientation = abs.(orientation)
for i = 1:3
idx = findmax(abs_orientation[i, rg])[2]
ras[i] = splice!(rg, idx)
Expand All @@ -172,27 +162,32 @@ function main()
dim_info = dim_info[output_axes]

# Flip signs and dimensions
flip_sign = find(sign[specified_axes] .!= specified_sign)
flip_sign = findall(sign[specified_axes] .!= specified_sign)
orientation[:, flip_sign] = -orientation[:, flip_sign]
for dim in flip_sign
raw = flipdim(raw, dim)
orientation[dim, 4] += (sign[specified_axes[dim]] ? 1 : -1) *
raw = reverse(raw, dims=dim)
orientation[dim, 4] += (sign[specified_axes[dim]] ? 1 : -1) *
(size(raw, dim) - 1)*voxel_size[dim]
end
end

# Create a directory for each protocol
protocol_dir = joinpath(cmd["niftidir"], replace(protocol_name, "/", ""))
protocol_dir = joinpath(cmd["niftidir"], replace(protocol_name, "/" => ""))

if !isdir(cmd["niftidir"])
mkdir(cmd["niftidir"])
end
if !isdir(protocol_dir)
mkdir(protocol_dir)
end

# Write NIfTI volumes
ni = NIVolume(raw; voxel_size=voxel_size,
orientation=orientation, dim_info=dim_info,
time_step=time_step != false && !isempty(time_step.data) ? time_step.data[1] : 0f0)
# convert to Float32 as mricron and fsleyes apparently do not support UInt16
ni = NIVolume(Array{Float32, 3}(raw); voxel_size=voxel_size,
orientation=convert(Array{Float32, 2}, orientation), dim_info=dim_info,
time_step=time_step != false ? time_step : 0f0)
niwrite(joinpath(protocol_dir, "$(series_number).nii"), ni)
end
end

main()
main()
2 changes: 1 addition & 1 deletion src/NIfTI.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -363,7 +363,7 @@ function NIVolume(
slice_start, (qfac, voxel_size..., time_step, 0, 0, 0), 352,
scl_slope, scl_inter, slice_end, slice_code,
xyzt_units, cal_max, cal_min, slice_duration,
toffset, glmax, glmin, string_tuple(descrip, 80), string_tuple(aux_file, 24), (method2 || method3),
toffset, glmax, glmin, string_tuple(descrip, 80), string_tuple(aux_file, 24), method2,
method3, quatern_b, quatern_c, quatern_d,
qoffset_x, qoffset_y, qoffset_z, (orientation[1, :]...,),
(orientation[2, :]...,), (orientation[3, :]...,), string_tuple(intent_name, 16), NP1_MAGIC), extensions, raw)
Expand Down