diff --git a/julia/FreeSurfer/Manifest.toml b/julia/FreeSurfer/Manifest.toml new file mode 100644 index 00000000000..34622d4ca8f --- /dev/null +++ b/julia/FreeSurfer/Manifest.toml @@ -0,0 +1,341 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.7.1" +manifest_format = "2.0" + +[[deps.AbstractFFTs]] +deps = ["ChainRulesCore", "LinearAlgebra"] +git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4" +uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c" +version = "1.1.0" + +[[deps.Adapt]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "3.3.3" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "7dd38532a1115a215de51775f9891f0f3e1bac6a" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.12.1" + +[[deps.ChangesOfVariables]] +deps = ["ChainRulesCore", "LinearAlgebra", "Test"] +git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1" +uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" +version = "0.1.2" + +[[deps.ColorTypes]] +deps = ["FixedPointNumbers", "Random"] +git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597" +uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +version = "0.11.0" + +[[deps.ColorVectorSpace]] +deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"] +git-tree-sha1 = "3f1f500312161f1ae067abe07d13b40f78f32e07" +uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4" +version = "0.9.8" + +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40" +uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" +version = "0.12.8" + +[[deps.Compat]] +deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] +git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "3.41.0" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" + +[[deps.Crayons]] +git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15" +uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f" +version = "4.1.1" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DelimitedFiles]] +deps = ["Mmap"] +uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.8.6" + +[[deps.Downloads]] +deps = ["ArgTools", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" + +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.4" + +[[deps.Graphics]] +deps = ["Colors", "LinearAlgebra", "NaNMath"] +git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2" +uuid = "a2bd30eb-e257-5431-a919-1863eab51364" +version = "1.1.1" + +[[deps.ImageBase]] +deps = ["ImageCore", "Reexport"] +git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835" +uuid = "c817782e-172a-44cc-b673-b171935fbb9e" +version = "0.1.5" + +[[deps.ImageCore]] +deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"] +git-tree-sha1 = "9a5c62f231e5bba35695a20988fc7cd6de7eeb5a" +uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534" +version = "0.9.3" + +[[deps.ImageInTerminal]] +deps = ["Crayons", "ImageBase", "ImageCore", "Requires"] +git-tree-sha1 = "2bd30630d417699c925fa7f16d51282d7343e3eb" +uuid = "d8c32880-2388-543b-8c61-d9f865259254" +version = "0.4.8" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.InverseFunctions]] +deps = ["Test"] +git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.2" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.1.1" + +[[deps.JLLWrappers]] +deps = ["Preferences"] +git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.4.1" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" + +[[deps.LibGit2]] +deps = ["Base64", "NetworkOptions", "Printf", "SHA"] +uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.LogExpFunctions]] +deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.6" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.MappedArrays]] +git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142" +uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900" +version = "0.4.1" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.MosaicViews]] +deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"] +git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d" +uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389" +version = "0.3.3" + +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" + +[[deps.NaNMath]] +git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "0.3.7" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" + +[[deps.OffsetArrays]] +deps = ["Adapt"] +git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.10.8" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] +git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.5+0" + +[[deps.PaddedViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "03a7a85b76381a3d04c7a1656039197e70eda03d" +uuid = "5432bcbf-9aad-5242-b902-cca2824c8663" +version = "0.5.11" + +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] +uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.2.3" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA", "Serialization"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SparseArrays]] +deps = ["LinearAlgebra", "Random"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[[deps.SpecialFunctions]] +deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "8d0c8e3d0ff211d9ff4a0c2307d876c99d10bdf1" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.1.2" + +[[deps.StackViews]] +deps = ["OffsetArrays"] +git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c" +uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15" +version = "0.1.1" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" + +[[deps.TensorCore]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6" +uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50" +version = "0.1.1" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl", "OpenBLAS_jll"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" diff --git a/julia/FreeSurfer/Project.toml b/julia/FreeSurfer/Project.toml new file mode 100644 index 00000000000..daea9b0900f --- /dev/null +++ b/julia/FreeSurfer/Project.toml @@ -0,0 +1,12 @@ +name = "FreeSurfer" +uuid = "cbb28e16-2113-4377-92c1-bba511c54513" +authors = ["Anastasia Yendiki "] +version = "0.1.0" + +[deps] +ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +ImageInTerminal = "d8c32880-2388-543b-8c61-d9f865259254" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/julia/FreeSurfer/README.md b/julia/FreeSurfer/README.md new file mode 100644 index 00000000000..d93d84a85ce --- /dev/null +++ b/julia/FreeSurfer/README.md @@ -0,0 +1,45 @@ +### Import this package + +Before starting julia, it is recommended to define the environment variable FREESURFER_HOME. This will be used, e.g., to load the FreeSurfer color look-up table automatically. + +```julia +julia> import FreeSurfer as fs +FREESURFER_HOME: /usr/local/freesurfer/dev +``` + +### Read .mgh, .mgz, .nii, .nii.gz volumes + +```julia +julia> aa = fs.mri_read("/usr/local/freesurfer/dev/subjects/fsaverage/mri/aparc+aseg.mgz"); + +julia> fa = fs.mri_read("/usr/local/freesurfer/dev/trctrain/hcp/MGH35_HCP_FA_template.nii.gz"); +``` + +### Show volume and header summary info + +```julia +julia> fs.show(aa) + +julia> fs.show(fa) +``` + +### Write .mgh, .mgz, .nii, .nii.gz volumes + +```julia +julia> fs.mri_write(aa, "/tmp/aparc+aseg.nii.gz") + +julia> fs.mri_write(fa, "/tmp/MGH35_HCP_FA_template.mgz") +``` + +### Read a .trk tractography streamline file + +```julia +julia> tr = fs.trk_read("/usr/local/freesurfer/dev/trctrain/hcp/mgh_1001/syn/acomm.bbr.prep.trk"); +``` + +### Write a .trk tractography streamline file + +```julia +julia> fs.trk_write(tr, "/tmp/acomm.trk") +``` + diff --git a/julia/FreeSurfer/src/FreeSurfer.jl b/julia/FreeSurfer/src/FreeSurfer.jl new file mode 100644 index 00000000000..295abf6a9af --- /dev/null +++ b/julia/FreeSurfer/src/FreeSurfer.jl @@ -0,0 +1,28 @@ +#= + Original Author: Anastasia Yendiki + + Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH" + + Terms and conditions for use, reproduction, distribution and contribution + are found in the 'FreeSurfer Software License Agreement' contained + in the file 'LICENSE' found in the FreeSurfer distribution, and here: + + https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense + + Reporting: freesurfer@nmr.mgh.harvard.edu +=# + +module FreeSurfer + +include("mri.jl") +include("trk.jl") +include("view.jl") +include("dti.jl") + +function __init__() + + println("FREESURFER_HOME: " * (haskey(ENV, "FREESURFER_HOME") ? + ENV["FREESURFER_HOME"] : "not defined")) +end + +end # module diff --git a/julia/FreeSurfer/src/dti.jl b/julia/FreeSurfer/src/dti.jl new file mode 100644 index 00000000000..97e25af1f2e --- /dev/null +++ b/julia/FreeSurfer/src/dti.jl @@ -0,0 +1,147 @@ +#= + Original Author: Anastasia Yendiki + + Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH" + + Terms and conditions for use, reproduction, distribution and contribution + are found in the 'FreeSurfer Software License Agreement' contained + in the file 'LICENSE' found in the FreeSurfer distribution, and here: + + https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense + + Reporting: freesurfer@nmr.mgh.harvard.edu +=# + +using LinearAlgebra, Statistics + +export DTI, dti_fit, dti_write + + +"Container for outputs of a DTI fit" +struct DTI + eigval1::MRI + eigval2::MRI + eigval3::MRI + eigvec1::MRI + eigvec2::MRI + eigvec3::MRI + rd::MRI + md::MRI + fa::MRI +end + + +""" + dti_fit_ls(dwi::MRI, mask:MRI) + +Fit tensors to DWIs and return a `DTI` structure. +""" +function dti_fit(dwi::MRI, mask::MRI) + dti_fit_ls(dwi::MRI, mask::MRI) +end + + +""" + dti_fit_ls(dwi::MRI, mask:MRI) + +Perform least-squares fitting of tensors from DWIs and return a `DTI` structure. +""" +function dti_fit_ls(dwi::MRI, mask::MRI) + + if isempty(dwi.bval) + error("Missing b-value table from input volume") + end + + if isempty(dwi.bvec) + error("Missing gradient table from input volume") + end + + idwi = dwi.bval[:,1] .> 0 + b = dwi.bval[idwi] + g = dwi.bvec[idwi, :] + + bB = b .* hcat(g[:,1].^2, 2*g[:,1].*g[:,2], 2*g[:,1].*g[:,3], + g[:,2].^2, 2*g[:,2].*g[:,3], g[:,3].^2) + + pbB = pinv(bB) + + Eval1 = MRI(mask, 1) + Eval2 = MRI(mask, 1) + Eval3 = MRI(mask, 1) + Evec1 = MRI(mask, 3) + Evec2 = MRI(mask, 3) + Evec3 = MRI(mask, 3) + + for iz in 1:size(dwi.vol, 3) + for iy in 1:size(dwi.vol, 2) + for ix in 1:size(dwi.vol, 1) + if mask.vol[ix, iy, iz] > 0 + s = dwi.vol[ix, iy, iz, idwi] + s0 = mean(dwi.vol[ix, iy, iz, .!idwi]) + + D = -pbB * log.(abs.(s) / abs(s0)) + + E = eigen([D[1] D[2] D[3]; + D[2] D[4] D[5]; + D[3] D[5] D[6]]) + + Eval1.vol[ix, iy, iz] = E.values[3] + Eval2.vol[ix, iy, iz] = E.values[2] + Eval3.vol[ix, iy, iz] = E.values[1] + Evec1.vol[ix, iy, iz, :] = E.vectors[:, 3] + Evec2.vol[ix, iy, iz, :] = E.vectors[:, 2] + Evec3.vol[ix, iy, iz, :] = E.vectors[:, 1] + end + end + end + end + + return DTI(Eval1, Eval2, Eval3, Evec1, Evec2, Evec3, + dti_maps(Eval1, Eval2, Eval3)...) +end + +""" + dti_maps(eigval1::MRI, eigval2::MRI, eigval3::MRI) + +Compute radial diffusivity (RD), mean diffusivity (MD), and fractional +anisotropy (FA) maps from the 3 eigenvalues the diffusion tensors. + +Return RD, MD, and FA maps as MRI structures. +""" +function dti_maps(eigval1::MRI, eigval2::MRI, eigval3::MRI) + + rd = MRI(eigval1) + md = MRI(eigval1) + fa = MRI(eigval1) + + imask = (eigval1.vol .!= 0) + + rd.vol[imask] = eigval1.vol[imask] + eigval2.vol[imask] + md.vol[imask] = (rd.vol[imask] + eigval3.vol[imask]) / 3 + rd.vol[imask] /= 2 + + fa.vol[imask] = sqrt.(((eigval1.vol[imask] - md.vol[imask]).^2 + + (eigval2.vol[imask] - md.vol[imask]).^2 + + (eigval3.vol[imask] - md.vol[imask]).^2) ./ + (eigval1.vol[imask].^2 + + eigval2.vol[imask].^2 + + eigval3.vol[imask].^2) * 3/2) + + return rd, md, fa +end + + +""" + dti_write(dti::DTI, basename::String) + +Write the volumes from a `DTI` structure that was created by `dti_fit()` +to files whose names start with the specified base name. +""" +function dti_write(dti::DTI, basename::String) + + for var in fieldnames(DTI) + mri_write(getfield(dti, var), basename * "_" * string(var) * ".nii.gz") + end +end + + diff --git a/julia/FreeSurfer/src/mri.jl b/julia/FreeSurfer/src/mri.jl new file mode 100644 index 00000000000..d821b7b21d3 --- /dev/null +++ b/julia/FreeSurfer/src/mri.jl @@ -0,0 +1,1772 @@ +#= + Original Author: Anastasia Yendiki + + Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH" + + Terms and conditions for use, reproduction, distribution and contribution + are found in the 'FreeSurfer Software License Agreement' contained + in the file 'LICENSE' found in the FreeSurfer distribution, and here: + + https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense + + Reporting: freesurfer@nmr.mgh.harvard.edu +=# + +#= + File i/o based on MATLAB code by Doug Greve, Bruce Fischl: + MRIfspec.m + MRIread.m + MRIwrite.m + load_mgh.m + load_nifti.m + load_nifti_header.m + save_nifti.m + fsgettmppath.m + vox2ras_0to1.m + vox2ras_tkreg.m + vox2rasToQform.m +=# + +using LinearAlgebra, Printf, DelimitedFiles + +export MRI, NIfTIheader, get_tmp_path, mri_filename, mri_read, mri_write, + mri_read_bfiles, mri_read_bfiles! + + +"Container for header and image data of a volume stored in NIfTI format" +mutable struct NIfTIheader + # NIfTI standard header fields + sizeof_hdr::Int32 + data_type::Vector{UInt8} + db_name::Vector{UInt8} + extents::Int32 + session_error::Int16 + regular::UInt8 + dim_info::UInt8 + dim::Vector{UInt16} + intent_p1::Float32 + intent_p2::Float32 + intent_p3::Float32 + intent_code::Int16 + datatype::Int16 + bitpix::Int16 + slice_start::Int16 + pixdim::Vector{Float32} + vox_offset::Float32 + scl_slope::Float32 + scl_inter::Float32 + slice_end::Int16 + slice_code::Int8 + xyzt_units::Int8 + cal_max::Float32 + cal_min::Float32 + slice_duration::Float32 + toffset::Float32 + glmax::Int32 + glmin::Int32 + descrip::Vector{UInt8} + aux_file::Vector{UInt8} + qform_code::Int16 + sform_code::Int16 + quatern_b::Float32 + quatern_c::Float32 + quatern_d::Float32 + quatern_x::Float32 + quatern_y::Float32 + quatern_z::Float32 + srow_x::Vector{Float32} + srow_y::Vector{Float32} + srow_z::Vector{Float32} + intent_name::Vector{UInt8} + magic::Vector{UInt8} + + # Additional fields + do_bswap::UInt8 + sform::Matrix{Float32} + qform::Matrix{Float32} + vox2ras::Matrix{Float32} + vol::Array +end + + +""" + NIfTIheader() + +Return an empty `NIfTIheader` structure +""" +NIfTIheader() = NIfTIheader( + Int32(0), + Vector{UInt8}(undef, 0), + Vector{UInt8}(undef, 0), + Int32(0), + Int16(0), + UInt8(0), + UInt8(0), + Vector{UInt16}(undef, 0), + Float32(0), + Float32(0), + Float32(0), + Int16(0), + Int16(0), + Int16(0), + Int16(0), + Vector{Float32}(undef, 0), + Float32(0), + Float32(0), + Float32(0), + Int16(0), + Int8(0), + Int8(0), + Float32(0), + Float32(0), + Float32(0), + Float32(0), + Int32(0), + Int32(0), + Vector{UInt8}(undef, 0), + Vector{UInt8}(undef, 0), + Int16(0), + Int16(0), + Float32(0), + Float32(0), + Float32(0), + Float32(0), + Float32(0), + Float32(0), + Vector{Float32}(undef, 0), + Vector{Float32}(undef, 0), + Vector{Float32}(undef, 0), + Vector{UInt8}(undef, 0), + Vector{UInt8}(undef, 0), + + UInt8(0), + Matrix{Float32}(undef, 0, 0), + Matrix{Float32}(undef, 0, 0), + Matrix{Float32}(undef, 0, 0), + [] +) + + +"Container for header and image data of an MRI volume or volume series" +mutable struct MRI + vol::Array + ispermuted::Bool + niftihdr::NIfTIheader + + fspec::String + pwd::String + + flip_angle::Float32 + tr::Float32 + te::Float32 + ti::Float32 + + vox2ras0::Matrix{Float32} + volsize::Vector{Int32} + height::Int32 + width::Int32 + depth::Int32 + nframes::Int32 + + vox2ras::Matrix{Float32} + nvoxels::Int32 + xsize::Float32 + ysize::Float32 + zsize::Float32 + + x_r::Float32 + x_a::Float32 + x_s::Float32 + + y_r::Float32 + y_a::Float32 + y_s::Float32 + + z_r::Float32 + z_a::Float32 + z_s::Float32 + + c_r::Float32 + c_a::Float32 + c_s::Float32 + + vox2ras1::Matrix{Float32} + Mdc::Matrix{Float32} + volres::Vector{Float32} + tkrvox2ras::Matrix{Float32} + + bval::Vector{Number} + bvec::Matrix{Number} +end + + +""" + MRI() + +Return an empty `MRI` structure +""" +MRI() = MRI( + [], + false, + NIfTIheader(), + + "", + "", + + Float32(0), + Float32(0), + Float32(0), + Float32(0), + + Matrix{Float32}(undef, 0, 0), + Vector{Int32}(undef, 0), + Int32(0), + Int32(0), + Int32(0), + Int32(0), + + Matrix{Float32}(undef, 0, 0), + Int32(0), + Float32(0), + Float32(0), + Float32(0), + + Float32(0), + Float32(0), + Float32(0), + + Float32(0), + Float32(0), + Float32(0), + + Float32(0), + Float32(0), + Float32(0), + + Float32(0), + Float32(0), + Float32(0), + + Matrix{Float32}(undef, 0, 0), + Matrix{Float32}(undef, 0, 0), + Vector{Float32}(undef, 0), + Matrix{Float32}(undef, 0, 0), + + Vector{Number}(undef, 0), + Matrix{Number}(undef, 0, 0) +) + + +""" + MRI(ref::MRI, nframes::Integer=ref.nframes) + +Return an `MRI` structure whose header fields are populated based on a +reference `MRI` structure `ref`, and whose image array are populated with +zeros. + +Optionally, the new `MRI` structure can be created with a different number of +frames (`nframes`) than the reference MRI structure. +""" +function MRI(ref::MRI, nframes::Integer=ref.nframes) + + mri = MRI() + + for var in fieldnames(MRI) + any(var .== (:vol, :fspec, :bval, :bvec)) && continue + setfield!(mri, var, getfield(ref, var)) + end + + mri.nframes = nframes + mri.vol = zeros(Float32, ref.volsize..., nframes) + + return mri +end + + +""" + get_tmp_path(tmpdir::String="") + +Return path to a directory where temporary files can be stored. + +Search for candidate directories in the following order: + 1. `\$``TMPDIR`: Check if environment variable is defined and directory exists + 2. `\$``TEMPDIR`: Check if environment variable is defined and directory exists + 3. `/scratch`: Check if directory exists + 4. `/tmp`: Check if directory exists + 5. `tmpdir`: Check if `tmpdir` argument was passed and directory exists + +If none of the above exist, use current directory (`./`) and print warning. +""" +function get_tmp_path(tmpdir::String="") + + if haskey(ENV, "TMPDIR") + tmppath = ENV["TMPDIR"] + if isdir(tmppath) + return tmppath + end + end + + if haskey(ENV, "TEMPDIR") + tmppath = ENV["TEMPDIR"] + if isdir(tmppath) + return tmppath + end + end + + tmppath = "/scratch" + if isdir(tmppath) + return tmppath + end + + tmppath = "/tmp" + if isdir(tmppath) + return tmppath + end + + tmppath = tmpdir + if isdir(tmppath) + return tmppath + end + + tmppath = "./" + println("WARNING: get_tmp_path could not find a temporary folder, " * + "using current folder") + return tmppath +end + + +""" + vox2ras_0to1(M0::Matrix) + +Convert a 0-based vox2ras matrix `M0` to a 1-based vox2ras matrix such that: + +Pxyz = M_0 * [c r s 1]' = M_1 * [c+1 r+1 s+1 1]' +""" +function vox2ras_0to1(M0::Matrix) + + if size(M0) != (4,4) + error("Input must be a 4x4 matrix") + end + + Q = zeros(4, 4) + Q[1:3, 4] = ones(3, 1) + + M1 = inv(inv(M0)+Q) + + return M1 +end + + +""" + vox2ras_tkreg(voldim::Vector, voxres::Vector) + +Return a 0-based vox2ras transform of a volume that is compatible with the +registration matrix produced by tkregister. May work with MINC xfm. + +# Arguments +- voldim = [ncols; nrows; nslices ...] +- volres = [colres; rowres; sliceres ...] +""" +function vox2ras_tkreg(voldim::Vector, voxres::Vector) + + if length(voldim) < 3 | length(voxres) < 3 + error("Input vectors must have at least 3 elements") + end + + T = zeros(4,4) + T[4,4] = 1 + + T[1,1] = -voxres[1] + T[1,4] = voxres[1] * voldim[1]/2 + + T[2,3] = voxres[3] + T[2,4] = -voxres[3] * voldim[3]/2 + + T[3,2] = -voxres[2] + T[3,4] = voxres[2] * voldim[2]/2 + + return T +end + + +""" + vox2ras_to_qform(vox2ras::Matrix) + +Convert a vox2ras matrix to NIfTI qform parameters. The vox2ras should be 6 DOF. + +Return the following NIfTI header fields: +- hdr.quatern_b +- hdr.quatern_c +- hdr.quatern_d +- hdr.qoffset_x +- hdr.qoffset_y +- hdr.qoffset_z +- hdr.pixdim[1] + +From DG's vox2rasToQform.m: + This code mostly just follows CH's mriToNiftiQform() in mriio.c +""" +function vox2ras_to_qform(vox2ras::Matrix) + + if size(vox2ras) != (4, 4) + error("vox2ras size=" * string(size(vox2ras)) * ", must be (4, 4)") + end + + x = vox2ras[1,4] + y = vox2ras[2,4] + z = vox2ras[3,4] + + d = sqrt.(sum(vox2ras[:, 1:3].^2; dims=1)) + Mdc = vox2ras[1:3, 1:3] ./ repeat(d, 3) + if det(Mdc) == 0 + error("vox2ras determinant is 0") + end + + r11 = Mdc[1,1] + r21 = Mdc[2,1] + r31 = Mdc[3,1] + r12 = Mdc[1,2] + r22 = Mdc[2,2] + r32 = Mdc[3,2] + r13 = Mdc[1,3] + r23 = Mdc[2,3] + r33 = Mdc[3,3] + + if det(Mdc) > 0 + qfac = 1.0 + else + r13 = -r13 + r23 = -r23 + r33 = -r33 + qfac = -1.0 + end + + # From DG's vox2rasToQform.m: "following mat44_to_quatern()" + a = r11 + r22 + r33 + 1.0 + if a > 0.5 + a = 0.5 * sqrt(a) + b = 0.25 * (r32-r23) / a + c = 0.25 * (r13-r31) / a + d = 0.25 * (r21-r12) / a + else + xd = 1.0 + r11 - (r22+r33) + yd = 1.0 + r22 - (r11+r33) + zd = 1.0 + r33 - (r11+r22) + if xd > 1 + b = 0.5 * sqrt(xd) + c = 0.25 * (r12+r21) / b + d = 0.25 * (r13+r31) / b + a = 0.25 * (r32-r23) / b + elseif yd > 1 + c = 0.5 * sqrt(yd) + b = 0.25 * (r12+r21) / c + d = 0.25 * (r23+r32) / c + a = 0.25 * (r13-r31) / c + else + d = 0.5 * sqrt(zd) + b = 0.25 * (r13+r31) / d + c = 0.25 * (r23+r32) / d + a = 0.25 * (r21-r12) / d + end + if a < 0 + a = -a + b = -b + c = -c + d = -d + end + end + + return [b, c, d, x, y, z, qfac] +end + + +""" + mri_filename(fstring::String, checkdisk::Bool=true) + +Return a valid file name, file stem, and file extension, given a string +`fstring` that can be either a file name or a file stem. + +Valid extensions are: mgh, mgz, nii, nii.gz. Thus a file name is expected to +have the form stem.{mgh, mgz, nii, nii.gz}. + +If `fstring` is a file name, then the stem and extension are determined from +`fstring`. + +If `fstring` is a file stem and `checkdisk` is true (default), then the +file name and extension are determined by searching for a file on disk +named `fstring`.{mgh, mgz, nii, nii.gz}, where the possible extensions are +searched for in this order. If no such file is found, then empty strings are +returned. +""" +function mri_filename(fstring::String, checkdisk::Bool=true) + + fname = "" + fstem = "" + fext = "" + + # List of supported file extensions + extlist = ["mgh", "mgz", "nii", "nii.gz"] + + idot = findlast(isequal('.'), fstring) + + if isnothing(idot) && checkdisk + # Filename has no extension, check if file exists with a supported extension + for ext in extlist + name = fstring * '.' * ext + + if isfile(name) + fname = name + fstem = fstring + fext = ext + end + end + else + # Filename has an extension, check if it is one of the supported ones + ext = lowercase(fstring[(idot+1):end]); + + if cmp(ext, "gz") == 0 + idot = findprev(isequal('.'), fstring, idot-1) + + if !isnothing(idot) + ext = lowercase(fstring[(idot+1):end]) + end + end + + if any(cmp.(ext, extlist) .== 0) + fname = fstring + fstem = fstring[1:idot-1] + fext = ext + end + end + + return fname, fstem, fext +end + + +""" + mri_read(infile::String, headeronly::Bool=false, permuteflag::Bool=false) + +Read a volume from disk and return an `MRI` structure similar to the FreeSurfer +MRI struct defined in mri.h. + +# Arguments +- `infile::String`: Path to the input file, which can be: + 1. An MGH file, e.g., f.mgh or f.mgz + 2. A NIfTI file, e.g., f.nii or f.nii.gz + 3. A file stem, in which case the format and full file name are determined + by finding a file on disk named `infile`.{mgh, mgz, nii, nii.gz} + +- `headeronly::Bool=false`: If true, the pixel data are not read in. + +- `permuteflag::Bool==false`: If true the first two dimensions of the image + volume are permuted in the .vol, .volsize, and .volres fields of the output + structure (this is the default behavior of the MATLAB version). The + `permuteflag` will not affect the vox2ras fields. + +In the output structure, times are in ms and angles are in radians. The +.vox2ras0 field contains the vox2ras matrix that converts 0-based (column, row, +slice) indices to x, y, z coordinates. The .vox2ras1 field contains the same +for 1-based (column, row, slice) indices. + +If the input file is NIfTI, the output `MRI` structure contains a `NIfTIheader` +object in its .niftihdr field. +""" +function mri_read(infile::String, headeronly::Bool=false, permuteflag::Bool=false) + + mri = MRI() + + (fname, fstem, fext) = mri_filename(infile) + if isempty(fname) + error("Cannot determine format of " * infile) + end + + if any(cmp.(fext, ["mgh", "mgz"]) .== 0) #-------- MGH --------# + (mri.vol, M, mr_parms, volsz) = + load_mgh(fname, nothing, nothing, headeronly) + + if isempty(M) + error("Loading " * fname * " as MGH") + end + + if isempty(mr_parms) + tr = flip_angle = te = ti = 0.0 + else + (tr, flip_angle, te, ti) = mr_parms + end + elseif any(cmp.(fext, ["nii", "nii.gz"]) .== 0) #------- NIfTI -------# + hdr = load_nifti(fname, headeronly) + + if !headeronly && isempty(hdr.vol) + error("Loading " * fname * " as NIfTI") + end + + # Compatibility with MRIread.m: + # When data have > 4 dims, put all data into dim 4. + volsz = hdr.dim[2:end] + volsz = volsz[volsz.>0] + volsz = Int.(volsz) + if headeronly + mri.vol = [] + else + mri.vol = hdr.vol + + if length(volsz) > 4 + mri.vol = reshape(mri.vol, volsz[1], volsz[2], volsz[3], :) + end + end + + tr = hdr.pixdim[5] # Already msec + flip_angle = 0 + te = 0 + ti = 0 + hdr.vol = [] # Already have it above, so clear it + M = hdr.vox2ras + + mri.niftihdr = hdr + else + error("File extension " * fext * " not supported") + end + + if permuteflag + mri.vol = permutedims(mri.vol, [2; 1; 3:ndims(mri.vol)]) + mri.ispermuted = true + end + + mri.fspec = fname + mri.pwd = pwd() + + mri.flip_angle = flip_angle + mri.tr = tr + mri.te = te + mri.ti = ti + + # Assumes indices are 0-based. See vox2ras1 below for 1-based. Note: + # mri_write() derives all geometry information (ie, direction cosines, + # voxel resolution, and P0 from vox2ras0. If you change other geometry + # elements of the structure, it will not be reflected in the output + # volume. Also note that vox2ras always maps Col-Row-Slice, even if + # the vol matrix has been permuted and is therefore Row-Col-Slice. + mri.vox2ras0 = M; + + # Dimensions not redundant when using header only + for k in (length(volsz)+1):4 # Make sure all 4 dims are represented + push!(volsz, 1) + end + mri.width = volsz[1] + mri.height = volsz[2] + mri.depth = volsz[3] + mri.nframes = volsz[4] + + if mri.ispermuted + mri.volsize = volsz[[2,1,3]] + else + mri.volsize = volsz[1:3] + end + + #--------------------------------------------------------------------# + # Everything below is redundant in that they can be derived from + # stuff above, but they are in the MRI struct defined in mri.h, so I + # thought I would add them here for completeness. Note: mri_write() + # derives all geometry information (ie, direction cosines, voxel + # resolution, and P0 from vox2ras0. If you change other geometry + # elements below, it will not be reflected in the output volume. + + mri.vox2ras = mri.vox2ras0 + + mri.nvoxels = mri.width * mri.height * mri.depth # Number of voxels + mri.xsize = sqrt(sum(M[:,1].^2)) # Col + mri.ysize = sqrt(sum(M[:,2].^2)) # Row + mri.zsize = sqrt(sum(M[:,3].^2)) # Slice + + mri.x_r = M[1,1]/mri.xsize # Col + mri.x_a = M[2,1]/mri.xsize + mri.x_s = M[3,1]/mri.xsize + + mri.y_r = M[1,2]/mri.ysize # Row + mri.y_a = M[2,2]/mri.ysize + mri.y_s = M[3,2]/mri.ysize + + mri.z_r = M[1,3]/mri.zsize # Slice + mri.z_a = M[2,3]/mri.zsize + mri.z_s = M[3,3]/mri.zsize + + ic = [mri.width/2; mri.height/2; mri.depth/2; 1] + c = M*ic + mri.c_r = c[1] + mri.c_a = c[2] + mri.c_s = c[3] + #--------------------------------------------------------------------# + + #--------------------------------------------------------------------# + # The stuff here is for convenience + + # 1-based vox2ras: Good for doing transforms on julia indices + mri.vox2ras1 = vox2ras_0to1(M) + + # Matrix of direction cosines + mri.Mdc = M[1:3, 1:3] * Diagonal(1 ./ [mri.xsize, mri.ysize, mri.zsize]) + + # Vector of voxel resolutions + if mri.ispermuted # Row-Col-Slice + mri.volres = [mri.ysize, mri.xsize, mri.zsize] + mri.tkrvox2ras = vox2ras_tkreg(mri.volsize[[2,1,3]], mri.volres[[2,1,3]]) + else # Col-Row-Slice + mri.volres = [mri.xsize, mri.ysize, mri.zsize] + mri.tkrvox2ras = vox2ras_tkreg(mri.volsize, mri.volres) + end + + #--------------------------------------------------------------------# + # Optional DWI tables + + bfile = fstem * ".bvals" + + if ~isfile(bfile) + bfile = fstem * ".bval" + + if ~isfile(bfile) + bfile = "" + end + end + + gfile = fstem * ".bvecs" + + if ~isfile(gfile) + gfile = fstem * ".bvec" + + if ~isfile(gfile) + gfile = "" + end + end + + if ~isempty(bfile) && ~isempty(gfile) + (b, g) = mri_read_bfiles(bfile, gfile) + + if length(b) == mri.nframes + mri.bval = b + mri.bvec = g + end + end + #--------------------------------------------------------------------# + + return mri +end + +""" + load_mgh(fname::String, slices::Union{Vector{Unsigned}, Nothing}=nothing, frames::Union{Vector{Unsigned}, Nothing}=nothing, headeronly::Bool=false) + +Load a .mgh or .mgz file from disk. + +Return: +- The image data as an array + +- The 4x4 vox2ras matrix that transforms 0-based voxel indices to coordinates + such that: [x y z]' = M * [i j k 1]' + +- MRI parameters as a vector: [tr, flip_angle, te, ti] + +- The image volume dimensions as a vector (useful when `headeronly` is true, + in which case the image array will be empty) + +# Arguments +- `fname::String`: Path to the .mgh/.mgz file + +- `slices::Vector`: 1-based slice numbers to load (default: read all slices) + +- `frames::Vector`: 1-based frame numbers to load (default: read all frames) + +- `headeronly::Bool=false`: If true, the pixel data are not read in. +""" +function load_mgh(fname::String, slices::Union{Vector{Unsigned}, Nothing}=nothing, frames::Union{Vector{Unsigned}, Nothing}=nothing, headeronly::Bool=false) + + vol = M = mr_parms = volsz = [] + + # Unzip if it is compressed + ext = lowercase(fname[(end-1):end]) + + if cmp(ext, "gz") == 0 + # Create unique temporary file name + tmpfile = tempname(get_tmp_path()) * ".load_mgh.mgh" + gzipped = true + + if Sys.isapple() + cmd = `gunzip -c $fname` + else + cmd = `zcat $fname` + end + run(pipeline(cmd, stdout=tmpfile)) + else + tmpfile = fname + gzipped = false + end + + io = open(tmpfile, "r") + + v = ntoh(read(io, Int32)) + ndim1 = ntoh(read(io, Int32)) + ndim2 = ntoh(read(io, Int32)) + ndim3 = ntoh(read(io, Int32)) + nframes = ntoh(read(io, Int32)) + type = ntoh(read(io, Int32)) + dof = ntoh(read(io, Int32)) + + if !isnothing(slices) && any(slices .> ndim3) + error("Some slices=" * string(slices) * " exceed nslices=" * string(dim3)) + end + + if !isnothing(frames) && any(frames .> nframes) + error("Some frames=" * string(frames) * " exceed nframes=" * string(nframes)) + end + + UNUSED_SPACE_SIZE = 256 + USED_SPACE_SIZE = 3*4+4*3*4 # Space for RAS transform + + unused_space_size = UNUSED_SPACE_SIZE-2 + ras_good_flag = ntoh(read(io, Int16)) + + if ras_good_flag > 0 + delta = ntoh.(read!(io, Vector{Float32}(undef, 3))) + Mdc = ntoh.(read!(io, Vector{Float32}(undef, 9))) + Mdc = reshape(Mdc, (3,3)) + Pxyz_c = ntoh.(read!(io, Vector{Float32}(undef, 3))) + + D = Diagonal(delta) + + Pcrs_c = [ndim1; ndim2; ndim3]/2 # Should this be kept? + + Pxyz_0 = Pxyz_c - Mdc*D*Pcrs_c + + M = [Mdc*D Pxyz_0; 0 0 0 1] + ras_xform = [Mdc Pxyz_c; 0 0 0 1] + unused_space_size = unused_space_size - USED_SPACE_SIZE + end + + skip(io, unused_space_size) + nv = ndim1 * ndim2 * ndim3 * nframes + volsz = [ndim1 ndim2 ndim3 nframes] + + MRI_UCHAR = 0 + MRI_INT = 1 + MRI_LONG = 2 + MRI_FLOAT = 3 + MRI_SHORT = 4 + MRI_BITMAP = 5 + MRI_USHRT = 10 + + # Determine number of bytes per voxel + if type == MRI_FLOAT + nbytespervox = 4 + dtype = Float32 + elseif type == MRI_UCHAR + nbytespervox = 1 + dtype = UInt8 + elseif type == MRI_SHORT + nbytespervox = 2 + dtype = Int16 + elseif type == MRI_USHRT + nbytespervox = 2 + dtype = UInt16 + elseif type == MRI_INT + nbytespervox = 4 + dtype = Int32 + end + + if headeronly + skip(io, nv*nbytespervox) + if !eof(io) + mr_parms = ntoh.(read!(io, Vector{Float32}(undef, 4))) + end + close(io) + + if gzipped # Clean up + cmd = `rm -f $tmpfile` + run(cmd) + end + + return vol, M, mr_parms, volsz + end + + if isnothing(slices) && isnothing(frames) # Read the entire volume + vol = ntoh.(read!(io, Array{dtype}(undef, Tuple(volsz)))) + else # Read a subset of slices/frames + isnothing(frames) && (frames = 1:nframes) + isnothing(slices) && (slices = 1:ndim3) + + nvslice = ndim1 * ndim2 + nvvol = nvslice * ndim3 + filepos0 = position(io) + + vol = zeros(dtype, ndim1, ndim2, length(slices), length(frames)) + + for iframe in 1:length(frames) + frame = frames(iframe) + + for islice in 1:length(slices) + slice = slices(islice) + + filepos = ((frame-1)*nvvol + (slice-1)*nvslice)*nbytespervox + filepos0 + seek(io, filepos) + + vol[:, :, islice, iframe] = + ntoh.(read!(io, Array{dtype}(undef, Tuple(volsz[1:2])))) + end + end + + # Seek to just beyond the last slice/frame + filepos = (nframes*nvvol)*nbytespervox + filepos0; + seek(io, filepos) + end + + if !eof(io) + mr_parms = ntoh.(read!(io, Vector{Float32}(undef, 4))) + end + + close(io) + + if gzipped # Clean up + cmd = `rm -f $tmpfile` + run(cmd) + end + + return vol, M, mr_parms, volsz +end + + +""" + hdr = load_nifti_hdr(fname::String) + +Load the header of a .nii volume from disk. + +Return a `NIfTIheader` structure, where units have been converted to mm and +msec, the sform and qform matrices have been computed and stored in the .sform +and .qform fields, and the .vox2ras field has been set to the sform (if valid), +then the qform (if valid). + +Assume that the input file is uncompressed (compression is handled in the +wrapper load_nifti()). + +Handle data structures with more than 32k cols by setting the .glmin field to +ncols when hdr.dim[2] < 0. This is FreeSurfer specific, for handling surfaces. +""" +function load_nifti_hdr(fname::String) + + hdr = NIfTIheader() + + io = open(fname, "r") + + hdr.sizeof_hdr = read(io, Int32) + hdr.data_type = read!(io, Vector{UInt8}(undef, 10)) + hdr.db_name = read!(io, Vector{UInt8}(undef, 18)) + hdr.extents = read(io, Int32) + hdr.session_error = read(io, Int16) + hdr.regular = read(io, UInt8) + hdr.dim_info = read(io, UInt8) + hdr.dim = read!(io, Vector{Int16}(undef, 8)) + hdr.intent_p1 = read(io, Float32) + hdr.intent_p2 = read(io, Float32) + hdr.intent_p3 = read(io, Float32) + hdr.intent_code = read(io, Int16) + hdr.datatype = read(io, Int16) + hdr.bitpix = read(io, Int16) + hdr.slice_start = read(io, Int16) + hdr.pixdim = read!(io, Vector{Float32}(undef, 8)) + hdr.vox_offset = read(io, Float32) + hdr.scl_slope = read(io, Float32) + hdr.scl_inter = read(io, Float32) + hdr.slice_end = read(io, Int16) + hdr.slice_code = read(io, Int8) + hdr.xyzt_units = read(io, Int8) + hdr.cal_max = read(io, Float32) + hdr.cal_min = read(io, Float32) + hdr.slice_duration = read(io, Float32) + hdr.toffset = read(io, Float32) + hdr.glmax = read(io, Int32) + hdr.glmin = read(io, Int32) + hdr.descrip = read!(io, Vector{UInt8}(undef, 80)) + hdr.aux_file = read!(io, Vector{UInt8}(undef, 24)) + hdr.qform_code = read(io, Int16) + hdr.sform_code = read(io, Int16) + hdr.quatern_b = read(io, Float32) + hdr.quatern_c = read(io, Float32) + hdr.quatern_d = read(io, Float32) + hdr.quatern_x = read(io, Float32) + hdr.quatern_y = read(io, Float32) + hdr.quatern_z = read(io, Float32) + hdr.srow_x = read!(io, Vector{Float32}(undef, 4)) + hdr.srow_y = read!(io, Vector{Float32}(undef, 4)) + hdr.srow_z = read!(io, Vector{Float32}(undef, 4)) + hdr.intent_name = read!(io, Vector{UInt8}(undef, 16)) + hdr.magic = read!(io, Vector{UInt8}(undef, 4)) + + close(io) + + # If numbers have the wrong endian-ness, reverse byte order + if hdr.sizeof_hdr == bswap(Int32(348)) + for var in fieldnames(NIfTIheader) + setfield!(hdr, var, bswap.(getfield(hdr, var))) + end + + hdr.do_bswap = 1 + end + + # This is to accomodate structures with more than 32k cols + # FreeSurfer specific. See also mriio.c. + if hdr.dim[2] < 0 + hdr.dim[2] = hdr.glmin + hdr.glmin = 0 + end + + # Look at xyz units and convert to mm if needed + xyzunits = hdr.xyzt_units & Int8(7) # Bitwise AND with 00000111 + if xyzunits == 1 + xyzscale = 1000.000 # meters + elseif xyzunits == 2 + xyzscale = 1.000 # mm + elseif xyzunits == 3 + xyzscale = .001 # microns + else + println("WARNING: xyz units code " * string(xyzunits) * + " is unrecognized, assuming mm") + xyzscale = 1 # just assume mm + end + + hdr.pixdim[2:4] = hdr.pixdim[2:4] * xyzscale + hdr.srow_x = hdr.srow_x * xyzscale + hdr.srow_y = hdr.srow_y * xyzscale + hdr.srow_z = hdr.srow_z * xyzscale + + # Look at time units and convert to msec if needed + tunits = hdr.xyzt_units & Int8(56) # Bitwise AND with 00111000 + if tunits == 8 + tscale = 1000.000 # seconds + elseif tunits == 16 + tscale = 1.000 # msec + elseif tunits == 32 + tscale = .001 # microsec + else + tscale = 0 # no time scale + end + + hdr.pixdim[5] = hdr.pixdim[5] * tscale + + # Change value in xyzt_units to reflect scale change + hdr.xyzt_units = Int8(2) | Int8(16) # Bitwise OR of 2=mm, 16=msec + + # Sform matrix + hdr.sform = [hdr.srow_x'; + hdr.srow_y'; + hdr.srow_z'; + 0 0 0 1] + + # Qform matrix + # (From DG's load_nifti_hdr.m: not quite sure how all this works, + # mainly just copied CH's code from mriio.c) + b = hdr.quatern_b + c = hdr.quatern_c + d = hdr.quatern_d + x = hdr.quatern_x + y = hdr.quatern_y + z = hdr.quatern_z + a = 1.0 - (b*b + c*c + d*d) + if abs(a) < 1.0e-7 + a = 1.0 / sqrt(b*b + c*c + d*d) + b = b*a + c = c*a + d = d*a + a = 0.0 + else + a = sqrt(a) + end + r11 = a*a + b*b - c*c - d*d + r12 = 2.0*b*c - 2.0*a*d + r13 = 2.0*b*d + 2.0*a*c + r21 = 2.0*b*c + 2.0*a*d + r22 = a*a + c*c - b*b - d*d + r23 = 2.0*c*d - 2.0*a*b + r31 = 2.0*b*d - 2*a*c + r32 = 2.0*c*d + 2*a*b + r33 = a*a + d*d - c*c - b*b + if hdr.pixdim[1] < 0.0 + r13 = -r13 + r23 = -r23 + r33 = -r33 + end + qMdc = [r11 r12 r13; r21 r22 r23; r31 r32 r33] + D = Diagonal(hdr.pixdim[2:4]) + P0 = [x y z]' + hdr.qform = [qMdc*D P0; 0 0 0 1] + + if hdr.sform_code != 0 + # Use sform first + hdr.vox2ras = hdr.sform + elseif hdr.qform_code != 0 + # Then use qform first + hdr.vox2ras = hdr.qform + else + println("WARNING: neither sform or qform are valid in " * fname) + D = Diagonal(hdr.pixdim[2:4]) + P0 = [0 0 0]' + hdr.vox2ras = [D P0; 0 0 0 1] + end + + return hdr +end + + +""" + load_nifti(fname::String, hdronly::Bool=false) + +Load a NIfTI (.nii or .nii.gz) volume from disk and return a `NIfTIheader` +structure. + +Handle compressed NIfTI (nii.gz) by issuing an external Unix call to +uncompress the file to a temporary file, which is then deleted. + +The output structure contains: +- the image data, in the .vol field +- the units for each dimension of the volume [mm or msec], in the .pixdim field +- the sform and qform matrices, in the .sform and .qform fields +- the vox2ras matrix, which is the sform (if valid), otherwise the qform, in + the .vox2ras field + +Handle data structures with more than 32k cols by looking at the .dim field. +If dim[2] = -1, then the .glmin field contains the numbfer of columns. This is +FreeSurfer specific, for handling surfaces. When the total number of spatial +voxels equals 163842, then reshape the volume to 163842x1x1xnframes. This is +for handling the 7th order icosahedron used by FS group analysis. +""" +function load_nifti(fname::String, hdronly::Bool=false) + + # Unzip if it is compressed + ext = lowercase(fname[(end-1):end]) + + if cmp(ext, "gz") == 0 + # Create unique temporary file name + tmpfile = tempname(get_tmp_path()) * ".load_nifti.nii" + gzipped = true + + if Sys.isapple() + cmd = `gunzip -c $fname` + else + cmd = `zcat $fname` + end + run(pipeline(cmd, stdout=tmpfile)) + else + tmpfile = fname + gzipped = false + end + + # Read NIfTI header + hdr = load_nifti_hdr(tmpfile) + + if isempty(hdr.vox2ras) + if gzipped # Clean up + cmd = `rm -f $tmpfile` + run(cmd) + end + return hdr + end + + # Check for ico7 + nspatial = prod(Int32.(hdr.dim[2:4])) + IsIco7 = (nspatial == 163842) + + # If only header is desired, return now + if hdronly + if gzipped # Clean up + cmd = `rm -f $tmpfile` + run(cmd) + end + if IsIco7 # Reshape + hdr.dim[2] = 163842 + hdr.dim[3] = 1 + hdr.dim[4] = 1 + end + return hdr + end + + # Get volume dimensions + dim = hdr.dim[2:end] + while (dim[end] == 0) && !isempty(dim) + pop!(dim) + end + + # Open to read the pixel data + io = open(tmpfile, "r") + + # Get past the header + seek(io, Int64(round(hdr.vox_offset))) + + if hdr.datatype == 2 + dtype = UInt8 + elseif hdr.datatype == 4 + dtype = Int16 + elseif hdr.datatype == 8 + dtype = Int32 + elseif hdr.datatype == 16 + dtype = Float32 + elseif hdr.datatype == 64 + dtype = Float64 + elseif hdr.datatype == 256 + dtype = Int8 + elseif hdr.datatype == 512 + dtype = UInt16 + elseif hdr.datatype == 768 + dtype = UInt32 + else + close(io) + if gzipped # Clean up + cmd = `rm -f $tmpfile` + run(cmd) + end + error("Data type " * string(hdr.datatype) * " not supported") + end + + hdr.vol = read!(io, Array{dtype}(undef, Tuple(dim))) + + close(io) + if gzipped # Clean up + cmd = `rm -f $tmpfile` + run(cmd) + end + + # Check if end-of-file was reached + if !eof(io) + error(tmpfile * ", read a " * string(size(hdr.vol)) * + " volume but did not reach end of file") + end + + # If needed, reverse order of bytes to correct endian-ness + if hdr.do_bswap == 1 + hdr.vol = bswap.(hdr.vol) + end + + if IsIco7 + hdr.dim[2] = 163842 + hdr.dim[3] = 1 + hdr.dim[4] = 1 + dim = hdr.dim[2:end] + + hdr.vol = reshape(hdr.vol, Tuple(dim)) + end + + if hdr.scl_slope!=0 && !(hdr.scl_inter==0 && hdr.scl_slope==1) + # Rescaling is not needed if the slope==1 and intersect==0, + # skipping this preserves the numeric class of the data + hdr.vol = Float64(hdr.vol) .* hdr.scl_slope .+ hdr.scl_inter + end + + return hdr +end + + +""" + mri_write(mri::MRI, outfile::String, datatype::DataType=Float32) + +Write an MRI volume to disk. Return true is an error occurred (i.e., the +number of bytes written were not as expected based on the size of the volume). + +# Arguments +- `mri::MRI`: A structure like that returned by `mri_read()`. The geometry + (i.e., direction cosines, voxel resolution, and P0) are all recomputed from + mri.vox2ras0. So, if a method has changed one of the other fields, e.g., + mri.x_r, this change will not be reflected in the output volume. + +- `outfile::String`: Path to the output file, which can be: + 1. An MGH file, e.g., f.mgh or f.mgz (uncompressed or compressed) + 2. A NIfTI file, e.g., f.nii or f.nii.gz (uncompressed or compressed). + +- `datatype::DataType=Float32`: Only applies to NIfTI and can be UInt8, Int16, + Int32, Float32, Float64, Int8, UInt16, UInt32. +""" +function mri_write(mri::MRI, outfile::String, datatype::DataType=Float32) + + err = true + + if isempty(mri.vol) + error("Input structure has empty vol field") + end + + vsz = collect(size(mri.vol)) + nvsz = length(vsz) + if nvsz < 4 + vsz = [vsz; ones(eltype(vsz), 4-nvsz)] + end + + if isempty(mri.volsize) + mri.volsize = vsz[1:3] + end + + if mri.nframes == 0 + mri.nframes = vsz[4] + end + + if isempty(mri.vox2ras0) + mri.vox2ras0 = Diagonal(ones(4)) + end + + if isempty(mri.volres) + mri.volres = sqrt.((ones(3)' * (mri.vox2ras0[1:3,1:3].^2))') + end + + (fname, fstem, fext) = mri_filename(outfile, false) # false = no checkdisk + if isempty(fname) + error("Cannot determine format of " * outfile) + end + + if any(cmp.(fext, ["mgh", "mgz"]) .== 0) #-------- MGH --------# + M = mri.vox2ras0 + mr_parms = [mri.tr, mri.flip_angle, mri.te, mri.ti] + + if mri.ispermuted + err = save_mgh(permutedims(mri.vol, [2; 1; 3:ndims(mri.vol)]), + fname, M, mr_parms) + else + err = save_mgh(mri.vol, fname, M, mr_parms) + end + elseif any(cmp.(fext, ["nii", "nii.gz"]) .== 0) #------- NIfTI -------# + hdr = NIfTIheader() + + hdr.db_name = zeros(UInt8, 18) + hdr.data_type = zeros(UInt8, 10) + + if mri.ispermuted + hdr.dim = vcat(mri.volsize[[2,1,3]], mri.nframes) + else + hdr.dim = vcat(mri.volsize[1:3], mri.nframes) + end + + if datatype == UInt8 + hdr.datatype = 2 + hdr.bitpix = 8 + elseif datatype == Int16 + hdr.datatype = 4 + hdr.bitpix = 16 + elseif datatype == Int32 + hdr.datatype = 8 + hdr.bitpix = 32 + elseif datatype == Float32 + hdr.datatype = 16 + hdr.bitpix = 32 + elseif datatype == Float64 + hdr.datatype = 64 + hdr.bitpix = 64 + elseif datatype == Int8 + hdr.datatype = 256 + hdr.bitpix = 8 + elseif datatype == UInt16 + hdr.datatype = 512 + hdr.bitpix = 16 + elseif datatype == UInt32 + hdr.datatype = 768 + hdr.bitpix = 32 + else + error("Data type " * string(datatype) * " not supported") + end + + if mri.ispermuted + hdr.pixdim = vcat(0, mri.volres[[2,1,3]], mri.tr) # Physical units + else + hdr.pixdim = vcat(0, mri.volres[1:3], mri.tr) # Physical units + end + + hdr.vox_offset = 348 # Will be set again + hdr.scl_slope = mri.niftihdr.scl_slope + hdr.scl_inter = mri.niftihdr.scl_inter + + hdr.xyzt_units = Int8(2) | Int8(16) # Bitwise OR of 2=mm, 16=msec + + hdr.cal_max = maximum(mri.vol) + hdr.cal_min = minimum(mri.vol) + hdr.descrip = collect(@sprintf("%-80s","FreeSurfer julia")) + hdr.aux_file = [UInt8(0)] + hdr.qform_code = 1 # 1=NIFTI_XFORM_SCANNER_ANAT + hdr.sform_code = 1 # 1=NIFTI_XFORM_SCANNER_ANAT + + # Qform (must have 6 DOF) + (b, c, d, x, y, z, qfac) = vox2ras_to_qform(mri.vox2ras0) + hdr.pixdim[1] = qfac + hdr.quatern_b = b + hdr.quatern_c = c + hdr.quatern_d = d + hdr.quatern_x = x + hdr.quatern_y = y + hdr.quatern_z = z + # Sform (can be any affine) + hdr.srow_x = mri.vox2ras0[1,:] + hdr.srow_y = mri.vox2ras0[2,:] + hdr.srow_z = mri.vox2ras0[3,:] + hdr.intent_name = collect("huh?") + hdr.magic = collect("n+1") + + if mri.ispermuted + hdr.vol = permutedims(mri.vol, [2; 1; 3:ndims(mri.vol)]) + else + hdr.vol = mri.vol + end + + err = save_nifti(hdr, fname) + else + error("File extension " * fext * " not supported") + end + + if err + println("WARNING: Problem saving " * outfile) + end + + return err +end + + +""" + save_mgh(vol::Array, fname::String, M::Matrix=Diagonal(ones(4)), mr_parms::Vector=zeros(4)) + +Write an MRI volume to a .mgh or .mgz file. Return true is an error occurred +(i.e., the number of bytes written were not as expected based on the size of +the volume). + +# Arguments +- `vol::Array`: the image data + +- `fname::String`: path to the output file + +- `M::Matrix`: the 4x4 vox2ras transform such that [x y z]' = M * [i j k 1]', + where the voxel indices (i, j, k) are 0-based. + +- `mr_parms::Vector`: a vector of MRI parameters, [tr, flip_angle, te, ti] +""" +function save_mgh(vol::Array, fname::String, M::Matrix=Diagonal(ones(4)), mr_parms::Vector=zeros(4)) + + if size(M) != (4, 4) + error("M size=" * string(size(M)) * ", must be (4, 4)") + end + + if length(mr_parms) != 4 + error("mr_parms length=" * string(length(mr_parms)) * ", must be 4") + end + + # Not all of these appear to be used + MRI_UCHAR = 0 + MRI_INT = 1 + MRI_LONG = 2 + MRI_FLOAT = 3 + MRI_SHORT = 4 + MRI_BITMAP = 5 + MRI_TENSOR = 6 + + io = open(fname, "w") + + (ndim1, ndim2, ndim3, frames) = size(vol) + nb = 0 + + # Write everything as big-endian + nb += write(io, hton(Int32(1))) # magic number + nb += write(io, hton(Int32(ndim1))) + nb += write(io, hton(Int32(ndim2))) + nb += write(io, hton(Int32(ndim3))) + nb += write(io, hton(Int32(frames))) + if ndims(vol) == 5 + is_tensor = true + nb += write(io, hton(Int32(MRI_TENSOR))) # type = MRI_TENSOR + else + is_tensor = false + nb += write(io, hton(Int32(MRI_FLOAT))) # type = MRI_FLOAT + end + + nb += write(io, hton(Int32(1))) # dof (not used) + + UNUSED_SPACE_SIZE = 256 + USED_SPACE_SIZE = (3*4+4*3*4) # Space for RAS transform + + MdcD = M[1:3,1:3] + + delta = sqrt.(sum(MdcD.^2; dims=1)) + Mdc = MdcD ./ repeat(delta, 3) + + Pcrs_c = [ndim1/2, ndim2/2, ndim3/2, 1] + Pxyz_c = M*Pcrs_c + Pxyz_c = Pxyz_c[1:3] + + nb += write(io, hton(Int16(1))) # ras_good_flag = 1 + nb += write(io, hton.(Float32.(delta))) + nb += write(io, hton.(Float32.(Mdc))) + nb += write(io, hton.(Float32.(Pxyz_c))) + + unused_space_size = UNUSED_SPACE_SIZE-2 + unused_space_size = unused_space_size - USED_SPACE_SIZE + nb += write(io, UInt8.(zeros(unused_space_size))) + + nb += write(io, hton.(Float32.(vol))) + + nb += write(io, hton.(Float32.(mr_parms))) + + close(io) + + err = (nb != (sizeof(Int32) * 7 + + sizeof(Int16) + + sizeof(UInt8) * unused_space_size + + sizeof(Float32) * (19 + length(vol)))) + + ext = lowercase(fname[(end-1):end]) + + if cmp(ext, "gz") == 0 + cmd = `gzip -f $fname` + run(cmd) + cmd = `mv $fname.gz $fname` + run(cmd) + end + + return err +end + + +""" + save_nifti(hdr::NIfTIheader, fname::String) + +Write an MRI volume to a .nii or .nii.gz file. Return true is an error occurred +(i.e., the number of bytes written were not as expected based on the size of +the volume). + +# Arguments +- `hdr::NIfTIheader`: a NIfTI header structure that contains the image data in + its .vol field + +- `fname::String`: path to the output file + +Handle data structures with more than 32k cols by setting hdr.dim[2] = -1 and +hdr.glmin = ncols. This is FreeSurfer specific, for handling surfaces. The +exception to this is when the total number of spatial voxels equals 163842, +then the volume is reshaped to 27307x1x6xnframes. This is for handling the 7th +order icosahedron used by FS group analysis. +""" +function save_nifti(hdr::NIfTIheader, fname::String) + + ext = lowercase(fname[(end-1):end]) + if cmp(ext, "gz") == 0 + gzip_needed = true + fname = fname[1:(end-3)] + else + gzip_needed = false + end + + # Check for ico7 + sz = size(hdr.vol) + if sz[1] == 163842 + dim = (27307, 1, 6, size(hdr.vol,4)) + hdr.vol = reshape(hdr.vol, dim) + end + + io = open(fname, "w") + + hdr.data_type = vcat(hdr.data_type, + repeat([UInt8(0)], 10-length(hdr.data_type))) + hdr.db_name = vcat(hdr.db_name, + repeat([UInt8(0)], 18-length(hdr.db_name))) + + hdr.dim = ones(8) + if size(hdr.vol, 4) > 1 + hdr.dim[1] = 4 + else + hdr.dim[1] = 3 + end + for k in (2:5) + hdr.dim[k] = size(hdr.vol, k-1) + end + + # This is to accomodate structures with more than 32k cols + # FreeSurfer specific. See also mriio.c. + if hdr.dim[2] > 2^15 + hdr.glmin = hdr.dim[2] + hdr.dim[2] = -1 + end + + hdr.pixdim = vcat(hdr.pixdim, zeros(8-length(hdr.pixdim))) + + hdr.descrip = vcat(hdr.descrip, + repeat([UInt8(0)], 80-length(hdr.descrip))) + + hdr.aux_file = vcat(hdr.aux_file, + repeat([UInt8(0)], 24-length(hdr.aux_file))) + + hdr.intent_name = vcat(hdr.intent_name, + repeat([UInt8(0)], 16-length(hdr.intent_name))) + + hdr.magic = vcat(hdr.magic, + repeat([UInt8(0)], 4-length(hdr.magic))) + + hdr.vox_offset = 352 # not 348 + + nb = 0 + + nb += write(io, Int32(348)) + nb += write(io, UInt8.(hdr.data_type)) + nb += write(io, UInt8.(hdr.db_name)) + nb += write(io, Int32(hdr.extents)) + nb += write(io, Int16(hdr.session_error)) + nb += write(io, UInt8(hdr.regular)) + nb += write(io, UInt8(hdr.dim_info)) + nb += write(io, Int16.(hdr.dim)) + nb += write(io, Float32(hdr.intent_p1)) + nb += write(io, Float32(hdr.intent_p2)) + nb += write(io, Float32(hdr.intent_p3)) + nb += write(io, Int16(hdr.intent_code)) + nb += write(io, Int16(hdr.datatype)) + nb += write(io, Int16(hdr.bitpix)) + nb += write(io, Int16(hdr.slice_start)) + nb += write(io, Float32.(hdr.pixdim)) + nb += write(io, Float32(hdr.vox_offset)) + nb += write(io, Float32(hdr.scl_slope)) + nb += write(io, Float32(hdr.scl_inter)) + nb += write(io, Int16(hdr.slice_end)) + nb += write(io, Int8(hdr.slice_code)) + nb += write(io, Int8(hdr.xyzt_units)) + nb += write(io, Float32(hdr.cal_max)) + nb += write(io, Float32(hdr.cal_min)) + nb += write(io, Float32(hdr.slice_duration)) + nb += write(io, Float32(hdr.toffset)) + nb += write(io, Int32(hdr.glmax)) + nb += write(io, Int32(hdr.glmin)) + nb += write(io, UInt8.(hdr.descrip)) + nb += write(io, UInt8.(hdr.aux_file)) + nb += write(io, Int16(hdr.qform_code)) + nb += write(io, Int16(hdr.sform_code)) + nb += write(io, Float32(hdr.quatern_b)) + nb += write(io, Float32(hdr.quatern_c)) + nb += write(io, Float32(hdr.quatern_d)) + nb += write(io, Float32(hdr.quatern_x)) + nb += write(io, Float32(hdr.quatern_y)) + nb += write(io, Float32(hdr.quatern_z)) + nb += write(io, Float32.(hdr.srow_x)) + nb += write(io, Float32.(hdr.srow_y)) + nb += write(io, Float32.(hdr.srow_z)) + nb += write(io, UInt8.(hdr.intent_name)) + nb += write(io, UInt8.(hdr.magic)) + + # Pad to get to 352 bytes (header size is 348) + nb += write(io, UInt8.(zeros(4))) + + if hdr.datatype == 2 + dtype = UInt8 + elseif hdr.datatype == 4 + dtype = Int16 + elseif hdr.datatype == 8 + dtype = Int32 + elseif hdr.datatype == 16 + dtype = Float32 + elseif hdr.datatype == 64 + dtype = Float64 + elseif hdr.datatype == 256 + dtype = Int8 + elseif hdr.datatype == 512 + dtype = UInt16 + elseif hdr.datatype == 768 + dtype = UInt32 + else + println("WARNING: data type " * string(hdr.datatype) * + " not supported, but writing as float") + dtype = Float32 + end + + nb += write(io, dtype.(hdr.vol)) + + close(io) + + err = (nb != (sizeof(Float32) * 36 + + sizeof(Int32) * 4 + + sizeof(Int16) * 16 + + sizeof(Int8) * 2 + + sizeof(UInt8) * 158 + + sizeof(dtype) * length(hdr.vol))) + + if gzip_needed + cmd = `gzip -f $fname` + run(cmd) + end + + return err +end + + +""" + mri_read_bfiles(infile1::String, infile2::String) + +Read a DWI b-value table and gradient table from text files `infile1` and +`infile2`. The two input files can be specified in any order. The gradient +table file must contain 3 times as many entries as the b-value table file. + +Return the b-value table as a vector of size n and the gradient table as a +matrix of size (n, 3). +""" +function mri_read_bfiles(infile1::String, infile2::String) + + tab = [] + + for infile in (infile1, infile2) + if ~isfile(infile) + error("Could not open " * infile) + end + + push!(tab, readdlm(infile)) + + if ~all(isa.(tab[end], Number)) + error("File " * infile * " contains non-numeric entries") + end + + if size(tab[end], 2) > size(tab[end], 1) + tab[end] = permutedims(tab[end], 2:1) + end + + if size(tab[end], 2) == 1 + tab[end] = tab[end][:,1] + end + end + + if size(tab[1], 1) != size(tab[2], 1) + error("Dimension mismatch between tables in " * + infile1 * " " * string(size(tab[1])) * " and " * + infile2 * " " * string(size(tab[2]))) + end + + if (size(tab[1], 2) == 1 && size(tab[2], 2) == 3) || + (size(tab[1], 2) == 3 && size(tab[2], 2) == 1) + return tab[1], tab[2] + else + error("Wrong number of entries in tables " * + infile1 * " " * string(size(tab[1])) * " and " * + infile2 * " " * string(size(tab[2]))) + end +end + + +""" + mri_read_bfiles!(dwi::MRI, infile1::String, infile2::String) + +Set the .bval and .bvec fields of the MRI structure `dwi`, by reading a DWI +b-value table and gradient table from the text files `infile1` and `infile2`. +The two input files can be specified in any order. The gradient table file +must contain 3 times as many entries as the b-value table file. + +Return the b-value table as a vector of size n and the gradient table as a +matrix of size (n, 3). +""" +function mri_read_bfiles!(dwi::MRI, infile1::String, infile2::String) + + (tab1, tab2) = mri_read_bfiles(infile1, infile2) + + if size(tab1, 1) != size(dwi.vol, 4) + error("Number of frames in volume (" * string(size(dwi.vol, 4)) * + ") does not match dimensions of table in " * + infile1 * " " * string(size(tab1))) + end + + if size(tab1, 2) == 1 + dwi.bval = tab1 + dwi.bvec = tab2 + else + dwi.bval = tab2 + dwi.bvec = tab1 + end + + return tab1, tab2 +end + + diff --git a/julia/FreeSurfer/src/trk.jl b/julia/FreeSurfer/src/trk.jl new file mode 100644 index 00000000000..6bf6d092542 --- /dev/null +++ b/julia/FreeSurfer/src/trk.jl @@ -0,0 +1,311 @@ +#= + Original Author: Anastasia Yendiki + + Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH" + + Terms and conditions for use, reproduction, distribution and contribution + are found in the 'FreeSurfer Software License Agreement' contained + in the file 'LICENSE' found in the FreeSurfer distribution, and here: + + https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense + + Reporting: freesurfer@nmr.mgh.harvard.edu +=# + +using LinearAlgebra + +export Tract, trk_read, trk_write + + +"Container for header and streamline data stored in .trk format" +mutable struct Tract + # Header fields (.trk format version 2) + id_string::Vector{UInt8} + dim::Vector{Int16} + voxel_size::Vector{Float32} + origin::Vector{Float32} + n_scalars::Int16 + scalar_name::Matrix{UInt8} + n_properties::Int16 + property_name::Matrix{UInt8} + vox_to_ras::Matrix{Float32} + reserved::Vector{UInt8} + voxel_order::Vector{UInt8} + voxel_order_original::Vector{UInt8} + image_orientation_patient::Vector{Float32} + pad1::Vector{UInt8} + invert_x::UInt8 + invert_y::UInt8 + invert_z::UInt8 + swap_xy::UInt8 + swap_yz::UInt8 + swap_zx::UInt8 + n_count::Int32 + version::Int32 + hdr_size::Int32 + + # Streamline data fields + npts::Vector{Int32} + properties::Matrix{Float32} + xyz::Vector{Any} + scalars::Vector{Any} +end + + +""" + Tract() + +Return an empty `Tract` structure +""" +Tract() = Tract( + Vector{UInt8}(undef, 0), + Vector{Int16}(undef, 0), + Vector{Float32}(undef, 0), + Vector{Float32}(undef, 0), + Int16(0), + Matrix{UInt8}(undef, 0, 0), + Int16(0), + Matrix{UInt8}(undef, 0, 0), + Matrix{Float32}(undef, 0, 0), + Vector{UInt8}(undef, 0), + Vector{UInt8}(undef, 0), + Vector{UInt8}(undef, 0), + Vector{Float32}(undef, 0), + Vector{UInt8}(undef, 0), + UInt8(0), + UInt8(0), + UInt8(0), + UInt8(0), + UInt8(0), + UInt8(0), + Int32(0), + Int32(0), + Int32(0), + + Vector{Int32}(undef, 0), + Matrix{Float32}(undef, 0, 0), + [], + [] +) + + +""" + Tract(ref::MRI) + +Return a `Tract` structure whose header fields are populated based on the +reference `MRI` structure `ref` +""" +function Tract(ref::MRI) + + tr = Tract() + + # In the following I must take into account that mri_read() + # has reversed the first 2 dimensions of volsize and volres, + # but it has not changed the vox2ras matrix + + # Find orientation of image coordinate system + orient = Vector{Char}(undef, 3) + + for idim in 1:3 + (amax, imax) = findmax(abs.(ref.vox2ras[idim, 1:3])) + if imax == 1 + if ref.vox2ras[idim, imax] > 0 + orient[idim] = 'R' + else + orient[idim] = 'L' + end + elseif imax == 2 + if ref.vox2ras[idim, imax] > 0 + orient[idim] = 'A' + else + orient[idim] = 'P' + end + else + if ref.vox2ras[idim, imax] > 0 + orient[idim] = 'S' + else + orient[idim] = 'I' + end + end + end + + # Find patient-to-scanner coordinate transform: + # Take x and y vectors from vox2RAS matrix, convert to LPS, + # divide by voxel size + p2s = [-1 0 0; 0 -1 0; 0 0 1] * ref.vox2ras[1:3, 1:2] * + Diagonal([1, 1]./ref.volres[2,1]) + + tr.id_string = UInt8.(collect("TRACK\0")) + + tr.dim = Int16.(ref.volsize[[2,1,3]]) + tr.voxel_size = Float32.(ref.volres[[2,1,3]]) + tr.origin = Float32.(zeros(3)) + + tr.n_scalars = Int16(0) + tr.scalar_name = UInt8.(zeros(10, 20)) + tr.n_properties = Int16(0) + tr.property_name = UInt8.(zeros(10, 20)) + + tr.vox_to_ras = ref.vox2ras + tr.reserved = UInt8.(zeros(444)) + tr.voxel_order = vcat(UInt8.(collect(orient)), UInt8(0)) + tr.voxel_order_original = tr.voxel_order + tr.image_orientation_patient = Float32.(p2s[:]) + tr.pad1 = UInt8.(zeros(2)) + + tr.invert_x = UInt8(0) + tr.invert_y = UInt8(0) + tr.invert_z = UInt8(0) + tr.swap_xy = UInt8(0) + tr.swap_yz = UInt8(0) + tr.swap_zx = UInt8(0) + tr.n_count = Int32(0) + tr.version = Int32(2) + tr.hdr_size = Int32(1000) + + return tr +end + + +""" + trk_read(infile::String) + +Read tractography streamlines from `infile` and return a `Tract` structure + +Input file must be in .trk format, see: +http://trackvis.org/docs/?subsect=fileformat +""" +function trk_read(infile::String) + + tr = Tract() + + io = open(infile, "r") + + # Read .trk header + tr.id_string = read!(io, Vector{UInt8}(undef, 6)) + + tr.dim = read!(io, Vector{Int16}(undef, 3)) + tr.voxel_size = read!(io, Vector{Float32}(undef, 3)) + tr.origin = read!(io, Vector{Float32}(undef, 3)) + + tr.n_scalars = read(io, Int16) + tr.scalar_name = read!(io, Matrix{UInt8}(undef, 20, 10)) + tr.scalar_name = permutedims(tr.scalar_name, [2,1]) + + tr.n_properties = read(io, Int16) + tr.property_name = read!(io, Matrix{UInt8}(undef, 20, 10)) + tr.property_name = permutedims(tr.property_name, [2,1]) + + tr.vox_to_ras = read!(io, Matrix{Float32}(undef, 4, 4)) + tr.vox_to_ras = permutedims(tr.vox_to_ras, [2,1]) + tr.reserved = read!(io, Vector{UInt8}(undef, 444)) + tr.voxel_order = read!(io, Vector{UInt8}(undef, 4)) + tr.voxel_order_original = read!(io, Vector{UInt8}(undef, 4)) + tr.image_orientation_patient = read!(io, Vector{Float32}(undef, 6)) + tr.pad1 = read!(io, Vector{UInt8}(undef, 2)) + + tr.invert_x = read(io, UInt8) + tr.invert_y = read(io, UInt8) + tr.invert_z = read(io, UInt8) + tr.swap_xy = read(io, UInt8) + tr.swap_yz = read(io, UInt8) + tr.swap_zx = read(io, UInt8) + tr.n_count = read(io, Int32) + tr.version = read(io, Int32) + tr.hdr_size = read(io, Int32) + + # Read streamline data + tr.npts = Vector{Int32}(undef, tr.n_count) + tr.properties = Matrix{Float32}(undef, tr.n_properties, tr.n_count) + + for istr in 1:tr.n_count # Loop over streamlines + tr.npts[istr] = read(io, Int32) + + push!(tr.xyz, Matrix{Float32}(undef, 3, tr.npts[istr])) + push!(tr.scalars, Matrix{Float32}(undef, tr.n_scalars, tr.npts[istr])) + + for ipt in 1:tr.npts[istr] # Loop over points on a streamline + # Divide by voxel size and make voxel coordinates 0-based + tr.xyz[istr][:, ipt] = read!(io, Vector{Float32}(undef, 3)) ./ + tr.voxel_size .- .5 + + tr.scalars[istr][:, ipt] = read!(io, Vector{Float32}(undef, tr.n_scalars)) + end + + tr.properties[:, istr] = read!(io, Vector{Float32}(undef, tr.n_properties)) + end + + close(io) + + return tr +end + +""" + trk_write(tr::Tract, outfile::String) + +Write a `Tract` structure to a file in the .trk format + +Return true if an error occurred (i.e., the number of bytes written was not the +expected based on the size of the `Tract` structure) +""" +function trk_write(tr::Tract, outfile::String) + + io = open(outfile, "w") + + nb = 0 + + # Write .trk header + nb += write(io, UInt8.(tr.id_string)) + + nb += write(io, Int16.(tr.dim)) + nb += write(io, Float32.(tr.voxel_size)) + nb += write(io, Float32.(tr.origin)) + + nb += write(io, Int16(tr.n_scalars)) + nb += write(io, UInt8.(permutedims(tr.scalar_name, [2,1]))) + nb += write(io, Int16(tr.n_properties)) + nb += write(io, UInt8.(permutedims(tr.property_name, [2,1]))) + + nb += write(io, Float32.(permutedims(tr.vox_to_ras, [2,1]))) + nb += write(io, UInt8.(tr.reserved)) + nb += write(io, UInt8.(tr.voxel_order)) + nb += write(io, UInt8.(tr.voxel_order_original)) + nb += write(io, Float32.(tr.image_orientation_patient)) + nb += write(io, UInt8.(tr.pad1)) + + nb += write(io, UInt8(tr.invert_x)) + nb += write(io, UInt8(tr.invert_y)) + nb += write(io, UInt8(tr.invert_z)) + nb += write(io, UInt8(tr.swap_xy)) + nb += write(io, UInt8(tr.swap_yz)) + nb += write(io, UInt8(tr.swap_zx)) + nb += write(io, Int32(tr.n_count)) + nb += write(io, Int32(tr.version)) + nb += write(io, Int32(tr.hdr_size)) + + # Write streamline data + for istr in 1:tr.n_count # Loop over streamlines + nb += write(io, Int32(tr.npts[istr])) + + for ipt in 1:tr.npts[istr] # Loop over points on a streamline + # Make voxel coordinates .5-based and multiply by voxel size + nb += write(io, Float32.((tr.xyz[istr][:, ipt] .+ .5) .* tr.voxel_size)) + + nb += write(io, Float32.(tr.scalars[istr][:, ipt])) + end + + nb += write(io, Float32.(tr.properties[:, istr])) + end + + close(io) + + err = (nb != sizeof(UInt8) * 866 + + sizeof(Int32) * (3 + length(tr.npts)) + + sizeof(Int16) * 5 + + sizeof(Float32) * (28 + sum(length.(tr.xyz)) + + sum(length.(tr.scalars)) + + length(tr.properties))) + + return err +end + diff --git a/julia/FreeSurfer/src/view.jl b/julia/FreeSurfer/src/view.jl new file mode 100644 index 00000000000..39186472591 --- /dev/null +++ b/julia/FreeSurfer/src/view.jl @@ -0,0 +1,229 @@ +#= + Original Author: Anastasia Yendiki + + Copyright © 2022 The General Hospital Corporation (Boston, MA) "MGH" + + Terms and conditions for use, reproduction, distribution and contribution + are found in the 'FreeSurfer Software License Agreement' contained + in the file 'LICENSE' found in the FreeSurfer distribution, and here: + + https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense + + Reporting: freesurfer@nmr.mgh.harvard.edu +=# + +using ColorTypes, DelimitedFiles, ImageInTerminal + +export LUT, color_lut, show, view + +"Container for segmentation and tract look-up tables" +mutable struct LUT + id::Vector{Int} + name::Vector{String} + rgb::Vector{RGB} +end + + +""" + LUT() + +Return an empty `LUT` structure +""" +LUT() = LUT( + Vector{Int}(undef, 0), + Vector{String}(undef, 0), + Vector{RGB}(undef, 0) +) + + +""" + LUT(infile::String) + +Read a look-up table from `infile` and return a `LUT` structure + +The input file is assumed to have the format of FreeSurferColorLUT.txt +""" +function LUT(infile::String) + + lut = LUT() + + if ~isfile(infile) + error(infile * "is not a regular file") + end + + tab = readdlm(infile; comments=true, comment_char='#') + + lut.id = tab[:,1] + lut.name = tab[:,2] + lut.rgb = RGB.(tab[:,3]/255, tab[:,4]/255, tab[:,5]/255) + + return lut +end + + +"The FreeSurfer color look-up table" +const global color_lut = haskey(ENV, "FREESURFER_HOME") ? + LUT(ENV["FREESURFER_HOME"]*"/FreeSurferColorLUT.txt") : + LUT() + + +""" + vol_to_rgb(vol::Array) + +Convert an image array to an RGB/Gray array for display. + +Determine how the image should be displayed: +- If all the values are IDs in the FreeSurfer color look-up table, + the image is assumed to be is a segmentation map and is converted + to RGB based on the FreeSurfer color look-up table. +- If the image has size 3 in any dimension, and the sum of squares + of the values in that dimension is approx. 1 everywhere, the image + is assumed to be a vector map and is converted to RGB based on + vector orientation. +- Otherwise, the image is assumed to be a generic intensity map and + is converted to grayscale. + +Return an array of RGB/Gray values the same size as the input vol. +""" +function vol_to_rgb(vol::Array) + + if isempty(color_lut.id) + error("FreeSurfer color look-up table is undefined") + end + + if ~any(isnothing.(indexin(unique(vol), color_lut.id))) + # Assume the input is a segmentation map, get RGB of labels from LUT + return color_lut.rgb[indexin(vol, color_lut.id)] + end + + dim3 = findall(size(vol) .== 3) + + for idim in dim3 + if all(isapprox.(sum(vol.^2; dims=idim), 1) .|| all(vol.==0, dims=idim)) + # Assume the input is a vector map, get RGB based on vector orientation + return RGB.(abs.(selectdim(vol,idim,1)), + abs.(selectdim(vol,idim,2)), + abs.(selectdim(vol,idim,3))) + end + end + + # Otherwise, assume the input is a generic intensity map + return Gray.(vol / maximum(vol)) +end + + +""" + show(mri::MRI, mrimod::Union{MRI, Nothing}=nothing) + +Show an `MRI` structure (an image slice and a summary of header info) in the +terminal window + +# Arguments: +- `mri::MRI`: the main image to display +- `mrimod:MRI`: an optional image to modulate the main image by (e.g., an FA + map to modulate a vector map) +""" +function show(mri::MRI, mrimod::Union{MRI, Nothing}=nothing) + + # Find non-empty slices in z dimension + iz = any(mri.vol .!= 0; dims=([1:2; 4:ndims(mri.vol)])) + iz = reshape(iz, length(iz)) + iz = findall(iz) + + # Keep only the middle non-empty slice in z dimension + iz = iz[Int(round(end/2))] + + # Find non-empty slices in x, y dimensions + ix = any(mri.vol[:,:,iz,:] .!= 0; dims=(2:ndims(mri.vol))) + ix = reshape(ix, length(ix)) + ix = findall(ix) + + iy = any(mri.vol[:,:,iz,:] .!= 0; dims=([1; 3:ndims(mri.vol)])) + iy = reshape(iy, length(iy)) + iy = findall(iy) + + # Keep only the area containing non-empty slices in x, y dimensions + ix = ix[1]:ix[end] + iy = iy[1]:iy[end] + + # Subsample to fit display in x, y dimensions + nsub = Int(ceil(length(iy) ./ displaysize()[2])) + + ix = ix[1:nsub:end] + iy = iy[1:nsub:end] + + # Convert image to RGB/Gray array + rgb = vol_to_rgb(mri.vol[ix, iy, iz, :]) + + # Keep first frame only + rgb = rgb[:, :, 1] + + # Optional intensity modulation + if ~isnothing(mrimod) + if size(mrimod.vol)[1:3] != size(mri.vol)[1:3] + error("Dimension mismatch between main image " * + string(size(mri.vol)[1:3]) * + " and modulation image " * + string(size(mrimod.vol)[1:3])) + end + + rgbmod = mrimod.vol[ix, iy, iz, 1] / maximum(mrimod.vol) + end + + if hasfield(eltype(rgb), :r) # RGB + # Add α channel to make zero voxels transparent + rgb = RGBA.(getfield.(rgb, :r), getfield.(rgb, :g), getfield.(rgb, :b), + Float64.(rgb .!= RGB(0,0,0))) + + # Modulate intensity with (optional) second image + if ~isnothing(mrimod) + rgb = RGBA.(getfield.(rgb, :r) .* rgbmod, + getfield.(rgb, :g) .* rgbmod, + getfield.(rgb, :b) .* rgbmod, + getfield.(rgb, :alpha)) + end + else # Gray + # Add α channel to make zero voxels transparent + rgb = GrayA.(getfield.(rgb, :val), Float64.(rgb .!= Gray(0))) + + # Modulate intensity with (optional) second image + if ~isnothing(mrimod) + rgb = GrayA.(getfield.(rgb, :val) .* rgbmod, + getfield.(rgb, :alpha)) + end + end + + # Display image + if mri.ispermuted + imshow(rgb) + else + imshow(permutedims(rgb, [2; 1])) + end + + # Print image info + println() + if ~isempty(mri.fspec) + println("Read from: " * mri.fspec) + end + println("Volume dimensions: " * string(collect(size(mri.vol)))) + println("Spatial resolution: " * string(Float64.(mri.volres))) + if ~isempty(mri.bval) + println("b-values: " * string(Float64.(unique(mri.bval)))) + end + println("Intensity range: " * string(Float64.([minimum(mri.vol), + maximum(mri.vol)]))) +end + + +""" + view(mri::MRI) + +View an `MRI` structure in a slice viewer +""" +function view(mri::MRI) + +#= +=# + +end +