AP186 Activity 8 Blog Report - Morphological Operations - Part 2

Here comes the application part! Ma'am Jing provided us with two images of scattered punched paper, one with circles of the same size (relatively), and another with larger circles among the regular circles to simulate cancer cells.

It's 10:00PM and I'm a bit tired, but so far, I've cropped seven images from the image of just regular circles. I did that using GIMP.

Fig. 1. Digital image of scattered punched paper simulating cells.

Fig. 2. Seven subimages covering 256 by 256 pixel areas of all cells in Fig. 1.

Fig. 3. Digital image of simulated cells with some simulated cancer cells.

---

(11/23/2016, 9:55AM PhST) I was working on my code yesterday afternoon as well as this morning, and so far I think I'm doing alright.

Anyway, here's a quick recap on my progress so far. I checked for the best segmentation threshold using the histogram of the first subimage, which looks like this:

Fig. 4. Number of pixels per grayscale value of the first subimage.

I settled with the threshold of 200, making all pixels with values greater than 200 as 1 (white) and others 0 (black), since I found that a threshold higher than 200 gets some information from the cells lost, while the artifacts can be removed by morphological operations.

So for all seven subimages, I segmented using a grayscale value of 200. I then tested the effects of a few morphological operations, including CloseImage, OpenImage, and TopHat, and with different circular structuring elements of varying radii. I also saw that a 'circle' structuring element is one of the possible parameter values in place of 'custom' from this source:

http://forge.scilab.org/index.php/p/IPD/source/tree/9/macros/CreateStructureElement.sci

I picked:

SE = CreateStructureElement('circle',8);
I(:,:,i) = OpenImage(BW, SE);

And from the first subimage, it seems to successfully remove artifacts and round out some of the broken circles. I then used SearchBlobs, which assigns a number to each blob:

I(:,:,i) = SearchBlobs(I(:,:,i));

I then use this method to make the separate blobs clear through different shades of gray:

imwrite((255.0/double(max(I(:,:,i))))*I(:,:,i),'S_C1_0'+string(i)+'_sb.png');

I tried the entire process for all, and these are the results:

Fig. 5. Segmented subimages (left) using t=200 and the results (right) after OpenImage and SearchBlobs, with each numbered blob made to show at a different shade of gray.

For the third and sixth subimages, the threshold of 200 had a lot of artifacts. However, the morphological operation seems to have removed the artifacts aptly.

I then got the areas by iterating between the slices of my multimatrix containing all labelled SearchBlobs. The maximum of each slice should also give the total number of blobs for that slice, so using that as the end point of my iterating index should allow me to access each numbered blob from 1 to that number. Accessing them and placing the calculated area in an empty array I called "areas" gives all the values I need.

for i = 1:7
    for j = 1:max(I(:,:,i))
        area = size(find(I(:,:,i) == j));
        areas($+1) = area(2);
    end
end

 I checked the values compiled in the "areas" array. It seems that there are blobs with area of zero, and this I simply removed by making a new list of areas "areas2" and finding all that aren't zero:

for i = find(areas > 0)
    areas2($+1) = areas(i);
end
histplot(100,areas2);

The histogram, with 100 bins, of all the pixel areas of each blob turned out like this:

Fig. 6. Histogram of the areas, with x-axis being number of pixels.

I am still trying to interpret the physical equivalent of the above histogram's y-axis. This is because histplot doesn't seem to give the number of blobs within that bin. I will also get the mean and standard deviations of "areas2" shortly.

However, what I know so far is that this new array with those odd zeros removed, "areas2" lists 74 areas when I checked its size. That means there are 74 blobs total for all subimages.

