Skip to content

Commit

Permalink
Merge pull request #41 from notZaki/apply_slope
Browse files Browse the repository at this point in the history 
Function for rescaling pixel data
  • Loading branch information
@notZaki
notZaki authored Dec 10, 2019
2 parents 20d1607 + bf75ffc commit 8d94ebf
Show file tree
Hide file tree
Showing 2 changed files with 92 additions and 53 deletions.
72 changes: 46 additions & 26 deletions src/DICOM.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module DICOM

export dcm_parse, dcm_write, lookup, lookup_vr
export dcm_parse, dcm_write, lookup, lookup_vr, rescale!
export @tag_str

"""
Expand Down Expand Up @@ -103,6 +103,23 @@ function lookup(d::Dict{Tuple{UInt16,UInt16},Any}, fieldnameString::String)
return (get(d, fieldname_dict[fieldnameString], nothing))
end

function rescale!(dcm::Dict{Tuple{UInt16,UInt16},Any}, direction = :forward)
if !haskey(dcm, tag"Rescale Intercept") || !haskey(dcm, tag"Rescale Slope")
return dcm
end
if direction == :forward
dcm[tag"Pixel Data"] =
@. dcm[tag"Pixel Data"] * dcm[tag"Rescale Slope"] + dcm[tag"Rescale Intercept"]
else
pixel_data = dcm[tag"Pixel Data"]
@. pixel_data -= dcm[tag"Rescale Intercept"]
@. pixel_data /= dcm[tag"Rescale Slope"]
dtype = determine_dtype(dcm)
dcm[tag"Pixel Data"] = @. convert(dtype, round(pixel_data))
end
return dcm
end

if ENDIAN_BOM == 0x04030201
order(x, endian) = endian == :little ? x : bswap.(x)
else
Expand Down Expand Up @@ -393,28 +410,7 @@ end

# always little-endian, "encapsulated" iff sz==0xffffffff
function pixeldata_parse(st::IO, sz, vr::String, dcm, endian)
# (0x0028,0x0103) defines Pixel Representation
isSigned = false
f = get(dcm, (0x0028, 0x0103), nothing)
if f !== nothing
# Data is signed if f==1
isSigned = f == 1
end
# (0x0028,0x0100) defines Bits Allocated
bitType = 16
f = get(dcm, (0x0028, 0x0100), nothing)
if f !== nothing
bitType = Int(f)
else
f = get(dcm, (0x0028, 0x0101), nothing)
bitType = f !== nothing ? Int(f) : vr == "OB" ? 8 : 16
end
if bitType == 8
dtype = isSigned ? Int8 : UInt8
else
dtype = isSigned ? Int16 : UInt16
end

dtype = determine_dtype(dcm)
yr = 1
zr = 1
# (0028,0010) defines number of rows
Expand Down Expand Up @@ -454,7 +450,7 @@ function pixeldata_parse(st::IO, sz, vr::String, dcm, endian)
data_dims = data_dims[data_dims.>1]
data = Array{dtype}(undef, data_dims...)
read!(st, data)
# Permute because Julia is column-major while DICOM is row-major
# Permute because Julia is column-major while DICOM is row-major
numdims = ndims(data)
if numdims == 2
perm = (2, 1)
Expand Down Expand Up @@ -486,6 +482,31 @@ function pixeldata_parse(st::IO, sz, vr::String, dcm, endian)
return order.(data, endian)
end

function determine_dtype(dcm)
# (0x0028,0x0103) defines Pixel Representation
is_signed = false
f = get(dcm, (0x0028, 0x0103), nothing)
if f !== nothing
# Data is signed if f==1
is_signed = f == 1
end
bit_type = 16
# (0x0028,0x0100) defines Bits Allocated
f = get(dcm, (0x0028, 0x0100), nothing)
if f !== nothing
bit_type = Int(f)
else
f = get(dcm, (0x0028, 0x0101), nothing)
bit_type = f !== nothing ? Int(f) : vr == "OB" ? 8 : 16
end
if bit_type == 8
dtype = is_signed ? Int8 : UInt8
else
dtype = is_signed ? Int16 : UInt16
end
return dtype
end

function undefined_length(st, vr)
data = IOBuffer()
w1 = w2 = 0
Expand Down Expand Up @@ -526,7 +547,6 @@ function dcm_write(
write(st, "DICM")
end
(is_explicit, endian) = determine_explicitness_and_endianness(dcm)

for gelt in sort(collect(keys(dcm)))
write_element(st, gelt, dcm[gelt], is_explicit, aux_vr)
end
Expand Down Expand Up @@ -656,7 +676,7 @@ function pixeldata_write(st, d, evr)
vr = nt === UInt8 || nt === Int8 ? "OB" :
nt === UInt16 || nt === Int16 ? "OW" :
nt === Float32 ? "OF" : error("dicom: unsupported pixel format")
# Permute because Julia is column-major while DICOM is row-major
# Permute because Julia is column-major while DICOM is row-major
# !warn! This part assumes that Planar Configuration (tag: 0x0028, 0x0006) is not 0
numdims = ndims(d)
perm = numdims == 2 ? (2, 1) : numdims == 3 ? (2, 1, 3) : (2, 1, 3, 4)
Expand Down
73 changes: 46 additions & 27 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,24 @@ if !isdir(data_folder)
end

const dicom_samples = Dict(
"CT_Explicit_Little.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/CT_Explicit_Little.dcm",
"CT_Implicit_Little_Headless_Retired.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/CT_Implicit_Little_Headless_Retired.dcm",
"MG_Explicit_Little.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MG_Explicit_Little.dcm",
"MR_Explicit_Little_MultiFrame.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MR_Explicit_Little_MultiFrame.dcm",
"MR_Implicit_Little.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MR_Implicit_Little.dcm",
"MR_UnspecifiedLength.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MR_UnspecifiedLength.dcm",
"OT_Implicit_Little_Headless.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/OT_Implicit_Little_Headless.dcm",
"US_Explicit_Big_RGB.dcm" => "https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/US_Explicit_Big_RGB.dcm",
"DX_Implicit_Little_Interleaved.dcm" => "https://github.com/OHIF/viewer-testdata/raw/master/dcm/zoo-exotic/5.dcm"
"CT_Explicit_Little.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/CT_Explicit_Little.dcm",
"CT_Implicit_Little_Headless_Retired.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/CT_Implicit_Little_Headless_Retired.dcm",
"MG_Explicit_Little.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MG_Explicit_Little.dcm",
"MR_Explicit_Little_MultiFrame.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MR_Explicit_Little_MultiFrame.dcm",
"MR_Implicit_Little.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MR_Implicit_Little.dcm",
"MR_UnspecifiedLength.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/MR_UnspecifiedLength.dcm",
"OT_Implicit_Little_Headless.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/OT_Implicit_Little_Headless.dcm",
"US_Explicit_Big_RGB.dcm" =>
"https://github.com/notZaki/DICOMSamples/raw/master/DICOMSamples/US_Explicit_Big_RGB.dcm",
"DX_Implicit_Little_Interleaved.dcm" =>
"https://github.com/OHIF/viewer-testdata/raw/master/dcm/zoo-exotic/5.dcm",
)

function download_dicom(filename; folder = data_folder)
Expand Down Expand Up @@ -73,18 +82,20 @@ end
outMG2 = joinpath(data_folder, "outMG2.dcm")

# Write DICOM files
outIO = open(outMR1, "w+"); dcm_write(outIO, dcmMR); close(outIO)
outIO = open(outCT1, "w+"); dcm_write(outIO, dcmCT); close(outIO)
outIO = open(outMG1, "w+"); dcm_write(outIO, dcmMG, aux_vr = vrMG); close(outIO)
dcm_write(outMG1b, dcmMG, aux_vr = vrMG)
dcm_write(outMR1, dcmMR)
dcm_write(outCT1, dcmCT)
dcm_write(outMG1, dcmMG; aux_vr = vrMG)
open(outMG1b, "w") do io
dcm_write(io, dcmMG; aux_vr = vrMG)
end
# Reading DICOM files which were written from previous step
dcmMR1 = dcm_parse(outMR1)
dcmCT1 = dcm_parse(outCT1)
(dcmMG1, vrMG1) = dcm_parse(outMG1, return_vr = true)
# Write DICOM files which were re-read from previous step
outIO = open(outMR2, "w+"); dcm_write(outIO, dcmMR1); close(outIO)
outIO = open(outCT2, "w+"); dcm_write(outIO, dcmCT1); close(outIO)
outIO = open(outMG2, "w+"); dcm_write(outIO, dcmMG1, aux_vr = vrMG1); close(outIO)
dcm_write(outMR2, dcmMR1)
dcm_write(outCT2, dcmCT1)
dcm_write(outMG2, dcmMG1; aux_vr = vrMG1)

# Test consistency of written files after the write-read-write cycle
@test read(outMR1) == read(outMR2)
Expand Down Expand Up @@ -125,18 +136,26 @@ end
(0x0020, 0x0070) => "LO",
(0x0028, 0x0005) => "US",
(0x0028, 0x0040) => "CS",
(0x0028, 0x0200) => "US")
dcmCTa = dcm_parse(fileCT, preamble = false, aux_vr = dVR_CTa);
(0x0028, 0x0200) => "US",
)
dcmCTa = dcm_parse(fileCT, preamble = false, aux_vr = dVR_CTa)
# 2b. Read with a master VR which skips elements
# Here we skip any element where lookup_vr() fails
# And we also force (0x0018,0x1170) to be read as float instead of integer
dVR_CTb = Dict( (0x0000, 0x0000) => "", (0x0018, 0x1170) => "DS")
dcmCTb = dcm_parse(fileCT, preamble = false, aux_vr = dVR_CTb);
dVR_CTb = Dict((0x0000, 0x0000) => "", (0x0018, 0x1170) => "DS")
dcmCTb = dcm_parse(fileCT, preamble = false, aux_vr = dVR_CTb)
@test dcmCTa[(0x0008, 0x0060)] == "CT"
@test dcmCTb[(0x0008, 0x0060)] == "CT"
@test haskey(dcmCTa, (0x0028, 0x0040)) # dcmCTa should contain retired element
@test !haskey(dcmCTb, (0x0028, 0x0040)) # dcmCTb skips retired elements

rescale!(dcmCTa)
@test minimum(dcmCTa[(0x7fe0, 0x0010)]) == -949
@test maximum(dcmCTa[(0x7fe0, 0x0010)]) == 1132
rescale!(dcmCTa, :backward)
@test minimum(dcmCTa[(0x7fe0, 0x0010)]) == minimum(dcmCTb[(0x7fe0, 0x0010)])
@test maximum(dcmCTa[(0x7fe0, 0x0010)]) == maximum(dcmCTb[(0x7fe0, 0x0010)])

# 3. DICOM file containing multiple frames
fileMR_multiframe = download_dicom("MR_Explicit_Little_MultiFrame.dcm")
dcmMR_multiframe = dcm_parse(fileMR_multiframe)
Expand All @@ -162,14 +181,14 @@ end
end

@testset "Test tag macro" begin
@test tag"Modality" === (0x0008, 0x0060) ===
DICOM.fieldname_dict["Modality"]
@test tag"Shutter Overlay Group" === (0x0018, 0x1623) ===
DICOM.fieldname_dict["Shutter Overlay Group"]
@test tag"Modality" === (0x0008, 0x0060) === DICOM.fieldname_dict["Modality"]
@test tag"Shutter Overlay Group" ===
(0x0018, 0x1623) ===
DICOM.fieldname_dict["Shutter Overlay Group"]
@test tag"Histogram Last Bin Value" === (0x0060, 0x3006)
DICOM.fieldname_dict["Histogram Last Bin Value"]
DICOM.fieldname_dict["Histogram Last Bin Value"]

# test that compile time error is thrown if tag does not exist
@test macroexpand(Main,:(tag"Modality")) === (0x0008, 0x0060)
@test_throws LoadError macroexpand(Main,:(tag"nonsense"))
@test macroexpand(Main, :(tag"Modality")) === (0x0008, 0x0060)
@test_throws LoadError macroexpand(Main, :(tag"nonsense"))
end

0 comments on commit 8d94ebf

Please sign in to comment.