× {{alert.msg}} Never ask again
Get notified about new tutorials RECEIVE NEW TUTORIALS

How to divide each lung's region into 4 sections using MATLAB

Ray Phan
Jan 30, 2015
<p>One possible algorithm I have is to take each region and find its bounding box. You would then split up this bounding box into a 2 x 2 grid and extract each subsection manually. Your image that you posted has an unnecessary white border, and so I'm going to use <a href="http://www.mathworks.com/help/images/ref/imclearborder.html" rel="nofollow"><code>imclearborder</code></a> to clear this up. </p> <p>Additionally, I'm going to make your regions completely binary and so I'm going to threshold the image. I chose intensity 32 to escape any quantization errors that resulted. Also, there is a bit of noise at the bottom of the right shape. As such, I did a morphological opening with <a href="http://www.mathworks.com/help/images/ref/imopen.html" rel="nofollow"><code>imopen</code></a> with a square structuring element of 5 x 5. Also, the uploaded image you have there is RGB, even though it appears grayscale. As such, I called <a href="http://www.mathworks.com/help/matlab/ref/rgb2gray.html" rel="nofollow"><code>rgb2gray</code></a> to convert this to grayscale.</p> <hr> <p>Now, back to your question, use <a href="http://www.mathworks.com/help/images/ref/regionprops.html" rel="nofollow"><code>regionprops</code></a> and specify the <code>filledImage</code> property so that you are able to extract the images where each shape is encapsulated in a bounding box. Once we have these images, we can just iterate over them, split up the image into a 2 x 2 grid, then put them into another array. Let's use a <code>cell</code> array for this. The key now is to figure out how to split up this grid. What I would recommend is use <a href="http://www.mathworks.com/help/matlab/ref/mat2cell.html" rel="nofollow"><code>mat2cell</code></a> so that you can essentially subdivide your image to a 2 x 2 grid then place each portion into its own place where each element is a cell array. This will be good because if you can't evenly divide your region into 4 equal regions, putting your portions in a cell array will handle any arbitrary size each region may be.</p> <p>To demonstrate that this works, I'm going to also write code where I generate a figure for each region, then inside the figure, each portion is placed in its corresponding spot in a 2 x 2 grid. Without further ado, here's the code.</p> <pre><code>% //Read in image im = rgb2gray(imread('http://i.stack.imgur.com/PQQox.jpg')); %// Threshold and clear border im = imclearborder(im &gt; 32); %// Do opening to clean up some noise im = imopen(im, strel('square', 5)); %//Perform regionprops fI = regionprops(im, 'FilledImage'); %//Initialize cell array out = cell(1, numel(fI)); %// For each region we have... for idx = 1 : numel(fI) %// Grab the region img = fI(idx).FilledImage; %// Get the dimensions of the region [rows, cols] = size(img); %// Split up the region into a 2 x 2 grid and place into individual %// cells top_half = floor(rows/2); bottom_half = rows - top_half; left_half = floor(cols/2); right_half = cols - left_half; out{idx} = mat2cell(img, [top_half bottom_half], [left_half right_half]); end %// Now create a figure for each region where we place a portion in its %// respective place in a 2 x 2 grid close all; for idx = 1 : numel(fI) figure; img = out{idx}.'; for idx2 = 1 : 4 subplot(2,2,idx2); imshow(img{idx2}); end end </code></pre> <p>The trick with the second <code>for</code> loop is that for each region that I'm looking at, I transpose the 2 x 2 cell array, because this will avoid me having to write 2 <code>for</code> loops to display the portions for each region. We can access elements in a 2D array using a single index, but these are <strong>column-major</strong> indices. The way we can populate the 2 x 2 grid of plots is to use <a href="http://www.mathworks.com/help/matlab/ref/subplot.html" rel="nofollow"><code>subplot</code></a> but this assumes that we are progressing in a row-major format (from left to right, and top to bottom). Column major assumes top to bottom, then left to right. If we transpose our 2 x 2 array of cells, we can populate the plot in row-major fashion using column-major indices.</p> <p>In any case, these are the two figures I get from the above code:</p> <p><img src="http://i.stack.imgur.com/BoVrH.png" alt="enter image description here"></p> <p><img src="http://i.stack.imgur.com/vqJuQ.png" alt="enter image description here"></p> <hr> <p>Note that there are some regions that have more pixels than others.... it's just the nature of how we split up the object within the bounding box of each region.</p> <hr> <p>Now, if you want to access a shape, you can use:</p> <pre><code>shp = out{idx}; </code></pre> <p><code>idx</code> would the shape that you want and <code>shp</code> would be the 2 x 2 grid of cells for each region within the shape. If you want to access a particular region in the shape, you can do:</p> <pre><code>reg = shp{row}{col}; </code></pre> <p><code>row</code> and <code>col</code> are the row and column location of the region that you want within the shape.... so <code>row</code> would be either 1 or 2 and <code>col</code> would also be either 1 or 2.</p> <p>This should hopefully get you started. Good luck!</p> <hr> <h1>Edit</h1> <p>Due to your comments in your post, you are assuming that each shape is roughly shaped like an ellipse, and you want to segment the ellipse into 4 unique regions. We can do this, but it's going to take a bit of hacking, as well as you having the Computer Vision toolbox. Basically, what you need to do is further extend <code>regionprops</code> so that you detected the major axis length, the minor axis length, the centroid, and the orientation. These properties measure elliptical properties of each shape and these tell you how wide and how tall the shape would be and how the shape is oriented in a certain way (using the orientation). We also want to determine where the rough centre of mass is for each shape in your image.</p> <p>After we determine all of these, for each shape that we see in the original image, I want to determine where each of the cardinal directions are located in this ellipse. The shape may be rotated, and so what I need to do is detect the orientation of the shape, then from this orientation we need to figure out where the cardinal directions occur in each shape. You need to know a bit of trigonometry regarding ellipses to know this, but I assume that you don't care about that and want the actual code. If you're curious, I used the standard elliptical equations to determine where the <code>x</code> and <code>y</code> co-ordinates are for the cardinal directions in the ellipse. You can read up about them here: <a href="http://en.wikipedia.org/wiki/Ellipse#Equations" rel="nofollow">http://en.wikipedia.org/wiki/Ellipse#Equations</a> . You start off with the base directions unrotated, then I simply rotate these co-ordinates so that we are in the same orientation as the ellipse. I offset by the centre of the ellipse so that we are in the right location of the image.</p> <p>You can read my comments regarding how I figure this out, but once I figure out these cardinal directions, I draw a black line from the north to south and again from west to east. This is with respect to the rotated co-ordinate system in the shape. To do this line drawing, I use <a href="http://www.mathworks.com/help/vision/ref/insertshape.html" rel="nofollow"><code>insertShape</code></a> from the Computer Vision toolbox. </p> <p>Once I draw these black lines, I erode the image using <a href="http://www.mathworks.com/help/images/ref/imerode.html" rel="nofollow"><code>imerode</code></a> to ensure a disconnection between all of the regions. I then call <code>regionprops</code> one more time with the <code>FilledImage</code> property to extract our 4 regions. I then put this into a cell array, and repeat the process for the other regions. If you want to access a shape and a region within the shape, it's still the same as above. Without further ado, here's the code. I've also modified the figures so that the top shows you what the shape looks like with the major and minor axes drawn in:</p> <pre><code>close all; % //Read in image im = rgb2gray(imread('http://i.stack.imgur.com/PQQox.jpg')); %// Threshold and clear border im = imclearborder(im &gt; 32); %// Do opening to clean up some noise im = imopen(im, strel('square', 5)); %//Perform regionprops and get the shapes s = regionprops(im, 'FilledImage'); %//Get angles for east, north, west and south cardinal_directions = [0 pi/2 pi 3*pi/2]; %//Make room to ensure that when we draw a blank line, it goes all the way %//out of the image buffer_room = 100; %//Store our regions like before out = cell(1, numel(s)); %//Store each shape with all of its subdivisions img_cell = cell(1, length(s)); %//For each shape we have... for k = 1:numel(s) %//Get filled image img = s(k).FilledImage; %// Extract major and minor axes as well as centroid s2 = regionprops(img, 'Centroid', 'Orientation', 'MajorAxisLength', 'MinorAxisLength'); % Convert to grayscale, so we can draw the major and minor axes in black %// Also, insertShape only accepts RGB or grayscale images img = im2uint8(img); %// Obtain the centroid xbar = s2.Centroid(1); ybar = s2.Centroid(2); %// Determine the length of the major and minor axis with respect %// to the centre of the ellipse. Add in buffer_room so that when %// we draw our lines, we make sure we cover the entire shape a = s2.MajorAxisLength/2 + buffer_room; b = s2.MinorAxisLength/2 + buffer_room; %// Create a rotation matrix so that we are rotating our base cardinal %//directions with respect to the major and minor axes theta = pi*s2.Orientation/180; R = [ cos(theta) sin(theta) -sin(theta) cos(theta)]; %// Obtain where each of the four cardinal regions are in an unrotated %//ellipse xy = [a*cos(cardinal_directions); b*sin(cardinal_directions)]; %// Now rotate, then offset by the centroid xy = R*xy; x = xy(1,:) + xbar; y = xy(2,:) + ybar; %// Draw the major and minor axes img = insertShape(img, 'Line', ... [x(1) y(1) x(3) y(3); x(2) y(2) x(4) y(4)], 'Color', 'black', 'SmoothEdges', false); %// Convert back to binary img = im2bw(img); %//Ensure disconnection img = imerode(img, strel('square', 3)); %// Save this for display later img_cell{k} = img; %// Now segment this further into 4 segments img_final = regionprops(img, 'FilledImage'); image_cell = cell(1,4); for l = 1 : 4 image_cell{l} = img_final(l).FilledImage; end out{k} = image_cell; end %// Now create a figure for each shape and display what each region looks %// like for idx = 1 : numel(s) figure; subplot(3,2,1:2); imshow(img_cell{idx}); img = out{idx}; for idx2 = 1 : 4 subplot(3,2,idx2+2); imshow(img{idx2}); end end </code></pre> <p>Here are the figures I get:</p> <p><img src="http://i.stack.imgur.com/aXqJq.png" alt="enter image description here"></p> <p><img src="http://i.stack.imgur.com/kvdU1.png" alt="enter image description here"></p> <hr> <p>Be aware that this is nowhere near perfect. There may be pixels that are outside the boundary of the ellipse, so you won't be able to get all of the pixels. That's kinda what you expect if you try and segment the shape into 4 equal parts like what you have shown above.</p> <p>This is something that you can fool around with though. Have fun!</p> <hr> <h1>Last Edit</h1> <p>You mentioned that you also wanted to grab the original pixels of the bounded regions instead of the binary versions. You can do that by using the <code>BoundingBox</code> property that's part of <code>regionprops</code>. I had to include some additional code so that once you determine where the shapes are located, I also had to determine the location of the bounding box in the image and use this to index into the original image. After this, once we determine where each of the 4 regions are for each shape, I use additional <code>BoundingBox</code> calls to extract where these locations are in the original image.</p> <p>I've added code to also display the binary regions as well as the original ones. The final iteration of the code is below. <code>orig_cell</code> is the variable where you can access the original regions from the original image, and you access it in the same style as <code>out</code>. I've also added in those figures to show you that this works.</p> <p>Good luck!</p> <pre><code>close all; clear all; % //Read in image im_orig = rgb2gray(imread('http://i.stack.imgur.com/PQQox.jpg')); im = im_orig; %// Threshold and clear border im = imclearborder(im &gt; 32); %// Do opening to clean up some noise im = imopen(im, strel('square', 5)); %//Perform regionprops and get the shapes s = regionprops(im, 'FilledImage', 'BoundingBox'); %//Get angles for east, north, west and south cardinal_directions = [0 pi/2 pi 3*pi/2]; %//Make room to ensure that when we draw a blank line, it goes all the way %//out of the image buffer_room = 100; %//Store our regions like before out = cell(1, numel(s)); %//Store each shape with all of its subdivisions img_cell = cell(1, numel(s)); orig_cell = cell(1, numel(s)); img_orig = cell(1, numel(s)); %//For each shape we have... for k = 1:numel(s) %//Get filled image img = s(k).FilledImage; %// Get the original image region bb = round(s(k).BoundingBox); img_orig{k} = im_orig(bb(2):bb(2)+bb(4)-1, bb(1):bb(1)+bb(3)-1); %// Extract major and minor axes as well as centroid s2 = regionprops(img, 'Centroid', 'Orientation', 'MajorAxisLength', 'MinorAxisLength'); % Convert to grayscale, so we can draw the major and minor axes in black %// Also, insertShape only accepts RGB or grayscale images img = im2uint8(img); %// Obtain the centroid xbar = s2.Centroid(1); ybar = s2.Centroid(2); %// Determine the length of the major and minor axis with respect %// to the centre of the ellipse. Add in buffer_room so that when %// we draw our lines, we make sure we cover the entire shape a = s2.MajorAxisLength/2 + buffer_room; b = s2.MinorAxisLength/2 + buffer_room; %// Create a rotation matrix so that we are rotating our base cardinal %//directions with respect to the major and minor axes theta = pi*s2.Orientation/180; R = [ cos(theta) sin(theta) -sin(theta) cos(theta)]; %// Obtain where each of the four cardinal regions are in an unrotated %//ellipse xy = [a*cos(cardinal_directions); b*sin(cardinal_directions)]; %// Now rotate, then offset by the centroid xy = R*xy; x = xy(1,:) + xbar; y = xy(2,:) + ybar; %// Draw the major and minor axes img = insertShape(img, 'Line', ... [x(1) y(1) x(3) y(3); x(2) y(2) x(4) y(4)], 'Color', 'black', 'SmoothEdges', false); %// Convert back to binary img = im2bw(img); %//Ensure disconnection img = imerode(img, strel('square', 3)); img_orig{k} = img_orig{k} .* uint8(img); %// Save this for display later img_cell{k} = img; % //Now segment this further into 4 segments img_final = regionprops(img, 'FilledImage', 'BoundingBox'); image_cell = cell(1,4); orig_img_cell = cell(1,4); for l = 1 : 4 image_cell{l} = img_final(l).FilledImage; bb = round(img_final(l).BoundingBox); orig_region = img_orig{k}(bb(2):bb(2)+bb(4)-1, bb(1):bb(1)+bb(3)-1); orig_img_cell{l} = uint8(img_final(l).FilledImage).*orig_region; end out{k} = image_cell; orig_cell{k} = orig_img_cell; end %// Now create a figure for each shape and display what each region looks %// like for idx = 1 : numel(s) figure; subplot(3,2,1:2); imshow(img_cell{idx}); img = out{idx}; for idx2 = 1 : 4 subplot(3,2,idx2+2); imshow(img{idx2}); end figure; subplot(3,2,1:2); imshow(img_orig{idx}); img = orig_cell{idx}; for idx2 = 1 : 4 subplot(3,2,idx2+2); imshow(img{idx2}); end end </code></pre> <p><img src="http://i.stack.imgur.com/YYGkI.png" alt="enter image description here"><br> <img src="http://i.stack.imgur.com/hTSau.png" alt="enter image description here"></p> <p>This tip was originally posted on <a href="http://stackoverflow.com/questions/26784968/How%20to%20divide%20each%20lung's%20region%20into%204%20sections%20using%20MATLAB/26785731">Stack Overflow</a>.</p>
comments powered by Disqus