(12:29PM PhST) I found out how to make the histogram not normalized, so that it gives the number of blobs within each bin. The default of plothist is normalized, so I just have to set the parameter for normalization as false instead of true. (Reference: https://help.scilab.org/docs/6.0.0/en_US/histplot.html)

Anyway, below is a better histogram, with 32 bins. Each bin has a 100 pixel range, and the x-axis ranges from areas of 0 pixels to 3200 pixels.

Fig. 7. Histogram of the areas with 32 bins, and the y-axis being the number of blobs within each range of area values.

From the above histogram, the mode should be at the 500 to 600 pixels range, with 39 out of 74 blobs falling within that range. But since I need to get the mean, I will calculate that using Scilab.

-->mean(areas2)
 ans  =

    782.3783783783783

-->stdev(areas2)
 ans  =

    646.4978721575065 

These values, particularly the standard deviation, are too skewed to the overlapping groups of cells forming large blobs. However, even with the given morphological operations, there doesn't seem to be a combination that eliminates those overlapping groups without eliminating the single cells first. So I don't think a change in morphological operation is the way to go here.

I think that, since the histogram already gives an idea on the outliers, I should calculate the mean and standard deviation with those outliers removed. Particularly, those above 650 pixels in area. The main goal of this application part is to identify which cells are the cancer cells, so getting a better set of data is key.

-->mean(areas3)
 ans  =

    481.3454545454546

-->stdev(areas3)
 ans  =

    108.5340394090154 

This range is much better for extracting cancer cells, I think. I will think more about this later, since I have class.

(4:27PM PhST) I will try using the same morphological operations on the image with cancer cells, then pick the blobs that are above the mean plus the standard deviation I obtained.

---

(11/25/2016, 10:08AM PhST) Using 590 (mean plus standard deviation) as the threshold and the same grayscale threshold and morphological operation as the subimages with just normal cells, these are the images I obtained:

Fig. 8. Image with simulated cancer cells among normal cells, segmented.

Fig. 9. Image with artifacts removed through OpenImage with a circular SE of radius 8.

Fig. 10. The previous image with different grayscale values for each numerically assigned blob.

Fig. 11. Blobs that exceed the range of areas within the mean plus standard deviation.

 The above image still has cells that overlap with each other. Thus, we need a method that removes those that aren't cancer cells. Since we still have the structuring element from before, we can use it for erosion twice and then dilation twice.

I do this through these lines of code:

cancer = bool2s(ErodeImage(bool2s(ErodeImage(cancer,SE)),SE));
cancer = bool2s(DilateImage(bool2s(DilateImage(cancer,SE)),SE));

with the matrix named "cancer" being Fig. 11 and SE being the same circular structuring element of radius 8 used in the previous morphological operations. This leaves us with the actual cancer cells:

Fig. 12. Filtered cancer cells from the previous image.

I'm not sure if there is a better way to extract the cancer cells in Fig. 12, but I am happy that I was able to pinpoint them from the rest.

Thank you to Albert Yumol, Zed Fernandez and Robertson Esperanza for helping answer questions.

Self-Evaluation (for both Activity 8 Part 1 and Part 2): 8/10

AP186 Activity 8 Blog Report - Morphological Operations - Part 1

Morphological operations seem to include logical operations on sets, but with new concepts such as erosion and dilation. For Part 1 of this activity's blog report, I will try to explain what happens to some logical operations. This would be based on how I understand what I've read from the introduction of this activity's manual.

∈: Something being an element of another. For instance, a ∈ A with a being a point and A being a shape means that a is inside A.

∉: Something not being an element of another. For instance, b ∈ A with b being a point and A being a shape means that b is outside of A.

---

(11/4/2016, 11:23AM Philippine Standard Time or PhST)
For the first part of the activity, we are tasked to hand-draw:

1. A 5×5 square
2. A triangle, base = 4 boxes, height = 3 boxes
3. A hollow 10×10 square, 2 boxes thick
4. A plus sign, one box thick, 5 boxes along each line.

We are then tasked to draw the result of a) erosion by the following structuring elements, and (b) dilation by the following structuring elements:

1. 2×2 ones
2. 2×1 ones
3. 1×2 ones
4. cross, 3 pixels long, one pixel thick.
5. A diagonal line, two boxes long, i.e. [0 1; 1 0]

I am currently in the process of drawing these, before I proceed with the next part of the activity. I will scan my drawings and edit them here.

---

(11/9/2016, 11;19AM PhST)

1. 1-5. (a) & (b)

2. 1-5. (a) & (b)

3. 1-5. (a) & (b)

The reference I used was:
https://www.mathworks.com/help/images/morphological-dilation-and-erosion.html

I will finish 4. 1-5. (a) & (b) shortly.

---

(11/11/2016, 9:43AM PhST) I am done with drawing! I will now proceed to installing IPD and checking if my drawings are correct through Scilab. Once I am also able to access our lab's scanner, I will also scan the predicted output I drew for each and include those images in this post.

