QMRITools: Processing Quantitative MRI Data: Live with the R&D Team

Hi everybody. My name is Martijf Guling and today I will be giving a talk about my toolkit, QMRI tools, which is intended to process and quantify quantitative MRI data. So, quick overview of the presentations, three parts. So why human right tools? Now it all started. I'll discuss a bit how I set up my packlet and how I use it in development and then of course how I use it in my own research. So why human eye tools? Well, my first experience with Mathematica was I think back in 2003 when I was doing a course at university, we're studying and that was taught in Mathematica. From that vision and multiscale image analysis, it was my first experience. And after that I continued my master, my PhD on the research topic of quantitative MRI of muscle and nerve. And before we're going to dive into all the Mathematica stuff, I quickly wanted to give you a quick overview about what we mean with quantitative MRI and especially of muscle and nerve. For that I'm going to quickly switch to a PowerPoint presentation. So in my research I do quantitative multi parametric muscle imaging. So for that I use an MRI scanner, which is basically a big magnet. And we place our subject into the MRI scanner and then we put something on top of the person. This is shown over here, which are actually radio frequency coils, antennas that we use to send out radio waves and we collect back and from that we generate MRI images within MRI machine spins that have a nuclei that have a spin will actually be parallel or anti parallel to magnetic fields. And in our body we have a lot of these nuclei and therefore we get a net magnetization moment. And that's actually what we're measuring in mri. And from that we can measure different contrasts and different properties of these magnetic spins. And in MRI we can actually do two kinds of imaging. So qualitative imaging, just taking a picture of the body. And then we can actually enhance certain structures by using different settings. But I'm focusing more on quantitative imaging. And there the MRI machine is more like a precision measurement device that you actually want to measure tissue properties with. And then all kinds of modeling and precision repeatability become more relevant. So my aim especially is to actually understand muscle microstructure. So we want to really understand something about what's happening in the human body on a cellular level with muscle and disease or injury. But for that we use MRI techniques, which is more macroscopic measurement technique in the order of millimeters. So we can use quantitative measurements that actually tell us something about the celerity of the tissue. And especially for my research topic, I really want to link quantitative MRI to muscle architecture and muscle function. So with MRI we measure a small volume of tissue, which is actually, in my case, muscle. And it can actually, when you're diseased, it can infiltrate with fat. And throughout my presentation I will be talking about three different measurement techniques. It's called water fat quantification using Dixon Imaging T2 mapping, which is a contrast of MRI and diffusion weighted imaging. And with all these different techniques, we actually assess different parts of the cell. I'm going to quickly explain these techniques. First. One is actually water fat imaging. Both water effect molecules have different resonance frequencies in the scanner and therefore we can actually tease them apart if we look at their relative frequencies. Because one will be fixed and the other will be oscillating. And when we apply the correct models. So this is the most simple model, we can actually try to estimate the water and fat fractions. But it's more complicated because we also have to take into account all kinds of imperfections of the MRI machine. But if we do that successfully, we can actually get maps, quantitative maps that actually says, well, at this location in the muscle, this is the percentage of fat present in this muscle. The other technique we're using is EPG based T2 mapping. T2 is actually relaxation time. It's a half life of the signal that we actually measure. Again, both water and fat have different half times of the signal. They coexist in the voxel. And we actually can measure that with special measurement techniques. So we can actually measure this decay curve, but it's always a mixture of the two. And because again, the machine is not perfect, it's not a perfect exponential decay curve which is more complicated. So therefore we use a different kind of model called EPG model where we can actually estimate the quantitative T2 relaxation of the water signal. And again that we can use to actually say something about on each location of the of the muscle, what is the quantitative T2 value of that specific muscle? And it can actually get be a marker for inflammation. So you see some curves that's a disease progression of different muscle diseases. The sham muscular dystrophy, backus muscular dystrophy and SMA in different muscle locations. And I can actually see the different disease types have different T2 time courses over disease progression. The last technique I will be showcasing is diffusion weighted imaging. So in every cell there's water and water moves around due to temperature. And in the cell membranes the diffusion is not isotropic in all directions. So Therefore we can get sort of a distribution of the likelihood water moves a certain part and then we can actually measure that with an MRI machine. And then by looking at how the diffusion of water actually happens, we can say something about the structure that water is actually L is it in an isotropic cell or is it in cylinders, or is it in sheet like structures? Again, for that we collect data typically in different directions and different diffusion weighing strengths. And again, we can apply a model to that to actually estimate the true diffusion constant and the directionality of diffusion in tissue. And again, from there we can make quantitative maps what's actually happening on a cellular level. But also the nice thing about this, because we can also get directionality, we can actually also make these fiber photography images like this, where we can actually study the macroscopic muscle architecture. So remember this image because I'm going to show you later on how I'm going to make this in Mathematica. So this was the PowerPoint part of the presentation. So back to Mathematica that you will know better. So how did it all start? So the code actually started when I was doing my Master's and my PhD back in 2008-2012 and I used a lot of the code I'm still using today to make my PhD thesis. That code actually has always stayed with me and I've been adding to that ever since I first posted it on GitHub in 2016. So if people are interested in the code, go to this GitHub link. It's been published in the Journal of Open Source Software. So if you want to see that publication, you can go here. And at the moment I'm also promoting the toolbox or sharing it with the community via two different websites.com and also a website related to my current research, the muscle atlas.org website. If you want to install it, you can directly install it from the GitHub releases page where you can actually use this packet. Install to actually install the packet on your own computer. So the packet also comes with some demo data. So let's start with loading the packlet. So currently we're at version 3.9 and it takes a while to load and I'll explain later why it takes so long to load. So it's a fully integrated packlet. Of course it has the shows up in the add on pages and also as fully documented. So we can open the documentation and we can actually see that for all the functions in there. Everything's documented and should work seamlessly with your Mathematica. So every Every function has help and documentation and options. So what is in the toolbox? So currently it has a lot of different sub packages, 25 to be exact, and exposes almost 500 functions to the front end. And these 500 functions of course has options and it's around 450. So which packages are in there? So this is the list of the the packages that actually come with the toolbox. So depending on my research topic, typically I group my functions into sort certain packages, for example different application when I did cardiac research like the diffusion tensor imaging, cryptography. So typically I try to group my functions based on topic. So if you want a full list of all functions, you can actually run this function that comes with the toolbox and it actually lists all functions alphabetically and also all options alphabetically. Or you can choose to actually have it listed per package where actually the functions per package are highlighted. So this is something I actually use to verify if all the functions show properly and if the functions are actually in the right package. So it's not only functions actually. The packlet also comes with some packed in external executables software that I use written by other people that I can call from my main toolbox so they don't have to leave Mathematica to run some of the external packages since this summer. It also comes with some pre trained neural nets that I will showing later some resources. So that's for example the demonstration data I will be using in this notebook, also a demonstration notebook that actually highlights all the functionality. So these are actually all linked in the extensions part of the packlet. So when I build my packlet everything is actually taken along. So you can see there's assets there like the demonstration notebook, some neural networks from language. NET format, but also a lot of like executables of external dependencies for different software operating systems. And I have some functions actually in my code. I can easily get the local location of these information. For example, on my computer I've installed everything in the repository and therefore these zip or executables can be found there. And it actually helps me to locate things on the computer when I need them in the functions. So locally I typically work with two versions. So the default Kymri tools that's installed in the packlet repository, that's the one that I actually make public. But typically I work in a development version and actually development version is just my local repository which you can load by using backlit directory load, backlet data rebuild. But doing that all the time, I got a bit bored by that so actually what I did is made a this file is placed in the base directory applications and then when I run this, this is actually reloading toolkit from the development directory. So I can actually switch back and forth between the two quite easy by just running one or the other. And to prevent shadowing, I'll show some tricks later. Most of my development and stuff I do in Visual City code at the moment with the Wolfram language extension. And I'm quite structured in my file layout. So let's look at one of these Wolfram language files that actually are in a toolbox. For example, the Elastics Tools toolbox. So I typically don't look at these files in Mathematica themselves because I mostly see them in Visual Studio code, but they're always starting with a title, the begin package. Then I group all the users meshes for the functions so they're always on top. Then I group all the options that come with the function, so I know which function, what's an option. And then I have a list of all the error messages for that specific file. So everything is contained in one file and it's always local. And then the private part begins. And because the declaration of the usage function is done here, all the definitions are in the private, so you can see that everything is nicely grouped. And I typically use headings to actually say what's in there. And I also quite consistent with trying to comment code because everything's online, I want other people to actually learn from it as well. So I'm typically trying to do quite elaborate code commenting. So all the files are the same layout. And then of course I add the package. So when I develop stuff, I typically do it just in Mathematica notebooks and play around with it until I'm happy with the functionality of the function and only then I move to Visual Studio code, copy the code there, make a good layout, add the documentation, specify the options and provide the error messages and then I post it to GitHub. So the main loader Qmarto loads all the other sub packages into the context path. So the list of sub packages is actually one of the variables within the toolbox. So these are currently the sub packages that expose to the user and that actually has this context. So it's all in the Kumari context. So after you have loaded the toolbox, actually there are a lot of different toolboxes are added to the context path. There are some tricks to like not have the context part be corrupted by everything, but for me this works best. So the loader itself is not just like get all the sub packages, but it also does some other smart things. So there's some like general parameters being defined, the sub packages are defined, the context of everything is defined. And I quietly already load the toolbox once. So I can get actually all the context names and all the function names that are located in the toolbox. So I do this quietly without the user seeing it. So I can actually destroy all system names that are already in the global context that are actually shadowing my toolbox. And therefore I can easily reload the toolbox or switch between two different toolbox because first I destroy everything and then I rebuild everything up again. So that's actually done here. In the end, that's the loader and I also have some commenting with that. So if I would load the toolbox with first specifying that this version is true, actually toolbox is loaded again, but actually it will echo back what it's doing. So it's loading version 3.90. There are these packages defined and these are the functions that are actually in those packages. Not really sure why it's not scrolling. There it is. So these are the lists of everything has been added. And then in the end you can actually see that it's removing all local and global definitions of every function that's in these lists. And then after it's done that for all the different sub tools, it will reload and protect all the definitions that are in there. So I can't override the definitions anymore. And this prevents from getting shadowing warnings. So all these different toolboxes are aware of each other. So every package begins with begin package the package name and then it has all this sub dependency. So basically the need call I always load the developer toolbox and then the entire context of the entire toolbox is locally loaded with the exception of the toolbox itself. And that makes it all the I don't really care about the hierarchy of when the in what order things are loaded, because every function in every context is available in all other contexts. So I also have some functions to actually visualize the function dependency. So for example, this is a function in my toolbox plot data and I can make function graphs of the plot data function that you can see here and which functions actually calls. So green is always a public set delayed function. So the double dot is sign the set elate and it calls actually all these blue functions and blue functions are local functions that are not exposed to user and then the yellow ones are actually functions that are public but outside its own context. So you can see that in here as a triangle, for example. This is a pure function. So color on is a pure function. So I can actually, by the shape I can see what kind of type of function it is, and by the color I can see where it lives. And of course I can do this for all functions. So that's nice thing about for example, manipulate, where I can actually say, well, let's look at some other functions, for example, and segmentation tools. Let's pick a function, segment data. So it's like these are quite large functions. So we can actually see that this is also a quite elaborate function. And then the red ones are actually private functions of other contexts and the stars are actually compiled functions. So we can look at some other things. For example, I know that in tensor tools, TensorCalc, this is one of the tools that does some heavy computing. And almost all the sub functions calls are compiled functions. So you can see all the stars here are all compiled functions. So compile is a function I use a lot. And when I actually do data analysis of one of my studies, When I have acquired the data, I have this one major overseeing script that actually calls everything and it actually a function is called Muscle Bits process and that actually looks in a directory, looks for a configuration file of that study, which actually specifies the folders, the options for the different processing methods it needs to use for a specific study. And once it knows the config file, it can actually run and process all data in all the subfolders in that folder. This is actually quite a big function which I think calls a lot of the functionality in my toolbox. So this is the network of that one. So you can see this is the main function. And in total there's like almost 900 both public and private functions. And this function actually calls almost a quarter of those. So when you build the toolbox, you also need to deploy it. I used to do a lot of manual steps in there, but over the years I've made a build script which can also be found on my GitHub that I use. And actually this has some like functions that are not in the toolbox but just specified in this local thing. And actually this helps me automatically deploy my toolbox so I can look what's currently installed on my system. I can actually check the documentation. Let's see if it actually works. So it's initialization cells and actually this will run through all the usage notes of my toolbox, finds what it actually is there. If these usage matches are correct and if some are missing and then the application build command actually allows me to build the application, including the extensions, rebuild the documentation if I want to. To generate automatically generate the packlet, actually push the release to GitHub so I can give a release title and a release note. So it actually will push the release to the release page of GitHub and will install a local version of the current release. So basically when I want to push a new release, I don't have to do anything. I can just run this function and everything happens automatically. And of course also for the HTML documentation, it can be generated with documentation tools. Nowadays in Mathematica, I can also build the HTML files. This function is taking a bit longer than I expected it to be, so probably we'll abort it. But this actually will give me an overview of all the documentation pages that are missing. So when I build everything, I build everything to my local workspace into a separate build folder. So it's actually this is the build folder. You can see the different packlets recently built with the usage notes that I use for GitHub releases, the full toolbox, and the HTML help. And HTML help is automatically then ported to my website. So actually when we go here, you can see that the same documentation that's available Mathematica, is also available on my website. So this took a while to build over the years, but I'm now quite happy that everything goes automatically. So I use it to process MRI data. MRI data is typically numeric data. And when I do a scan session with one of my volunteers or patients, I typically scan for around 45 minutes. And during that scan time I generate around 9,000 what math considers daikon images. But I would rather see it as 18 individual data sets, because a data set is a 3D representation of the reality and I take multiple versions of the data. So it's sort of like time dependent or different diffusion weighings. So it's 18 individual DICOM data sets that give a representation of the 3D anatomy and some other stuff, but that comes down to around 9,000 2D images. But you would never really look at them as images, only when you want to visualize it. And that's around 2 gigabytes of data. So for numeric computation, I always try to balance the memory dependency of all functions, the speed and the usability, so it can be really accurate. But if it takes a day to do one of these data sets, it's not really useful. And I also want everything to run on my local system, so I typically try to balance those things. So one of the things I did, I think for forever is set history length to zero, especially in the old days when my computers didn't have that much internal memory. And you're working with data sets of a few gigabytes, then this is typically a smart thing to do. So what does this data actually look like? Let's look some at some of the example data. I'm going to import one of the data sets. This is diffusion weighted data from the lower leg. And then of course I need to visualize the data. And for that we have to function plot data. It takes the data set itself and something called the voxel size. So MRI data is volume, so it's not pixels, but it's volumes. And these volumes can be anisotropic. So what do I mean with that? So this is the cross section of two legs. So we see this shin bone here, muscle data. We can scroll up and down. So we're going up to the knee, down to the ankle. And then because it's diffusion weighted data, it's also 42 different volumes with different diffusion weighing. So this is a data set that's actually 25 slices by 42 volumes. And in plane resolution is 160, 160. And we can also look at it from different angles. So this is the side angle, for example, the same 42 volume. So you can actually see really represents 3D anatomy. So and these voxels are actually anisotropic. So if you zoom in a bit, you can actually see that these things are not square anymore, but actually a little bit elongated. Another thing that I've written is something called memory usage. It's actually lists all the data sets above a specific size in my memory. And this is useful because in those cases I can actually sort of clear data sets if I don't want them anymore, or see what's actually eating up my memory and just clear those. So to keep everything efficient, most of the data sets I immediately convert to back numeric arrays. So here is a sort of a test case where I have generated a random real 100 by 20 by 100 by 100. And if you then remove the packed array, so you take away the packed array or actually even maybe replace one of the values, instead of being numeric valued for a integer value, you can see some performance differences. Oh, first let's generate the data. If I then would do some testing, for example, just transposing a data set or mapping over different levels of the data set, the performance difference between different types of the same data. So either being packed or non Packed numeric or not numeric can have a big impact in performance. For example, when I just transpose a packed array, it's faster than transposing a non packed array. So unless you do repeated timing, there's some caching going on I think because if you do repeated transposing that this is way faster. But one of the most important things for me actually why I would start doing this was if you just map on a very deep and big array and one of the values is actually an integer, the mapping operation becomes almost six seconds, while for a packed array of numeric values only half a second. And this can actually show up in your performance because I'm not doing anything to the data, just mapping over it. And that's the only thing I do. So a lot of things actually are packed arrays, numeric packed arrays or are the sparse arrays when I work with the data. So one of my favorite functions, block to keep everything private and manipulate of course to show stuff. So manipulate is typically when I need interactive user environments to optimize, for example sequences. So this is a tool that I use to optimize this water fat quantification where I can actually design the experiment by tweaking certain parameters. And I have more of those dynamic interfaces that actually help me optimize and design specific experiments with settings to see actually what's going on. So manipulate I use a lot and also the entire plot data. One thing I showed is actually one big manipulate window. Compile is something that if you work with numeric data, really speeds up the mathematical code if you know that you can actually compile the functions. So most of my compile calls actually always have the option parallelization, runtime option speed and runtime option listable. This listable is really valuable because I say well, I want to do the operation for example two vectors. And when you give the input it figures out how deep is your data and actually how should I list over it? And you don't really have to think about it. Same as mapping. But typically it also figures out automatically on which level it has to map to actually get to vectors. And this compilable functions function actually lists all the functions that Mathematica has in the compiler. Some code I found the Internet actually that list is. So if I write an algorithm and I see I'm using most only functions that are in this list, I know I can probably compile the function quite efficiently. And there are some like compiled values that actually doesn't really hurt of this in the compiled function, but it doesn't really speed up the performance. If you use compile a lot, it actually expects reals or integers. And then if you go to symbolic expressions, The symbolic expressions actually can hurt the compiler. It will typically given this thing like the compiler cannot do it anymore switching to normal code. But then the code becomes really slow again. So I do some tricks. So typically when a value in my data is 0, part 0, it's typically background and I don't care about it. I can never really measure a zero with Mr. Data. So I have some functions like divide, no zero log, no zero exponential zero mean of zero. It actually just ignores zeros. And if a zero is given then it returns a zero. So normally we'd get complex infinity or weird express or like symbolic expressions that the mathematical language really is good in handling. But the comparative function doesn't really like that. So there's a lot of functions that I use that actually say if zero in then zero out. It's of course mathematically not correct, but it really helps me. So other favorite functions of mine are, for example, that I misuse more than actually use a sparse array and unitize. So one of the examples I want to give here is that let's say I have an integer array and I want to replace 11 by 12 and the array is really sparse. I can use tricks like sparse array to actually speed this process up a lot, to only actually replace the explicit values that are not zero. And even faster is just with unitize, find the locations of things that are 11, make those 0, and then add 12 again on those locations. So again, let's do some tests. So this is a really sparse array with only a small block of values in the middle and the rest is all zeros. And this is a fully filled array from 100, 100 and doing the 11 by 12 replacements. And you can see that the third method is really fast. And for the sparse array, the sparse array trick is even way faster than a replace. But if it's not sparse, this actually slower than the original version. But this still is faster. So sometimes you have to switch between different data analysis methods. If you know that your data is sparse, you can switch to different tricks to actually speed up code. So why is this important for me? Well, if I look at my data, it's not really only the data. I also do data annotations. Like this muscle is the specific muscle. I label those with integer integers and I manipulate those a lot. So this is actually a proper data set. So that's actually not just a piece of it, but just A full leg. And on the right you see actually the labels that I use. So let's overlay those on the. No, that's the wrong one. Transparent overlay this on the data. So the right one is actually integers where we actually want to show which muscle is which. So we have integer labels that we can actually that we actually use to identify which muscle is which. And then if I want to manipulate these kind of data sets, they're actually quite big. For example, if I want to split those segmentations. So for this I use parse arrays because these data sets are quite big, 187 by 80 and it's actually would typically not really fit in your memory. But luckily these are really sparse. So low sparsity. And then I can see in the memory that the normal representation matrix reputation of the segmentations actually is 150 MBs in memory because it's not like split it up, but I split it up in all different labels and it's way bigger but actually lower in memory because it uses a sparse array head. So this really helps reducing the amount of data I need memory. I need some other custom functions is import nifty so that the data that I use as store as nifty format which is actually supported by Mathematica natively with the bioimaging format. So the import export registers are done locally in the toolbox so they don't conflict with the normal loading of Mathematica of nifty files. So this is one of the data sets in the demo data folder. So in the demo data folder I have a list of nifty data sets that actually used for this demonstration. And I can load it with the normal import function. And I can also get some meta information from these files, but I can also load them with my actual own importer. So we're going to look at the. So here actually the data is imported as images. So the data is imported as images, but they're. The image representation is not correctly if the rescale data to get to the actual data and it gives some meta information. But that's not the full story about it because here I actually import the same data set with a binary loader that actually gives back not images but just data. So you can see this takes around 12 seconds. This is only one second gives me back the dimensions of the voxels but also gives me the full header which actually is a lot of things saved in a nifty file that the default loader doesn't do. Another trick that I use a lot is we have the function automatically says, well, rotate left or rotate right. So I actually wrote a function that says rotate dimensions right and left. So basically I have a Data set that's 25 or 42. I can say rotate dimensions left twice. And then actually you can see that actually dimensions are rotating. Sometimes you want the looping index to be the first one, sometimes you want to be the last one. And actually it's just under the hood, a transpose, of course. But this is more convenient then masking data. So like I said, backgrounds I don't really want to process. So if there's data there, I typically want to ignore it. And then we can of course do with binarization. So masking is basically just binarization the data. So I'm not binarizing everything above 10. And then if I look at the data that I just made, and then the mask is actually nicely overlaying on the data, but there's some holes in there and it doesn't look really pretty. So this of course has a lot of extra options with it where I can say, well, I want to smooth my mask, I want to close holes, I want to dilate a bit. And then the mask looks way smoother. And then after I've done that, so I applied a mask to the data, I can actually crop out the background. So auto crop data will actually find how to crop properly crop any dimensional data set, remove the background, so it actually is smaller in memory because just remove zeros and it will process faster and will also go to functions faster. This of course can also be done in a visual way. For example, say, well, I just want to manually crop one leg and then I can do the same thing and then it will actually crop a different size. And of course we also need to be able to reverse these crops. So I showed you the plot data. There's also a 3D version of plot data in development. It actually shows all three different views at the same time and gives you interaction options where you can actually also scroll to the 3D slices like this. So you can actually also view the data in actually the correct way. And of course the image 3D also allows to render the volume, but now it becomes a bit sluggishly. So there's still something you have to optimize, especially once you want to scroll through the volumes while the image 3D is there. It's not really performing that well yet. But I think visualizing data is really important because you really need to see what kind of data you're working with. So there's a lot of different ways to actually visualize the data. In my toolkit you can change colors. You can actually do a lot of cool stuff with that. So in the beginning I already said, well, I'm always balancing between speed and memory. So one of the things I did for that was plot data. To actually make this work, I wrote my own functions to actually generate really fast color settings. So actually what I do is not really calculate the color when needed, but I pre calculate like sort of an integer array of colors and then I convert the data to integers. So I just have a lookup table of RGB colors with opacity. So I can actually generate images way quicker than with the normal image hat. So these are two different ways you can do it. Like you can do image and then colorize or you can make an image of a numeric array if those are color values. So in getting the color values by either using a color function or this can actually have a. So both methods are really fast. But if you really want an interactive environment, I don't really care that I have the exact color. I want to be close enough. And you can see this is the custom color map and this is the normal image. This is already really fast, but this one is actually 34 times faster than the one. And for small images it doesn't really matter. But if you go to large images this starts to take quite a long time. And in an interactive environment you don't really want this. Another location where I really have done some optimization for speedup as for example, when we do the fusion tensor imaging, you actually have to calculate eigensystems from tensors. So here I'm going to simulate 10,000 tensors, random ones, with some custom functions. And from each tensor I actually want to calculate. So this is the tensor, they're all symmetric and positively positive. Definite so means positive eigenvalues. And I'm actually only interested in the first eigenvector and the eigensystem is really fast, but I'll show later why I need an even faster version. So I actually do a direct method of calculating the eigenvectors with the compile which is. So I'm now doing the timing experiment. So the normal function does 10,000 in 0.1 second, which is really fast. But this one is actually faster and I do lose some accuracy. But for my application it really doesn't matter. So in the last part, last few minutes, I want to give you some demonstrations of toolbox functionality. So for that I'm going to use the Demo data. And this demo data is also highlighted in this demonstration notebook that's packed with a toolbox. So this is a notebook that demonstrates most functionality within the toolbox. So I, I put it with the tool so people can actually play around with it. But also this is a different function for me that each time mathematics upgraded, I or I upgrade my own toolbox. I run this entire notebook to see if everything's still working as I expected. Sort of, sort of a test notebook for me as well. And it contains data, various data sets and labels. And by now we should have like getting the memory a little bit more full filled up. So I'm going to quickly clear some of the bigger ones because some of the data sets I'm loading in a bit need more memory. So I already talked a bit about data visualization. So the interactive things are nice to explore data, but of course sometimes you also want to automate things. So we have things forms that can actually automatically find interesting slices in data and then make images out of it with different color bars and imaging legends or do the same when you have segmentations, overlay those on images. So this actually when I'm processing data, I'm constantly making automated images of the data to actually track the processing pipeline. And of course when you go to publications or presentations, you also want fancy 3D visualizations. So for example, for the segmentations, we can actually make nice renders of our 3D segmentations in the toolbox as well. So now I'm going to show you some real data from the Dixon, which actually is so Mr. Data. Up to now I've shown grayscale images and that's actually, that's the magnitude part of the Mr. Signal because actually when we measure EMR signal, it's a complex signal, so it's consists of a real and imaginary part. So actually the data that comes from the scanner is not looking like the way I typically show up to now, but actually has a real component and imaginary component that looks more like this. And from each complex number you can also represent it of course the magnitude and the phase of the signal. And actually the phase is really interesting properties as well. And so here you can see this, that water and fat have different face properties. And this phase data has to be analyzed as well. So one of the system functions in my toolbox is phase unwrapping, which is really important for Mr. Data. So I'm going to simulate some phase images that have actually have reps in there and some noise and also some hard splits. And there's different ways of unwrapping these data sets, either with best path methods or with direct cosine. And again here, this one is a bit slower, but it's more accurate. So you can actually see here the rep data set, the true unwrap data set and the data after unwrapping. And you see the difference only in the noise, unless you have splits, because it's pretty difficult to estimate what's the offset. So this is quite slowish, this one is way faster, but actually it's less accurate, so sometimes you want to switch out accuracy. So let's look at some data. So here I'm actually going to process one of these data sets. So here I'm going to load in the data set, pre calculate some things I need, which amongst others is this phase map that I will unwrap with the unwrapping algorithm and then I will actually perform the, the reconstruction of this equation that I showed in the presentation earlier on this data set. So I'm doing two unwrapping sessions, one for each leg. And then here we can actually see how that looks. I should start this in the meantime. So this is before unwrapping and this is after unwrapping. And this data doesn't contain many reps, only in the edges you can see some errors and both these things are queued so that block window won't respond to any updates until this computation is done. So this will actually take this, this complex data, some of the pre calculations and actually fit all the things in this equation and give those as an output. So you can see the FAT fractions, the water fat maps, the phases that are here. And we can actually look at what then comes out of this algorithm. So you can see this, this was done in, in 32 seconds and that nonlinear system, so I, it's quite optimized. So here you see for example the phase map comes out of it and the T2 star map as we call it. And of course in the end we really interested in the FAT fraction maps. So you see the data set again and on the right you can see the water and FAT fractions of this data. So this is also the data set we use for, for, for annotation which muscle is which. And recently we've introduced a neural net for that. So neural net training is fine, but actually to make it usable, more is needed. So actually four neural nets, one that actually says, well, is it a left or right leg or is it both? And according to the classification of the data set, it will actually go to a different pipeline. Then we have to know is it upper leg or lower leg? Because we have to segment different muscles. And once we know that, we can actually have to patch the data into data into patches that actually can be done in the memory of the computer. And then we segment both the different legs and muscles and then we put everything together again. So for this we use pre trained models that come with the toolbox. For example, this unit that actually will do the thigh muscle segmentation. So this is a convolutional neural network that takes patches of data of a certain size and can actually generate the label outputs. And it's actually done with this function, segment data. And that segment data will actually go through this entire pipeline and actually put the logging of this function on so you can actually see what it's actually doing. So it's analyzing the data, it's segmenting the legs, so it's found a lower right leg and a lower left leg and it cropped it to the minimum needed size of that leg so it can actually have less data being ported to the network. And then it's now segmenting this right leg first and now it's segmenting the left leg and then it will put everything together again. And now we actually, this segment that we're going to show, I'm going to show the 3D version actually. So these segmentations have now been done by the neural network real time. Normally this would take me half a day to do. And here actually we have these segmentations done by the neural net in a few seconds. So that's actually saving me a lot of time at the moment. For the sake. So for the, the second part of the method, it's a 2, 2 mapping. We actually use the block equations that actually define the Mr. System, which we of course can simulate using Mathematica. And we do that because when we do a measurement, as I said before, we actually send out a radio frequency pulse to the system to actually excite a certain type amount of tissue should typically call a slice or volume. And this will be the envelope of the, the RF waveform and it will actually go inside with a certain slice profile, which is not perfect. Nothing is perfect in the real world. So we need to, to actually understand what's happening. And we're not giving one RF pulse, we're giving a series RF pulses of different frequencies and different shapes. And actually with this toolbox we can actually simulate what's going on. And then for the fitting, we actually need to know exactly how our signal behaves over time. So this easily, equations I showed before, this is actually the signal over time which we can actually generate with this function. So these images that are being made here, this function actually generates them. So you can see here and then if you want to generate one of these signals, we use this function. So as you can see, EPG signal is the same as we saw in the equation here. And although this is quite fast, if we use this in a minimization algorithm, it's too slow. So there's one thing I can actually use as a trick. So we can pre calculate a dictionary of all possible solutions in our solution space, which I can do with this function. So we're creating a dictionary of all possible solutions that we know that can happen within the search base. And we can do this really quickly with using parallelization and in this case with parallel table. So once we've done it once, we can actually use these solutions in our minimization algorithm rather than the full model. So of course we use some, we lose some accuracy, but we gain a lot of speed. And actually this, this calculation of this dictionary takes quite a lot of time. So I don't really want to do it each time I go to a data set with the same settings. So within the toolbox I use some tricks to actually cache these kind of large computations. So hopefully the prediction of this thing is accurate and we're almost done here. So now I've created this dictionary with these settings once and now if I would rerun the function again, it's instant, so it's already knows I already calculated this, so it's already in my memory. So how can you do this? Easily, in Mathematica, if you define a function with set delayed, you give the same function again, but then not without set delayed. So if you then define this function and you run it once with certain settings and you look at the definition, these things are actually stored in memory and you can just neatly access those. So I want to end with some examples of diffusion tensor imaging. So again, I'm going to start with importing some data sets. And as I mentioned before, with diffusion tensor imaging, we actually measure the diffusion in different spatial orientations, but also with different diffusion directions, because we need to solve this equation. So for this specific experiment, we did measurements with these directions spatially orientated on different diffusion strengths, different shells we call it. And to actually design these experiments, we have this interactive window where we can actually say, well, I want five shells with different B values and we can actually generate them evenly spaced. So we have the most Perfect measurement experiments. And we can of course visualize it in different ways. So this is another example of how we can use manipulate to actually help us design experiments. So for the fusion data we do a couple of processing steps, for example, data denoising. So when we have our data, I'm first gonna do some preparation and then one of the first things is that we want to remove noise from our data because it really helps our fitting algorithms. So so for that we use principal component analysis denoising. So I'm going to denoise this data set. Now actually how this works is that we have our Signal. It's a 3D space with a fourth dimension. We take a small patch out of it, we do a singular value decomposition, look at the eigenvalues, and only going to take along the eigenvalues that contain enough signal and the rest we throw away basically. And then we reconstruct the signal. So here I'm running the denoising on this data. So again, this is one of the functions that we optimized for speed and not accuracy. And this is taking longer probably because I'm streaming than normal. Once it's done, we're going to compare the two data sets. So it's really important, we don't really lose any of the information in the data and we only actually remove noise from our data set. So again, looking at data is really important. On the left side we have our non denoise data set and on the right we have our denoise data set. And if we go to different volumes, especially with the higher diffusion weighing, we can see some noise here and that's removed here. So to actually evaluate, we can look at the difference of these two data sets. So this is the difference image of these two data sets. And you can actually see what is removed in all different slices is actually noise and not structured. So how can we validate that? So let's subtract the two data sets and look at the distribution of what's removed between the two data sets. And here we can actually see that the thing we remove from a data set is a normally distributed signal, which means it's noise. The last thing we do to our data before we can actually process it is remove motion from it. So here I'm going to use this example data set and actually purposely rotate and translate it. And then we're going to compare the two data sets. You can actually see left is the original or right is moved. So of course we can also quickly look at how much is moved along and within the mathematical language you have the option to do, I think it's called Image Align. But for registration I use an external toolbox called Elastics, which gives me larger flexibility. So I have my original data set as a target and my warp data set. So now I've registered it and hopefully now the two are. The register set is aligned with the original data set perfectly. So this is actually an executable that I call externally. So I export a lot of the data to a temp folder and there the executable runs on it and then I import back the output of that executable into the mathematical language. And the nice thing about this, you can also get back the trans. So this was the imposed translation rotation. This is the actual translation rotation found by algorithm. So once we have everything fully processed, we can actually do cool stuff. So in the end we're going to use fiber triptography, as I promised on our data set. So our data actually after processing and fitting is actually represented as a tensor. So tensor is a three by three, in this case a three by three matrix, which a diagonal. On the diagonal we have the, the diffusion and on the cross diagonal we have some cross dependencies. And on that kind of data we can actually do tractography. And here we need this eigenvalue algorithm that I showed earlier, because within this data set I'm going to start at a specific point, calculate the first eigenvector of the tensor, follow the direction to the next step, get new location and and so on and so on. So I'm going to define some stock criteria and then actually do the fiber tractography. So this case I'm using an Euler method to do the stepping and zero order interpolation for speed. The fiber tracks I need to find are within this range and I can't make huge bands, so I can't make bands bigger than 25 degrees. So it's seeding 10,000 fiber tracks and there's a step size of 1.5 millimeter. This took 17 seconds and it generated around 500 tracks per second. And in the end it found 10,000 tracks. So how does this look like? Of course we have a plotting function for that. So I'm using the masks for segmentation the tracks I found and I'm going to label everything accordingly. So I'm plotting only 5,000 of the 10,000 tracks. So here we can see our scene where we have actually from this data set, not only the segmentation but also the muscle fibers. And you can see the back of the calf where the muscle fibers actually Go outward inward. And this I'm using for quantification of muscle force. So this is the toy example I'm going to stop with showing last thing. I want to show some really the data that actually I'm actually working with, and not the example data, but the full thing. So here I'm going to import all volume tracks of the entire leg and also the labels of the entire lag. And hopefully I can show everything. So why I said speed was important. So in the previous example, I generated around 10,000 fiber tracks. And after this is loaded, you can see when the data sets become bigger, the amount of tracks also become bigger. So in this specific data set, I have a quarter of a million fiber tracks generated with an average length of 180 steps. So actually I did the eigenvalue calculation within this tracking algorithm over 180 million times. So then this factor 10 speed up is actually quite useful. And you can also see that with actual data, the memory usage is quite more than the toy examples. So again, of course we need to evaluate the tracks so we can actually get some parameters out of these track maps. So for example, the seed density and the track density. So the seed density is where did I start? Each fiber track. And then the track density itself is how many tracks are at every different location. So you can see it's nicely evenly distributed. Make this a bit bigger on the right, let's change to a different color map. We can actually see that the fiber track density actually coincides with the muscle. But the more interesting things are actually the average track length. So how long are muscle fibers and what's their local angle compared to the human body? Because it actually is a determinant force. So again, the same maps on the left you actually have a track length map. So you can see that different muscles actually have different fiber lengths. So you can see, for example, the one here on the inside has quite long values. Actually you can put on the legend. So this is in millimeters. In this specific case, you can actually see that some muscle fibers are quite long, some are quite short. And on the right we can actually see the angle of our tracts. So you can see also that the angle of the different tracks actually is different for different muscles as well. And these are actually the values I'm interested in. Of course, we can also visualize all this. And I'm gonna just load the contours I need for this and the threats I'm gonna load. So these data sets typically are around 1 meter in real life size. So the, the length of Leg, of course we're on 1 meter and the resolution we're actually acquiring a data at is in the order of 6 by 1.5 by 1.5 millimeters. So we have around typically between 1 million 1.5 million data points in our lag where we actually sample everything. So in the end when we actually want to bring everything together. So I've loaded the tracts, the muscles, the bone segmentation, I know the dimensions of the data, the voxel size of the data. I'm going to try to plot the image that I showed earlier in the presentation. And this is always a bit scary. At some point my front end will probably hang. Everything's pre calculated. I'm just loading the graphics complex. I hope it will render. And with this. Yeah, there we go. Yeah, here we are. So here we have this image that I showed earlier, but now as a nice graphics in, in Mathematica, my chat window that I was trying to see the questions in froze a bit during the presentation. So I think I missed a few questions along the way. I think for the gradient optimization it's electrostatic repulsion. The spatial resolution indeed of the data set I sort of showed. So this Data set, it's 164 by 320 by 320 with a spatial resolution of 6 by 1.5 by 1.5 millimeters. And I see also some questions about other nuclei than, than, than protons. So of course MRI can also work on, on different nuclei for which we also have tools in place within the toolbox. So with this image I actually want to end my presentation if I don't see any other questions in the, in the chat. So I really want to thank you for, for tuning in and listening to what I had to tell. Hopefully you got some tricks out of it. And if you want to see more of the code, feel free to visit the GitHub and give me maybe some tips on where I can still improve on my code because I'm always willing to learn. So thank you all and with this nice image, I want to end the session.