Image J-ing. AP 186 Activity 4: Area estimation

To avoid further negativity for my blogs, let’s rename my posts to

starting to become a promising academic endeavor:

AP 186 edition.

In the third installment, I will be discussing about our activity called area estimation. In the age of computers, a lot of analysis in science can already be accomplished mainly on our devices. Gone are the days where the only option to perform graphical analyses is by hand.

I consider myself blessed since I was born in the age of technology, although I believe that not a lot of people are sufficiently aware of what the limits of our devices are (but then again, I may just be shined with a bit more light than others). At the very least, I hope that I may be able to shine a bit of light to people who might be able to read this.

Going back to the topic: area estimation.

 

Part 1. Green’s function

It’s best if we start with the discussion of Green’s functions for now. The simplified equation gives equation for obtaining an area given the Green’s function is

x and y are your friendly neighborhood Cartesian coordinates. We do not expect current programs to fully comprehend the complexity of continuous variables, so we recast the equation above to a digestible, discretized form,

Some things to note:

• N here would be the number of pixels that constitute the edge of the image, x and y would refer to the horizontal and vertical pixel coordinates, respectively. As with previous activities, the origin is at the top left corner of the image.
• The sum is cyclic, as the integral in the 1st equation is a closed contour integral. At the Nth iteration of i, the 1st pixel is designated as the (i+1) pixel.
• The order of the pixel in the iteration of the sum is important. The progression of the contour integral follows a specific trace of the contour (counter-clockwise positive, clockwise negative), so it would make sense for the discretized sum to be the same. The code later will give light on this.

Part 2. General code used

 image = imread("insert image filename here")

Nothing fancy here; just load the image to scilab. The image format used for now is monochrome bitmap to avoid having a color image.

//finds x,y coordinates of edges in the image (contour of image)
[x0,y0] = find(edge(image)) //for geometric figures
[x0,y0] = find(edge(image,'prewitt')) //for land area

Next part finds the edge contour of the shape needed to be analyzed. We only need the boundaries of the figure and nothing else. x0, y0 are the edge coordinates.

Note that the real land area uses a prewitt edge finding algorithm for the google map edge later. See discussion there, not here.

//corrects x and y coordinates with respect to the centroid
x_cent = sum(x0)/length(x0);
y_cent = sum(y0)/length(y0);
x = x0 - x_cent; y = y0 - y_cent;

Standard pixel coordinate systems use the top left corner of the images as the origin, so translations according to the centroid are performed. Given a plane of uniform density but arbitrary shape, the centroid is the point in the plane coincident to the center of mass, and can be expressed as a function of 2-D Cartesian coordinates as:

N is the number of pixels, and x_i, y_i are the coordinates in pixel space. We need to center the origin on the centroid for more accurate calculations of the area.

 

//generates theta list in polar coordinates
 r = [];
 theta = [];
 for i = 1:length(x)
 theta(i) = atan(y(i),x(i));
 r(i) = sqrt(y(i)**2 + x(i)**2)
 end

Remember that the contour integral in Green’s theorem requires that the progression of the sum of the integrands follow the trace of the contour integral. To ensure that this happens, proper order with respect to the angular variable, is required. For this purpose, the code above ensures that I first have the list of the values of and . The radius variable is more for completion and verification of the radius of certain geometric figures.

 

//obtains correct order of theta
 theta_corr = []
 [theta_corr, k] = gsort(theta, 'g', 'i')

This part of the code is pretty important, as it sorts the values from lowest to highest.

A little something about the gsort function (c/o scilab documentation):

gsort is a quick sorting algorithm. For my purposes…

inputs:

• list to sort – theta
• type of sorting to perform – ‘g’, or sort everything. (Note that gsort can also organize lists that have more than 1 dimension)
• order – ‘i’, or increasing

outputs:

• theta_corr – here be my list of proper order
• k – here be my list of indices of the previous unsorted list. Very useful for next step

 

//sorts x and y lists according to increasing theta and applies Green's formula
 area_green = 0;
 x_corr = x(k);
 y_corr = y(k);
 for i = 1:length(x)-1
 area_green = area_green + x_corr(i)*y_corr(i+1) - x_corr(i+1)*y_corr(i)
 end

area_green here is our sum variable, and will be responsible for giving us the sweet sweet area. x_corr and y_corr are the sorted lists of the Cartesian coordinates according to increasing theta (see how magical k is; make it an argument for the original list and it will sort accordingly. Thank you, Carlo Solibet for this piece of juicy info). I then proceed to sum the contributions of the integrands to the sum, but take note that I only iterate until the term due to the fact that the sum should be cyclic: term is the first term.