(10:01AM PhST) I have successfully installed IPD, thanks to Robertson Esperanza. He also helped me with a small code that can be placed at the start of every .sce file that uses IPD, basically a one-liner that executes the IPD to be commented out after one execution. I will now proceed with writing the code for checking erosion and dilation.

(10:43AM PhST) It turns out that my predictions are off. Way off.

My predictions did not take into account the origin used by erosion and dilation. I now realize that I was eroding more than necessary and dilating more than necessary. I thought that even if one tile of the structuring touched the base shape, that would be a valid overlay, and the entire area covered by the structuring element would become zero (erosion) or one (dilation) if the area has at least one zero.

It turns out that the overlaying works like this, I think:

For erosion, each "zero" tile of the base shape simply gets overlayed with the designated origin tile of the structuring element, and then that entire area becomes zero if the area has at least one "one".

For dilation, each "one" tile of the base shape simply gets overlayed with the designated origin tile of the structuring element, and then that entire area becomes one if the area has at least one "zero".

At least now I know how it works. 
EDIT (11/16/2016, 12:26PM PhST): The above is an incorrect interpretation.

I hope my drawings, although incorrect, fulfill the requirement to predict the output.

Anyway, the designation of the structuring element's origin tile by Scilab is not clear, even when I asked Robbie. However, for each structuring element, I observed a particular origin tile, which I will enumerate here later. For now, I will eat lunch.

1. 2×2 ones
2. 2×1 ones
3. 1×2 ones
4. cross, 3 pixels long, one pixel thick.
5. A diagonal line, two boxes long, i.e. [0 1; 1 0]

---

(11/16/2016, 9:28AM PhST)

After reviewing for the quiz as well as actually taking the quiz, I was informed of how erosion and dilation actually work. Erosion takes the origin points of the structuring element being moved around (but not rotated) to fit inside the base shape. I got that part of the quiz correct.

However, for dilation, it takes the origin points of the inversion of the structuring element (inverted based on its origin) being moved around to intersect with the base shape. For that, I failed to perform inversion before checking for the intersecting relocations of the structuring element.

For the structuring elements, this is how I implemented them in Scilab:

SE1 = CreateStructureElement('custom',[ %t %t; %t %t]);
SE2 = CreateStructureElement('custom',[ %t %t]);
SE3 = CreateStructureElement('custom',[ %t; %t]);
SE4 = CreateStructureElement('custom',[ %f %t %f; %t %t %t; %f %t %f]);
SE5 = CreateStructureElement('custom',[ %f %t; %t %f]);

The locations of the origins as I observe them are as follows:

2x2 ones: the origin seems to be at the lower-right corner "one".
2x1 ones: the origin seems to be at the "one" on the right.
1x2 ones: the origin seems to be at the "one" on the bottom.
cross, 3 pixels long, one pixel thick: the origin seems to be at the center.
A diagonal line, two boxes long, i.e. [0 1; 1 0]: the origin seems to be at the empty "zero" at the lower-right corner.

The last one is a problem since it uses a disjointed origin at a "zero" instead of an origin inside the structuring element.

To make the structuring elements consistent, Ma'am Jing suggested that the origins of the structuring element be placed in the upper-left corner point of the structuring element. This should apply for the structuring elements 1, 2, 3, and 5 (with the fifth being at the empty upper-left corner). For the fourth structuring element it should be okay that the origin is at the center.

To make structuring elements with those desired origins, I will use the method suggested by Ninya Zambale in class. It uses a 3x3 matrix of zeroes as the size of the structuring element, then the structuring element is positioned there such that the desired origin is at the center of the 3x3 matrix.

Through that method I have the following structuring elements:

SE1 = CreateStructureElement('custom',[ %f %f %f; %f %t %t; %f %t %t]);
SE2 = CreateStructureElement('custom',[ %f %f %f; %f %t %t; %f %f %f]);
SE3 = CreateStructureElement('custom',[ %f %f %f; %f %t %f; %f %t %f]);
SE4 = CreateStructureElement('custom',[ %f %t %f; %t %t %t; %f %t %f]);
SE5 = CreateStructureElement('custom',[ %f %f %f; %f %f %t; %f %t %f]);

