Python Image Segmentation Tutorial (2022)

so the other day i was reading a paper with research that i was doing which involved segmenting the lungs from a ct scan and i found this github profile here from a guy named sergey primakov so shout out to whoever you are from maastricht university and he had written a paper in nature communications by the way nature communications is a very highly repeatable journal and i noticed in the article that he actually provided the code from his project where he segmented the lungs publicly and it's actually right here and you can take a look at it and this is state-of-the-art stuff and he's using all the modern python libraries to take basically a 3d ct scan and just return a mask of the lungs and so that's what we're going to look at in this video if you enjoy this video please remember to subscribe as always join the discord server as well and finally i just started an instagram account for this channel as well so check that out anyways enjoy [Music] so the packages today are kind of important we have numpy and matplotlib as normal but the main image processing stuff that is state of the art in python that we'll be using is uh scikit image is a beautiful image and processing package it's got a little bit of machine learning stuff scikit learn is the machine learning library associated with it and then also scipy we have some stuff from the nd image library which is going to allow us to do make measurements of the image compute center of mass of various regions and some other functions that we'll get into in this video so the first thing is we're going to do is we're going to open our data i've put this data in the git repository so you should be able to access it it's just one directory back in the data repository it's a ct scan and we'll look at it and so here i just load it using numpy and the first thing to do and this is something i do always when i'm doing machine learning with images or any sort of array that i open the first step you should always get used to is looking at the shape of the image because this will tell you about what you just opened so i have no idea what to expect here well actually i do but you know you might not have any idea what to expect and you see the shape is 263 by 512 by 512. and the shape is the shape of the image and so what you have here is not a 2d image but a three-dimensional image and so what it is is the 263 is along the height of the person so from 0 to 263 and then the 512 by 512 is like you're taking slices of a person sort of like this it's called an axial slice and so you have a 5 512 by 512 image and you have 263 of them and they're sort of like a sliced ham of a person so we'll look at that and so what it means is we're looking at a 3d image and what we want to do is we want to process all of these slices and i'll let you know what our goal is pretty soon so let's look at a slice of the image here i'm going to take the 170th slice of the ct image and you can see that our image is here and what you're looking at here is well at the bottom you have the table and it looks like this and so the table is sort of whatever 500 is and then it's sort of an air gap in the center of the table and then you have like a slice of a person sort of directly like this it's like you cut someone up like a ham and you're looking in here and what you have is the two lungs and then you have all the tissue on the outside and you have this sort of arbitrary scale and you might not know quite what it is and we'll get into that so the first thing to do we're looking at a new type of data we should probably understand what the scale means and so the units correspond to houndsfield units and a houndsfield unit at every point x y so each pixel value has a houndsfield unit quality so i give an xy you'll tell me how many houndsfield units it is and that's the value given by the color bar here it's equal to a thousand times whatever the attenuation coefficient is at that location in space minus the attenuation coefficient of water divided by the difference between the attenuation coefficient of water and air and so you might think well why are you using this scale that seems a little bit bizarre and there's also the linear attenuation coefficient which is defined such that if i have a beam of intensity i not coming into a material of length x the intensity coming out is going to be i naught e to the negative mu x so it's how many photons get attenuated in the matter so the higher attenuating this is the higher the probability that an x-ray shining through your body gets absorbed by your body so these higher regions here correspond to higher probability of attenuation and it's sort of scaled like this and so the reason they define it like this is it gives a really nice scale for differentiating between different types of tissue so here i have for example air is minus 1 000 that's just something you remember a fat is minus 120 to minus 90 and then you got soft tissue 100 to 300 you have all these different values for bone blood um and sort of all this stuff here and you can determine exactly how many hounsfield units something is so looking at this image you can look at a particular color and it tells you something about what material the actual point in the image is made up of so maybe there's bone uh here in the lungs you can see it's close to minus 1000 because the lungs are mostly air the same thing on the outside and the table is pretty attenuating so it must be made of some sort of metal and then there's sort of an air gap in the center and so we're going to do image processing on this image and our goal is going to be to look through all the slices so here i'm looking at the 170th slice there's also the 160th slice and i can sort of go through all those slices from all the way from slice 0 which this is the legs of the person you can see the thigh bones here on the table all the way to the last slice of the image which is probably going to be the head and the shoulders like this and you can actually see the skull like this this higher attenuating region here and then the shoulders and the shoulder bones what we want to do is we want to get a mask of the lungs in the person we want to take this image and we want to find a corresponding 3d array where 1 means you're inside the lungs and 0 means you're not this is getting a lung mask from an image and of course it's not going to be easy to do because you can see that if i wanted to segment the lungs like this that's what we're doing it's segmentation here there's all this stuff in the center of the lungs too and then maybe there's air out here and how do we differentiate between air and all the other regions and so i'm going to show you how to use all these libraries to do a fairly complicated task the first thing we'll do is we'll create a threshold mask and so we're going to use -320 hounsfield units as a lower limit so this just gives us an array of boolean values i look for all the pixels where the image is less than minus 320 houndsfield units which if i look up here minus 320 is around here so that'll sort of select out the areas with air that's sort of like a preliminary step to segmenting the lungs and we'll go through a bunch of steps like this so this just gives me an array of trues and falses and i can set that equal to something called mask and of course then i can plot a slice of this mask at 170. so you see that all the errors in yellow here like this it's true when you're less than -320 hounsfield units which is getting close to air so you can see the lungs are segmented the background is all air here that's segmented and i can go through various slices i can look at mask 160. it's sort of the same thing here i can go to mask 10 it's segmenting all the air and all these various slices so that's our first step that we're going to do to segment the lungs of course here you have some air pockets and we'll deal with those so here's the lungs here again and now we're going to use the clear border function this is one of the functions that we actually imported up here from scikit image segmentation and what that's going to do is it's going to clear this outer border of in this case true values it's a boolean array so what we want to do is we want to apply this for all the slices in the image and so we're going to use np.vectorize here and this will come up a few times in the video what we're doing is if you look at the shape of mask mask dot shape we have a 263 by 5 12 by 5 12 image we want to do this segmentation we want to remove the border for each of these 2d slices so we want to loop through all 263 slices and apply clear border to each one of them so in order to do that i use numpy.vectorize which automatically creates a loop over whatever i feed into here the issue is that mask is a 3d array and we want to loop over 2d slices so what we tell the function in vectorize is we say look the thing that i want you to loop over is going to be a 2d array and it will know it'll say okay now i know i'm only looping over the first axis and i'm going to apply these to 512 by 512 slices that's sort of how you give a signature in np dot vectorize i'll make a video on this in the future because mp.vectorize is a very powerful function but if you have a 3d array and you give mp.vectorize it's going to try to call this on every element of the array you don't want that you want it on every 5 12 by 5 12 slice of the array and for that that's why you have to give this signature of saying we're giving it a 2d array and it's returning a 2d array in this clear border function so i'm vectorizing clear border which takes in a 2d array and it returns a 2d array and i feed in this 3d array mask it's going to do it for every slice in the 3d image which is the mask and it will return a corresponding 3d array where i've applied this function to each slice of my array and so i can do this and then i can look at the mask the same slice here and you can see that it's taken this outer border and it's removed it and that's exactly what we have here so now we're gonna do a little trick we have the lungs but you can see that the table also showed up that's because the table kind of looks like a lung and you might think well no it doesn't but to the computer it's it's a basically a hollowed out region of air surrounded by high attenuation coefficients so the table's a real pain actually for our challenge the reason why it's a pain is in some slices if i compare the area of the lung to the table the lung is going to have a larger area than the table but if i look at a different slice for example here so here the lung might have a larger area than the table but in a different slice you know for example this is part of the lung but it's clearly smaller than the table so what we want to do is we want to create different regions for each of the separate volumes in this image that's going to be the first step to just selecting the lungs so for example this might be region one region two region three and for this we have the label function which we imported from um scikit image dot measure here it's gonna take a boolean array where you either have ones or zeros like this and it's gonna take each separate distinct region and give it its own integer value so what i'm doing is i'm calling it on the 3d array mask remember that mask is a 3d array still and i'm saying okay you're going to apply this function label it works on 2d slices it's going to work on each 5 12 by 5 12 slice it's going to loop over all these slices and return a corresponding 3d array which is mask labeled and so i do this and you can see that i have different regions here and i'll give a little color bar just to show what's going on here so you can see that like for example this has the integer value seven this has the integer value six there's a few other little areas that you can see here that have also been segmented and the tables like this so i have all my regions in the image given separate integers and that allows me to differentiate between okay this is long one long two the table i now have integers corresponding to the different um segments of the image now the reason i'm doing this is because i want to keep the largest three areas in every slice so i can delete these little areas here and these little areas here and i'm going to keep lung lung table and you might think well why would i do that why would i keep the largest three areas well remember i'm doing this for all slices and like i said before in some slices one of the lungs is bigger than the table for example in this slice but in another slice maybe down here you're at the bottom of the lung the lung is going to be smaller than the table so we'll deal with the table later but the purpose of taking the largest three areas is i just want to get rid of any little areas the table will be left but it will guarantee that i keep the two lungs and then also unfortunately the table at this point so how are we going to do that well the first thing we need to do is we need to compute the areas of each of these separate regions and so the way you do that and i sort of go over that here so here i'm taking a slice that's basically the slice that i have plotted here so i call that slice slice is a 2d array and then i can use this function called region props which i imported above region props is from scikit image.measure library so what region props will do is it will give the properties of each of the regions in this slice so i have seven regions of this slice and region props is going to give me the properties of these so i can just going through this slowly here's my slice my slice is a 2d array my 2d array is separated into different regions here and if i call region props on this slice i have you know a list and then there's different properties for each of these regions and if i want the area of each region for example i can call for example here's my first region prop that's corresponds to the region for the integer 0 for example and i can call dot area it will give me the area in pixels of this particular region of this slice so i'll iterate through all the region props remember i have an array like this so i use list comprehension i take r dot area for each of the r in the region props and that gives me my areas so here i can set region props is the region props with the slice i get my areas i can look like this and i have my areas of all the different slices and pixels i have this sort of larger area and then these two sort of smaller areas here and then what i'm going to do remember this is the areas i'm going to take the areas from smallest to biggest so index three is two this is the smallest area remember zero one two three this gives me my smallest index uh four which corresponds to this here 0 1 2 3 4. 5 is the second smallest and then 34 is the next then 35 and you can see that the indices are in order and then what i'm going to do is i want from largest to smallest so i'm going to take this backwards that's why i give the two colons and the negative 1. so 0 is the largest area then 6 900 is the second biggest area so these are my list of my integers in order of biggest to smallest area so i have these indices and now i'm only going to consider the three largest areas by looping through the indices so the first thing i do is i create a new slice and i make it all zeros and i make it the same shape as slice so my new slice for example here this is just a dot shape it's a 2d array it's 5 12 by 5 12 for now i'm going to loop through these indices which i found above like this and my new slice here's that gets a little tricky i take my region props which i have here region props i take my ith region props so for example here for example i could take my zeroth one the zeroth region prop corresponds to the largest area so i go like this i can get the coordinates of this area like this so here's my coordinates now i want to get this in a way that i can index my new slice so i want to index my new slice which i created and i'm going to take out all these indices i'm going to set it equal to 0 or 1 or whatever the index is so here i have my things but i need to in order to do this i need to transpose it like this and then turn this into a tuple like this so i have a sort of indices that i can index another array then i take my new slice which is 512x512 and i can index using what i created here and now i have all my values at these points and if i want to set this equal to a particular value i'm going to set it equal to i plus 1. so here equals i plus 1. well here i don't have i defined but i could set it equal to one for example and now i can plot that p color mesh my new slice and you see that i've taken this region and i've set all the values equal to one and i'm going to do that for the three largest areas here so i create my new slice i go through the indices corresponding to the three largest areas which is going to be 0 6 and 5 this this and this each setting them to the corresponding separate distinct integer and that's what i do here and i don't have my indices to find i have to run this this and this and then i can plot my new slice and you can see that i have three different regions here each with its own color and i've removed all the smaller little areas again we went through a lot of code there you should sort of go through this slowly break it apart like i did and i think you'll understand exactly what's going on here so now what we're going to do is we're going to do this we did this for one slice we want to do this for all the slices in our image so here's a function called keep top 3. it's basically a summary of all the code we did above we create a new slice which is the same shape as the old slice set equal to zero i take the region props for that slice the areas and the indices and then i loop through the indices corresponding to the three largest areas and that new slice which is all zeros i set each of those equal to integer values like i did above so this takes in a slice the old slice which has all the different areas on it and it returns a new slice which corresponds to separate integer values the top three areas because i have a function that does this i can then vectorize it like i did before well keep top three only takes in a 2d slice that's how i defined it here the top three top three areas takes into 2d slice but i'm feeding it a 3d array called mask labeled and again vectorize if i didn't give it the signature it would try to go through every element of the 3d array but i'm saying no i want you to take in a 2d array and return a 2d array and that's what you're vectorizing here you rotate you're vectorizing a function that takes in a 2d array and so it's an nm to nm that's 2d array in the signature and it takes in mask labeled it's using the keeptop3 function and i get a new mask labeled array and then for example this is this slice here and it does it for all the different slices so all the different slices i have these top three areas like this the next thing we're going to do is we're going to fill in all these small holes in the lungs they're kind of gross and it's easy enough to do with these functions so first of all i'm going to turn this back into a binary array where i just have separate areas in order to do that i just create a mask and that's where a mask labeled is greater than zero and of course if i plot dot if you color mesh this again this is a 3d array this is not going to work but i can give it the 170th slice and now you can see it's just a boolean array of ones and zeros and what we're going to do is we're going to fill in all the little holes of this boolean array which is just trues and falses for that i use this ndi dot binary fill holes function again it's a function that takes in a 2d image that has ones or zeros and it fills in the little holes this is sort of what the image looks like and i'm vectorizing it so it does it for each slice in the 3d array mask that i give it because remember that mask is still a 3d array and i want to apply this to every single slice in that 3d array so i create my mask and then i vectorize this operation which fills in the holes so this is the function i'm vectorizing binary fill holes again if i wanted to i could just call this function on a slice of the mask like mask 170 if i can type here mask 170 and it will do this automatically but here i'm saying okay i want you to vectorize this function and i want you to do it over the full 3d mask but only on 2d slices 2d slice by 2d slice so it'll go off the first axis there of mask so mask.shape it's going to do it for all 263 elements so i get my mask this will fill in all the holes and i can plot that color mesh it and you see that i've gone from an image that has sort of holes in the lungs up here and now the lungs are nice and filled in now in some slices actually you end up taking the trachea and it's sort of annoying so if i look at this for example here sort of i'm above the lungs but i'm looking at a slice corresponding to the trachea and i have a hole here and i don't want that hole so we need a way to remove the trachea so it turns out that in a 512 by 512 ct scan the trachea typically takes up less than 0.69 percent of the area so what we do is we delete all areas that have less than 0.69 percent so again i can use my region props function to take all these separate three biggest areas and delete any areas that are less than 0.69 percent of the total area of this slice so here's my function that removes the trachea i take my slice and c which is the sort of amount of the area that it needs to be less than is 0.0069 times the total size of the image 0.69 percent so here's my new slice i copy the original slice of the 2d array a slice like this then i label the slice this is like what i did before so i label the slice so each of the separate areas will have their own label to them then i get the region properties of each of those labels like before so there's three different areas here because i've enforced that it's only the top three so i have here it will give me the properties of each of these three regions i get the areas like i did before and then my indices now we're going to be well i want to look at the location where these areas are less than c so the 512 squared divided by the areas that's sort of the percent of the image and if it's less than 0.0069 that would be like the trachea then delete it so we're getting rid of any small areas here in the image and that's typically the trachea here so these are the indices corresponding to each of those separate regions then i loop through these indices and at that new slice the locations where corresponding to those areas we set those equal to zero and i return the new slice so going through this slowly here's mask negative 50. i can call the label function on mask negative 50 like before and i get this so i have region 1 region 2 region 3 then i call my region props on this like before and it should be an array of size three that gives me the properties of each of these three separate regions i can then get the separate areas for each of these regions and then i can look at the areas you can see that one of them is big that's the table the other two are small in total area of pixels those are the other two regions and then i can look for the indices where the area divided by 512 squared is less than 0.069 which is the percentage of the image that we're using as a cutoff so if i look at this we'll see is 0.0069 it says that at indices 1 and 2 this holds true and then i loop through these indices in my new slice set it equal to zero so these areas are manually going to be set equal to the background value which is zero and then i'm left getting rid of all any small areas in this image so here's my function i vectorize this function because again i'm going to pass the full 3d mask in but i'm only operating on each 2d slice of the image and it will return remember it takes in a 2d slice it returns a new 2d slice it's going to do this for all 263 slices and return a corresponding 3d array and so now i have a new mask and if i look at this you can see that my trachea is indeed removed in the slices that i want and everything else before is like normal so now it's time to remove the table and we're going to use a clever little trick here to do this now you'll note that each of these separate regions has its own center of mass this is something from physics and you could think of all these little pixels as being like little mass elements so this center of mass might be here this center of mass might be here and then there's one here so what we're going to do is we're going to delete the area that has the lowest center of mass in y that's sort of one way to think about this so let's take an example so here's my labels like i got before and i can plot that color mesh the labels like this and we can compute the center of masses for each of these separate regions so i can call the center of mass function i'm going to take one of these regions so if i look at the plot dot color bar here you can see that this is it corresponds to the integer three this is the integer two and this is the integer one the purple is zero the background so what i'm going to do is i'm gonna say okay let's look at the binary array where labels equals equals three i can plot dot p color mesh this remember this just gives a boolean array labels equals equals three is a boolean array so it'll be one or true here and then zero everywhere else so i can color mesh this and it gives me this thing here from this 2d slice i can compute the center of mass and it tells me where it is in y and x so 185 in x 273 in y here on this image so now we have the delete table function it takes in a 2d slice the first thing we do is we copy the slice then we label it like we did above the slice is going to come from the mask and so it's going to be either ones or zeros we'll label it to have separate regions like this i will take the unique indices that are not equal to zero so we go labels the ampu.unique of the labels here but indexing from one onwards because we don't want to take the background then we'll compute the center of mass for each of these separate regions and then what we do is we loop through the corresponding index so whatever value is given to that region and the center of mass for each of these things and i say if the center of mass and y is less than 30 percent the height up here delete it so that will delete the table and then also if the center of mass is greater than 60 percent the height delete it the reason for that is sometimes you get the nasal bridge which has air in it as well so you get a little air pocket up here so that will delete any nasal bridge sort of thing here just to keep the lungs and so i can run this function we can create a new mask from our old 3d mask then we can look at our new mask and as you can see the table is deleted and it is deleted for all slices here so this is beautiful we have a 3d mask of the lungs the final thing we can do is we can take these lungs which are segmented it's almost exactly where the air is but we can grow them by growing the border a little bit and for that we use the binary dilation function that takes an array of ones and zeros and it expands the border of the areas just a little bit and it depends on the number of iterations you use so here i'm using five iterations we'll use this on our mask new so i'm gonna feed in mask new here and i can plot this and you can see that the area will be just a little bit more expanded so that gives a little more leeway so that your mask just sort of not only takes in the lungs but it's a little bit bigger as well so finally we'll we'll create a 3d image in plotly and now this is a very resolute image so what we're going to do is we're going to reduce the quality by about 40 percent the way to do that is using the zoom function so that's what also what i imported from the nd image library uh sorry from the yeah right here i ported the zoom function what that will do is it will take an image and it will reduce the quality by 40 so that means that the resolution decreases so instead of being 512 maybe it's 200 and something and so uh that's what i do and i get this image here so now if i look up the shape of my image you'll see that it is no longer 512x512 but it's shrunken by 40 so it's the same sort of thing here you just lose resolution quality of the image then we're going to get arrays of x y and z so in a ct scan by the way the x and y uh distance here so in the 512 by 512 is about one-fourth the size of the distance between separate z slices so in order to manage that i create arrays of x y and z and i multiply z by four so if i look at the actual arrays so the x goes like this y goes like this and z is in steps of four this is going to be for plotting in plotly i can create a 3d mesh grid out of these three things again i have a video on this function so then what i do is i create a 3d plot in plotly which i can then open in an html browser and so the data in this is going to be go.volume i feed in my mesh grid here and i make sure to flatten each of the mesh grids so the mesh grid of course is a 3d array if i call flatten on it it will sort of go like that so this is what it originally looks like um i also feed in my uh image here making sure to transpose it a little bit so that it fits the right dimensions for the image i make sure to flatten this as well so all of these are 1d arrays and i can save it and so here by the way i had to open the html file in edge it seems like edge can do this better than chrome so points to edge there but i opened up this uh interactive html from plotly and as you can see we have the two lungs segmented beautifully and of course the main purpose of this video wasn't quite the visualization at the end this is more just a tool for me to show you like all the slices at once but the idea is that once you have a mask of a region of interest it could be a lung for example in a scan it could be the heart it could be any picture where maybe you have a car and you want to segment the car you can do things with those masks right you can multiply for example if i multiplied the whole ct image by this mask if i go back here and i multiply uh where's my mask uh this is my final mask nope it's called mask new well what happens if i multiply mask new times my original image uh let's look at that here my original image is just called image i believe yep if i take my mask new and i multiply it by image i get a new image so i'm going to call this image new and then i can plot that color mesh this my image new and i'm going to take the 170th slice because that's where the lungs mainly were and i spelt it wrong plot dot p color mesh and you can see that now our whole image is um basically uh sort of ignored on the outside and then we're really just focusing in on the lungs here so masks are really useful for things like this and if you're training an ai model for example maybe you want it to ignore everything on the outside and just focus on the lungs while masks are a beautiful way of doing this anyways i hope you enjoyed this video if you made it to the end hopefully there's a few things that you got out of this video for your own research remember to subscribe if you enjoyed this video and i'll see you next time

Channel: Mr. P Solver

Views: 68,877

Keywords: image processing, python, python image processing, image processing using python, digital image processing, image processing python, python image processing opencv, digital image processing with python, image processing tutorial using python, what is image processing, code of image processing, image processing projects, image processing crash course, python images, python image processing pillow, image processing with python, learn image processing using python, segmentation

Id: UIgaLDgb2fY

Length: 31min 49sec (1909 seconds)

Published: Wed Jun 29 2022