Skip to content

Commit

Permalink
update release 4.1.0
Browse files Browse the repository at this point in the history
  • Loading branch information
@mfroeling
mfroeling committed Dec 5, 2024
1 parent 9f72129 commit 72f1dea
Show file tree
Hide file tree
Showing 5 changed files with 273 additions and 348 deletions.
4 changes: 3 additions & 1 deletion .vscode/settings.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
"Mathematica",
"MUMRI",
"Nifti",
"onnx",
"Paclet",
"QMRI",
"Relaxivity",
Expand All @@ -40,6 +41,7 @@
"transmural",
"Tversky",
"UNET",
"Westin"
"Westin",
"wlnet"
]
}
22 changes: 12 additions & 10 deletions QMRITools/Kernel/MuscleBidsTools.wl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2221,7 +2221,7 @@ MuscleBidsTractographyI[foli_, folo_, datType_, allType_, verCheck_, met_]:=Bloc
(* Perform tractography *)
{tracts, seeds} = FiberTractography[tens, vox, stop,
InterpolationOrder -> 0, StepSize -> step, Method -> "RK4", MaxSeedPoints -> seed,
FiberLengthRange -> len, FiberAngle -> ang, TractMonitor -> False,
FiberLengthRange -> len, FiberAngle -> ang, TractMonitor -> True,
TensorFlips -> flip, TensorPermutations -> per, Parallelization -> True
];