I have a problem, though. The results seem to indeed place the origin at the correct point of the structuring element. However, even with Ninya's method, the structuring element seemingly is not inverted when it performs the dilation operation. This is problematic, since the simulated isn't the actual expected dilation result. I am still looking for a way to make the dilation correct, in that the structuring element is actually inverted. The results are also very small images, so before I painstakingly enlarge each of the 40 images manually, I have to make sure that they are correct first.

(12:17PM PhST) The new plan is to create the inverted counterparts of the structuring elements as their own structuring element, then use those inverted SEs in dilating the base shapes in Scilab. It seems that it is Scilab's problem in that Scilab does not perform inversion on the structuring element before it dilates the base shape. That is why the easiest fix would be to make those inverted structuring elements myself.

---

(11/21/2016, 7:10PM PhST) The new plan seems to have worked! I will try resizing them with a batch resizer so that I can post them here. But yeah, the dilation and erosion that resulted from the new plan follows the expected dilation and erosion operations from Ma'am Jing's PDF.

(7:24PM PhST) I used http://picresize.com/ and the resizing seems to have worked. The images seem to have lost their relative sizes to each other though, so one tile in one image is not the same pixel size as one tile in another image. That's because picresize only has percentage resizing for shrinking, and none for expanding. I'll try a different website.

(7:40PM PhST) I tried https://bulkresizephotos.com/ but the resized images turned out blurry. That's no good. Let's try something else.

(8:24PM PhST) It seems a more mathematical solution that actually uses Scilab was in order, haha. I'd like to thank the answerer, Julius, to this Stack Overflow question:
stackoverflow.com/questions/35540459/enlarge-matrix-proportionally-in-r

He gave me the idea to use this on my imwrite lines:

imwrite(R11A.*.ones(10,10),'E:\AP186 - Activity 8\R11A.png');

Basically, the unique properties of the Kronecker product are used to enlarge the matrix proportionally, which is perfect for my purposes. A .*. B in Scilab returns the Kronecker product of matrices A and B. So the resulting output image by imwrite for the above code would give the image ten times its original size, and no blurring! Yay!

I'm really happy with the result. Crisp, enlarged images of the results for erosion and dilation. Finally. Here they are, in all their non-blurry glory (each "one" or "zero" tile is enlarged to be 10 by 10 pixels):

Fig. 1. Base shapes upon which erosion and dilation were performed.

Fig. 2. Original structuring elements and their manually inverted versions, with origins at the center. 

Fig. 3. The first base shape when eroded (left) and dilated (right) with each of the five structuring elements.

Fig. 4. The second base shape when eroded (left) and dilated (right) using the five structuring elements.

Fig. 5. The third base shape when eroded (left) and dilated (right) using the five structuring elements.

Fig. 6. The fourth base shape when eroded (left) and dilated (right) using the five structuring elements.

Basically I can sum up erosion and dilation with this visualization:

Imagine stencils that outline the base shape on top of a clean sheet of paper. You have a wooden block the size and shape of the structuring element with paint at the origin. When you erode, you move the block around inside the stencil, and where paint from the origin leaves a trail, that's what's left of your shape. When you dilate, you get the inverted version of the block, and you stamp each position of the block where it still intersects with the stencil shape, and where paint from the origin gets stamped, that's your new shape.

Anyway, once I finally scan my failed drawings of how I expected erosion and dilation would look like, then that would be the end of Part 1. Stay tuned for those scans. Meanwhile, the first few logs of Part 2 should be posted shortly after I write this sentence.

---

(11/22/2016, 2:00PM PhST)

Here are my scanned predictions when I still had a different understanding of erosion and dilation, as promised. The order for each set of scans is, given the base shape, erosion using the first structuring element, dilation using the first structuring element, erosion using the second structuring element, and so on until dilation using the fifth structuring element. Whenever I go to the next sheet of graphing paper due to lack of space, I redraw the original base shape for reference, which is why it repeats in Fig. 9 and 10.

Fig. 7. Scanned predictions for the first base shape under erosion and dilation.

Fig. 8. Scanned predictions for the second base shape under erosion and dilation.

Fig. 9. Scanned predictions for the third base shape under erosion and dilation.

Fig. 10. Scanned predictions for the fourth base shape under erosion and dilation.