//adds the final term to the sum (when i = max value)
 area_green = area_green + x_corr(length(x))*y_corr(1) - x_corr(1)*y_corr(length(y))
 area_green = area_green*0.5
 area_green_scaled = area_green*(1/9.8)**2 //for area scale

Here we add the last term, and since the sum should be halved, I did so as well. Last term to be discussed later, as this applies to the area estimation in google maps part.

//generates relative deviation of calculated area to actual area
 area_actual = %pi*(mean(r))**2
rel_dev = 100*abs(area_green - area_actual)/area_actual

No new method testing goes to work not bringing its relative deviation values with him/her!

 

Part 2. Estimation of a geometric figure

It’s a circle.

Now, we use the code to test it on a circle, the official 2-D geometric mascot of physics. For this part, I created a 1000 x 1000 pixel grid of zeros, and drew a circle with a radius that is 60% of the length of half the side length. I centered the circle in the middle of the grid for simplicity, and the code is right here

//Generates the circle image
r = sqrt(X.^2+Y.^2);
A = zeros(nx,ny);
A(find(r<0.6)) = 1;

I’ll refer to the units of length here in px (for pixels), and the area px^2. Theoretically, the radius of the circle should be equal to 300 px, and the area should be equal to . But due to the resolution of the image, I took the mean of the r values (ha! it has uses), and used it to calculate the ‘actual’ area of the circle subjected to pixelation in a 1000 x 1000 grid and resulted to a value of  . For the area calculation using the discretized Green’s function, I obtain , resulting in a 0.0023% relative deviation with the actual area of the pixelated circle.

This method must be pretty neat. Now let’s try it with something a little bit more complicated. Real areas.

Part 3. Land area estimation

The next part of the activity deals with the estimation of the land area of a certain building. I originally aimed for estimating the UP Diliman Church of the Holy Sacrifice, but part of the prerequisites for choosing an area is that the geometry has to be a bit complicated.

 

Complicated land area filed under: wishful thinking.

I accomplished and submitted this activity late, so someone already chose the National Institute of Physics (Robbie…). So to please my girlfriend, who is most probably not going to read this blog anyway, I chose the UP College of Law!

It’s not a circle!

It is important that the scale is determined. How many meters to a pixel. From the tiny scale in the bottom right corner of the image, it was determined to be 1 meter per 9.8 pixels. Recall that the scaling is important in this part of the code:

 area_green_scaled = area_green*(1/9.8)**2 //for area scale

Modifications to the original image have to be made, namely to convert the image to a simple black-white (BW) image below. For this, MS Paint has my gratitude.

Primed for editing

Now, we find the edge coordinates. It wasn’t as simple as with the geometric figures due to the rough edges from the conversion to a BW image. This is why I used the prewitt algorithm to ensure that a full coverage of the edge. But I mostly found it via trial and error.

Unfortunately, I was not able to obtain the actual and area of Malcolm Hall (UP law), but I was able to obtain an area of . That’s a third of a hectare, which I find to be somewhat reasonable.

Part 4. Image J

Now we come to arguably the most perplexing section of the activity. Image J!

For future purposes of area or length estimation, Image J will be my best friend. It makes area or length estimation easy as you only need to tinker with the UI of the program. As with Scilab, the program is free.

Model for the day: beep card.

Our example for today will be the LRT/MRT beep card. More power to mass transit. I chose this card to make the example rather easy.

Setting up the scale.

First step is to measure the actual length of the card, which was 85.6 mm. Then use the Image J’s digital measuring tape to get its distance in pixels (3064). Set the unit of length to mm and voila, you get your scale: 35.7944 pixels/mm.

Featuring Knock Knock by Nikki Yanofsky.

For testing of the scale, I obtained the measurement of its diagonal. Actual measurements led to 99 mm, whereas the digital length was obtained to be 98.342. Pretty darn close. 0.66% close, in fact.

Part 5. Evaluation

9/10. I was able to fulfill all of the required tasks for this activity, although I wish I explored more geometric objects, or tested other features of Image J.

Part 6. Acknowledgements and references

1. Soriano, M. Length and Area Estimation in Images. Applied Physics 186 Activity Manual.
2. Equation and explanation for the centroid. Taken from http://www.ce.memphis.edu/3322/Pdfs/PaulsPDFs/Centroids%20and%20Moment%20of%20Inertia%20Calculation.pdf.
3. Photograph for the Chapel of the Holy Sacrifice. Taken from https://www.pinterest.com/pin/563231497122075996/.
4. Carlo Solibet, Mich Cirunay, and Roland Romero for their assistance in performing this activity.