Expand Down Expand Up @@ -2559,9 +2559,10 @@ MuscleBidsAnalysisI[foli_, folo_, datDis_, verCheck_, imOut_] := Block[{
crp = FindCrop[ref, CropPadding -> 10];
refC = ApplyCrop[ref, crp];
size = Dimensions[refC] vox;

pos = GetSlicePositions[GaussianFilter[refC, 15], vox, MakeCheckPlot -> False,
DropSlices -> {1, 1, 1}, PeakNumber -> {0, 1, 2}];
pos[[1]] = Reverse[Range[0., 1., 1/(Round[Divide @@ size[[;; 2]]] + 1)][[2 ;; -2]] size[[1]]];
pos[[1]] = Reverse[Range[0., 1., 1/(Ceiling[Divide @@ size[[;; 2]]] + 1)][[2 ;; -2]] size[[1]]];

(*Function to extract slice data for 2D images*)
sliceData = Block[{slDat},
Expand All @@ -2584,14 +2585,15 @@ MuscleBidsAnalysisI[foli_, folo_, datDis_, verCheck_, imOut_] := Block[{
]&;

(*2D image function*)
make2DImage = ImagePad[ImageAssemble[ImageResize[#, {Automatic, 600}] & /@ Join[
{ImageAssemble[Transpose@{ImagePad[#, -3] & /@ #[[1]]}, Spacings -> 20, Background -> White]},
ImagePad[#, -3] & /@ #[[2]]], Spacings -> 20, Background -> White,
ImageSize -> {Automatic, 1200}, ImageResolution -> 300], 20, White
]&;
make2DImage = With[{di = Max[ImageDimensions[#][[2]]&/@#[[2]]]/4},
ImagePad[ImageAssemble[ImageResize[#, {Automatic, di}] & /@ Join[
{ImageAssemble[Transpose@{ImagePad[#, -5] & /@ #[[1]]}, Spacings -> 20, Background -> White]},
ImagePad[#, -5] & /@ #[[2]]], Spacings -> 20, Background -> White,
ImageResolution -> 300], 20, White
]]&;

addLabel = ImageAssemble[{
{ImageCrop[#2,ImageDimensions[#2] - {0, 20}, {Left, Bottom}]},
{ImageCrop[#2, ImageDimensions[#2] - {0, 20}, {Left, Bottom}]},
{LegendImage[#1, First[ImageDimensions@#2], #3]}}
]&;

Expand Down Expand Up @@ -2624,7 +2626,7 @@ MuscleBidsAnalysisI[foli_, folo_, datDis_, verCheck_, imOut_] := Block[{
(*import data and make image and export*)
{imDat, voxi} = ImportNii[imFile];
img = make2DImage@MakeSliceImages[sliceData@imDat, voxi,
ColorFunction -> cFun, PlotRange -> ran, ClippingStyle -> clip, ImageSize -> 600];
ColorFunction -> cFun, PlotRange -> ran, ClippingStyle -> clip, ImageSize -> 1200];
img = If[lab =!= None, addLabel[cFun, img, lab], img];
Export[fileNameO[partsO]<>".jpg", img, ImageResolution -> 300];
]
Expand Down Expand Up @@ -2653,7 +2655,7 @@ MuscleBidsAnalysisI[foli_, folo_, datDis_, verCheck_, imOut_] := Block[{
partsO["suf"] = If[hasKey, anaSeg[[2;;4]], anaSeg[[;;3]] ];
Export[fileNameO[partsO]<>".jpg",
make2DImage@MakeSliceImages[sliceData@ref, {sliceData@seg, GetSegmentationLabels[seg]}, vox,
ColorFunction -> {"BlackToWhite","RomaO"}, PlotRange -> Automatic, ClippingStyle -> Automatic, ImageSize -> 600]
ColorFunction -> {"BlackToWhite","RomaO"}, PlotRange -> Automatic, ClippingStyle -> Automatic, ImageSize -> 1200]
, ImageResolution -> 300];

(*make the 3D segmentation image*)
Expand Down
19 changes: 9 additions & 10 deletions QMRITools/Kernel/SegmentationTools.wl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -712,7 +712,7 @@ ConvBlock[block_, feat_?IntegerQ, {act_, dim_, dil_}] := ConvBlock[block, {feat,
ConvBlock[block_, {featOut_, featInt_}, {act_, dim_}] := ConvBlock[block, {featOut, featInt}, {act, dim, 1}]

ConvBlock[block_, {featOut_, featInt_}, {act_, dim_, dil_}] := Block[{
blockType, type, blockSet, repBlock, dep, rep, nam, dilf, sclf, cons
blockType, type, blockSet, repBlock, dep, rep, nam, scaleF, cons
},

(*short notation for naming layers*)
Expand Down Expand Up @@ -767,21 +767,20 @@ ConvBlock[block_, {featOut_, featInt_}, {act_, dim_, dil_}] := Block[{

"U2Net",
{dep, type} = If[IntegerQ[blockSet], {blockSet, True}, blockSet];
dilf = If[! #1, 2^(#2 - 1), #3] &;
sclf = If[! #1, 1, 2] &;
scaleF = If[! #1, 1, 2] &;
NetGraph[
Join[
Table[debugUnet[nam["U2enc_", i]]; nam["U2enc_", i] -> MakeNode[
(*scale up -> never, scale down*)
{1, If[i==dep, 1, sclf[type]]},
{1, If[i==dep, 1, scaleF[type]]},
(*skip in -> only input from above, skip out -> always for enc*)
{0, True},
(*config*)
{"Conv", featInt, {act, dim}}
], {i, 1, dep}],
Table[debugUnet[nam["U2dec_", i]]; nam["U2dec_", i] -> MakeNode[
(*scale up -> always, scale down -> never*)
{sclf[type], 1},
{scaleF[type], 1},
(*skip in -> accepts one skip, skip out -> never for dec*)
{1, False},
(*config*)
Expand Down Expand Up @@ -2053,9 +2052,9 @@ Options[GetTrainData] = {

SyntaxInformation[GetTrainData] = {"ArgumentsPattern" -> {_, _, _, _., OptionsPattern[]}};

GetTrainData[datas_, nBatch_, patch_, opts:OptionsPattern[]]:=GetTrainData[datas, nBatch, patch, False, opts]
GetTrainData[dataSets_, nBatch_, patch_, opts:OptionsPattern[]]:=GetTrainData[dataSets, nBatch, patch, False, opts]

GetTrainData[datas_, nBatch_, patch_, nClass_, OptionsPattern[]] := Block[{
GetTrainData[dataSets_, nBatch_, patch_, nClass_, OptionsPattern[]] := Block[{
itt, datO, segO, dat, seg, vox, augI, aug, nSet, pad
},

Expand All @@ -2068,13 +2067,13 @@ GetTrainData[datas_, nBatch_, patch_, nClass_, OptionsPattern[]] := Block[{
(*get the correct number of sets*)
itt = Ceiling[nBatch/nSet];
Do[
dat = RandomChoice[datas];
dat = RandomChoice[dataSets];
If[StringQ[dat],
(*data is wxf file format*)
{dat, seg, vox} = Import[dat];
,
If[Length[dat]===2,
(*datas is list of nii files {dat.nii, seg.nii}*)
(*dataSets is list of nii files {dat.nii, seg.nii}*)
{seg, vox} = ImportNii[dat[[2]]];
{dat, vox} = ImportNii[dat[[1]]];
,
Expand Down Expand Up @@ -2211,7 +2210,7 @@ PrepareTrainingData[{labFol_?StringQ, datFol_?StringQ}, outFol_?StringQ, Options
]
, {sf, segFiles}];

(*export the overview of what has happend*)
(*export the overview of what has happened*)
legend = Grid[{{}, Join[{""}, Item[Style[#[[1]], White, Bold], Background -> #[[2]]] & /@ {{"hole & n > 1", Red}, {"n > 1", Purple}, {"hole", Blue}}, {""}], {}}, Spacings -> {1, 0.5}];
head = Style[#, Bold] & /@ {"#", "Name", "Labels"};

Expand Down
99 changes: 36 additions & 63 deletions QMRITools/Kernel/TractographyTools.wl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -577,13 +577,14 @@ FiberTractography[tensor_, voxi_, {par_?ArrayQ, {min_?NumberQ, max_?NumberQ}}, o
FiberTractography[tensor_, vox:{_?NumberQ,_?NumberQ,_?NumberQ}, inp : {{_, {_, _}} ...}, OptionsPattern[]] := Block[{
minLength, maxLength, maxAng, maxSeed, flip, per, int, stopT, step, tractF, vecF, trFunc, ran,
tens, tensMask, inpTr, stop, coors, vecInt, stopInt, ones, dim, crp, t2,
seedN, seedI, seedT, seeds, t1, tracts, iii, drop, maxStep, len, sel, mon
seedN, seedI, seeds, t1, tracts, iii, drop, maxStep, len, sel, mon
},

(*get the options*)
{{minLength, maxLength}, maxAng, maxSeed, flip, per, int, stopT, step, mon} = OptionValue[{
FiberLengthRange, FiberAngle, MaxSeedPoints, TensorFlips, TensorPermutations,
InterpolationOrder, StopThreshold, StepSize, TractMonitor}];
SeedRandom[1234];

step = N@If[NumberQ[step],step, Min[0.5 vox]];
maxStep = Ceiling[(maxLength/step)];
Expand All @@ -602,66 +603,40 @@ FiberTractography[tensor_, vox:{_?NumberQ,_?NumberQ,_?NumberQ}, inp : {{_, {_, _
];

(*make the random seed points*)
SeedRandom[1234];
seedN = Total@Flatten@stop;
maxSeed = Round@Which[
NumberQ[maxSeed], maxSeed,
Head[maxSeed]===Scaled, First[maxSeed] seedN,
maxSeed === Automatic, 0.3 seedN,
maxSeed === All, seedN,
Head[maxSeed]===Scaled, First[maxSeed] seedN
];

(*random seed point selection until maxSeed is found*)
seeds = {};
seedI = SparseArray[stop]["NonzeroPositions"];
stopInt = MakeInt[stop, vox, int];
While[Length[seeds] < maxSeed,
seedT = Flatten[RandomSample[seedI, #] & /@ Block[{i = maxSeed},
First@Last@Reap[While[i > 0, Sow[If[i > seedN, seedN, i]];
i = i - seedN;]]], 1];
seedT = # vox & /@ (seedT + RandomReal[{-0.999, 0}, {maxSeed, 3}]);
seedT = Pick[seedT, stopInt @@@ seedT, 1.];
seeds = Join[seeds, seedT];
seeds = seeds[[;; Min[{Length@seeds, maxSeed}]]]
];

(*finalize seed points*)
seeds = ToPackedArray@N@seeds;
NumberQ[maxSeed], maxSeed,
maxSeed === All, seedN];
seeds = SparseArray[stop]["NonzeroPositions"];
seedN = Length@seeds;
(*for randomization where maxSeeds > seeds be sure to fill homogeniously*)
seeds = Flatten[RandomSample[seeds, #] & /@ Join[ConstantArray[seedN, Quotient[maxSeed, seedN]], {Mod[maxSeed, seedN]}], 1];
seeds = # vox & /@ ToPackedArray@N[seeds + RandomReal[{-0.999, 0}, {maxSeed, 3}]];
seedN = Length@seeds;

(*start the tractography*)
If[mon,
Echo["Starting with "<>ToString[seedN]<>" seed points and step size "<>ToString[step]<>" mm"];
];
If[mon, Echo["Starting with "<>ToString[seedN]<>" seed points and step size "<>ToString[step]<>" mm"]];

vecInt = MakeInt[tens, vox, int];
stopInt = MakeInt[stop, vox, int];
trFunc = TractFunc[#, step, {maxAng, maxStep, stopT}, {vecInt, stopInt, tractF, vecF}]&;

(*check if parallel or normal computing is needed*)
If[OptionValue[Parallelization],
(*parallel tractography preparation*)
If[!OptionValue[Parallelization],
(*normal tractography*)
{t1, tracts} = AbsoluteTiming@If[mon,
ind = 0; Monitor[(ind++; trFunc[#]) & /@ seeds, ProgressIndicator[ind, {0, seedN}]],
trFunc/@seeds];
, (*parallel tractography*)
If[mon, Echo["Starting parallel preparation"]];
ran = Thread[{vox/2, vox Dimensions[stop] - If[int>0, vox/2, 0]}];
t2 = First@AbsoluteTiming[
DistributeDefinitions[step, maxAng, maxStep, stopT, tractF, vecF, seeds, tens, stop, vox, int, ran,
TractFunc, TractFuncI, MakeInt, EigVec, AlignVec, VecAng, Euler, RK2, RK4];
ParallelEvaluate[
vecInt = MakeInt[tens, vox, int];
stopInt = MakeInt[stop, vox, int];
trFunc = TractFunc[#, step, {maxAng, maxStep, stopT}, {vecInt, stopInt, tractF, vecF}]&;
, DistributedContexts -> None]];
If[mon, Echo["Parallel preparation time: "<>ToString[Round[t2, .1]]<>" seconds"]];
(*actual tractography for parallel with memory clear*)
{t1, tracts} = AbsoluteTiming@ParallelMap[trFunc, seeds,
Method -> "EvaluationsPerKernel" -> 10, ProgressReporting -> mon];
DistributeDefinitions[tens, stop, vox, int, step, maxAng, maxStep, stopT, tractF, vecF, seeds,
trFunc, TractFunc, TractFuncI, MakeInt, EigVec, AlignVec, VecAng, Euler, RK2, RK4];
(*Redifine int function on parallel kernels, see help distributeDefinitions*)
ParallelEvaluate[vecInt = MakeInt[tens, vox, int]; stopInt = MakeInt[stop, vox, int];];
{t1, tracts} = AbsoluteTiming@ParallelMap[trFunc, seeds, ProgressReporting -> mon];
ParallelEvaluate[Clear[vecInt, stopInt, trFunc, tens, stop]];
ClearDistributedDefinitions[];
,
(*normal tractography with monitoring*)
vecInt = MakeInt[tens, vox, int];
stopInt = MakeInt[stop, vox, int];
trFunc = TractFunc[#, step, {maxAng, maxStep, stopT}, {vecInt, stopInt, tractF, vecF}]&;
{t1, tracts} = AbsoluteTiming@If[mon,
Monitor[Table[trFunc@seeds[[iii]], {iii, 1, seedN}], ProgressIndicator[iii, {0, seedN}]],
Table[trFunc@seeds[[iii]], {iii, 1, seedN}]
];
];

(*select only tracts within correct range and clip tracts that are longer*)
Expand All @@ -676,7 +651,7 @@ FiberTractography[tensor_, vox:{_?NumberQ,_?NumberQ,_?NumberQ}, inp : {{_, {_, _
];

(*output tracts move them to coordinates before cropping*)
MoveTracts[#, vox (crp[[;; ;; 2]]-1)]&/@{tracts, seeds}
MoveTracts[#, vox (crp[[;; ;; 2]] - 1)]& /@ {tracts, seeds}
]


Expand Down Expand Up @@ -723,17 +698,15 @@ TractFunc[loc0_?VectorQ, h_, stp_, fun_] := Block[{dir0, out},
TractFuncI[{loci_, stepi_, h_}, {maxAng_, maxStep_, stop_}, {vecInt_, stopInt_, tractF_, vecF_}] := Block[{
loc1 = loci, step0 = stepi, step1
},
(*perform the tractography*)
Flatten[Last[Reap[Do[
(*direction at current location*)
step1 = tractF[loc1, step0, h, vecInt, vecF];
step1 = tractF[loc1, step0, h, vecInt, vecF];
Flatten[Last@Reap[Do[
(*check for stop*)
If[stopInt @@ loc1 < stop || VecAng[step0, step1] > maxAng || Norm[step1] < 0.8 h, Break[]];
(*sow valid point and update for next location*)
Sow[loc1];
If[stopInt @@ loc1 < stop || VecAng[step0, step1] > maxAng || Norm[step1] < 0.8 h, Break[], Sow[loc1]];
(*Update location*)
loc1 = loc1 + step1;
step0 = step1;
, maxStep]]], 1]
step1 = tractF[loc1, step0, h, vecInt, vecF];
, maxStep]], 1]
];


Expand Down Expand Up @@ -790,7 +763,7 @@ EigVec = Compile[{{tens, _Real, 1}, {vdir, _Real, 1}}, Block[{
dxx, dyy, dzz, dxy, dxz, dyz, dxy2, dxz2, dyz2,
i1, i2, i3, i, v, s, v2, vv2, l1, a, b, c, norm, vec
},
(*tens = 1000 t;*)

(*method https://doi.org/10.1016/j.mri.2009.10.001*)
vec = If[Total[Abs[tens]]<10.^-15, {0., 0., 0.},

Expand Down Expand Up @@ -860,7 +833,7 @@ FindTensorPermutation[tensor_, vox:{_?NumberQ,_?NumberQ,_?NumberQ}, opts : Optio
FindTensorPermutation[tensor_, vox:{_?NumberQ,_?NumberQ,_?NumberQ}, {par_?ArrayQ, {min_?NumberQ, max_?NumberQ}}, opts : OptionsPattern[]] := FindTensorPermutation[tensor, vox, {{par, {min, max}}}, opts]

FindTensorPermutation[tens_, vox:{_?NumberQ,_?NumberQ,_?NumberQ}, stop_, opts : OptionsPattern[]] := Block[{
step, tracts, perms, flips, lengths, i, j, pl, ind
step, tracts, perms, flips, lengths, i, j, pl, ind
},

perms = {{"x", "y", "z"}, {"x", "z", "y"}, {"y", "x", "z"}, {"y","z", "x"}, {"z", "x", "y"}, {"z", "y", "x"}};
Expand Down Expand Up @@ -1237,7 +1210,7 @@ PlotSegmentedTracts[tracts_, segmentIn_, bones_, dim_, vox:{_?NumberQ,_?NumberQ,
colListT = colListC = Reverse[ColorData[colFunc] /@ Rescale[ran]][[RandomSample[ran]]];
];

(*reference environement*)
(*reference environment*)
If[mon, Echo["Making muscle iso volumes"]];
ref = PlotTracts[tractsF, vox, dim, MaxTracts -> 1, Method -> "line", TractColoring -> RGBColor[0, 0, 0, 0]];
ref[[1]] = {};
Expand Down
Loading

0 comments on commit 72f1dea

Please sign in to comment.