The last activity for our Applied 186 course, this one let us take our first step into the world of video processing. We’re to observe a kinematic phenomenon on camera, and compute for a kinematic constant. I grouped up with some friends, and we decided to keep it simple by calculating the acceleration g due to gravity from a video of a falling ball.
I used scilab 5.5.2 for the coding, in conjunction with Paint, WPS Spreadsheets, and PhotoScape for animating the gifs.
Okay, so we took a video of a ball falling and broke it up into frames. It looked something like this.
Look at it go. The ratio of pixels to meters at the full resolution I operated on was 387 pixels for every 0.5 meters, as measured using Paint. The camera operated at 59 frames per second.
Anyway, my plan was simple:
Perform colour segmentation on each of the frames like in Activity 7 to isolate the ball,
Apply Activity 8′s morphological operations to clean up the resulting blobs,
Determine the centroid of the blob like in Activity 4 and get the y-coordinates of these centroids,
Get the instantaneous speeds of the ball using the change in y-coordinates between two frames,
Use these speeds to calculate the average acceleration of the ball,
Convert this into SI units (m/s^2) with the ratio I got earlier, and
Compare this to the theoretical value of 9.8 m/s^2.
Here’s an example of the segmentation and cleaning done on one of the frames.
This particular frame had the centroid of the ball blob located at a y-coordinate of 91.5.
Here’s how the final segmented frames looked all together.
Things went rather smoothly this time, and before long I has a list of y-coordinates ready for operating on. I worked using WPS Spreadsheets and got an average acceleration of 10.72 m/s^2, which differed from 9.8 m/s^2 by 9.4%.
Unacceptable, I said to myself. I then turned to knowledge gained in our Physics 191 class on how to discard outliers in our data. Testing for varying thresholds of deviation, I managed to discard two acceleration values. Computing for the new average acceleration gave me 9.80 m/s^2. Acceptable.
I had to run through this one quickly, since the end of the semester (and te whole year of 2016) meant many more deadlines, both academic and non-academic, were catching up to me. Still, I did as the activity asked, and managed to apply a lot of the knowledge we gained in this course, and even from Physics 191. I’ll give myself an 11 out of 10 for this.
And my final acknowledgements for this course:
Joshua Abuel, Anthony Fox, Kit Guial, and Roland Romero for allowing me to join them in recording the video for this activity
Dr. Maricor Soriano, for guiding us through the course; it was a blast
Dr. Wilson Garcia, for his Physics 191 class where we learned methods on dealing with data
We’ve used histograms previously in image segmentation, but in this activity I learned they can also be used to enhance both monochrome and coloured images. Now when I hear “enhanced images,” I can’t help but think of the crazy image enhancement stuff they used to do in classic CSI, so I was definitely looking forward to this activity.
Grayscale Enhancement
I know we’ve all experienced excessive shadows in pictures. They can come from the shot being in a poorly-lit area, from a side effect of using the flash without a reflector or diffuse, or from backlight. Here’s a photo illustrating the first example.
We may never know what lies in the void I was sticking my hand into; or so I thought at first, but it turns out those shadows still contain information. We just have to force it out, and that’s where histogram manipulation comes in.
First, we test it on the image in grayscale. Since we only have one channel here, the normalized histogram of the image corresponds to its probability distribution function or PDF.
From this, we use the cumulative sum function to allow us to obtain the image’s cumulative distribution function or CDF.
This is what we will use to modify the image’s histogram. We’ll want to pick a CDF to backproject onto the PDF. For this activity we’ll try both a linear CDF and a nonlinear sigmoid CDF, generated using an exponential Gaussian function.
And with no further ado, here are the results from the modified PDFs:
I’d say that while the linear CDF produces a brighter image, it seems a bit overwhelming in its brightness. The sigmoid CDF gives an image that doesn’t go as bright, but it reveals enough detail.
Colour Enhancement
Alright, let’s get to it. The modification we made to the image’s histogram on the grayscale level was basically modifying the “intensity” of each pixel, with high intensity being closer to white.
For coloured images, we use our knowledge of normalized colour coordinates (NCC) to obtain the intensity I, as well as the previously used r, g, and b values for each pixel in the image.
From here, it’s just a matter of applying histogram modification to the intensity matrix of the image, then using it to compute for the new RGB values of the matrix from NCCs r, g, and b. I used both the linear and sigmoid CDFs for this, and the results look pretty good.
I had to try it on something else, and so I found this picture of my niece’s furby I had for.... some reason.
Working on this was pretty straightforward. I had a problem when I tried to apply my grayscale code to the coloured image, but somehow I managed to fix it by starting the code again from scratch. I assume I just missed something while modifying it. I delivered what the activity required, and my images are all there so I say I get a 10 out of 10 for this one.
Here are my acknowledgements and references:
Roland Romero, for helping me troubleshoot my code
Last activity, we used morphological operations and the IPD kit by Dr. Galda to identify abnormal “cells” in an image, but did you know we can use this kit in a variety of other ways? In this activity, we were made to let scilab read and play sheet music. It wasn’t an easy task, but it was pretty fun to work on.
The sheet music I used for this activity is a piano adaptation of the song ECHO by Crusher-P. The piano piece can be found here, while the original song can be found here. For this activity, I only included the melody following the treble clef up until the 64th measure of the piece.
First of all, we should know that scilab can produce sounds from sinusoid functions.
The frequency f of the wave here determines the pitch of the sound, and the time t defines how long the sound plays.
Now the product of this is called pure tone, a continuous sound that doesn’t quite sound like a note played by a piano. To fix that, we applied an sound envelope or an ADSR (attack, decay, sustain, release) envelope to emulate the sound of an instrument.
Here’s the note function I used, based on the one provided by our professor, Dr. Soriano:
//note function
function n = note(f, t)
n = sin (2*%pi*f*t);
//defines the envelope
env1 = linspace(0, 1, round(0.025*length(t))); //attack
env2 = linspace(1, 1, round(0.05*length(t))); //sustain
env3 = linspace(1, 0.9, round(0.05*length(t))); //drop before decay
env4 = linspace(0.9, 0.45, round(0.85*length(t))); //decay
env5 = linspace(0.45, 0, length(t)-length([env1,env2,env3,env4])); //release
env = [env1, env2, env3, env4, env5];
//applies the envelope
n = n.* env;
endfunction;
So with that done I had two main objectives: identifying the kinds of notes (half, quarter, eighth, sixteenth), and their pitch, in addition to recognizing rests and dotted notes.
I stitched together the many lines of music to form one long continuous line from the first to 64th measure. This was my input image.
As I said before, I used Dr. Galda’s IPD kit to analyze the notes here. After reading the image, I realized I had a problem. If I applied morphological operations to isolate the quarter notes, I’d be removing the half notes. And even so, how would I differentiate between quarter, eight, and sixteenth notes?!
I ended up using three segmentation methods for the notes:
one to detect the “filled” note heads that quarter, eight, and sixteenth notes shared,
another to detect the “hollow” half note heads, and
the last to detect the bars and flags of eighth and sixteenth notes
Here’s an example of how the first segmentation process worked.
The other two used different thresholds and morphological operations.
With the notes indexed, I’d get the x and y-coordinates of the centroids of the note heads and used the x-coordinates to control the sequence at which the notes would be played, and the y-coordinates to determine the pitch of the note. Additionally, if the x-coordinate of a note head coincided with those included in a bar or flag, the note was identified as an eighth or sixteenth note.
For the rests and the dotted notes, I had to use two more segmentation processes.
Fortunately, all rests are located in the same spots on the musical staff, so if a blob was found to occupy that spot, I could mark it as a corresponding quarter or eighth rest.
For the dots, I had a bit of trouble since most of my notes were staccato notes. To remove the staccato marks, I only considered dots that didn’t line up vertically with other notes. After finding the notes closest to these dots, I applied the corresponding modification to their time (dotted notes are played at a total time equal to 1.5 times the time of the usual note).
I... ignored the staccato notes for now, as modifying the envelope either produced no noticeable effect, or too much of an effect on the final notes in the wave file.
With all of that done, I have my final output.
You can listen to it here. HEADPHONES WARNING!! It can get kinda loud.
That was a really tedious activity, but I can’t deny it was a lot of fun too. I met the requirements of the activity, plus adding in the sound envelope, together with recognizing rests and dotted notes. For that, I think I deserve a 12 out of 10.
Now I get to sing praises to the following as my acknowledgements:
Dr. Maricor Soriano, our professor who provided the note function
Joshua Abuel, for helping me with the process of detecting notes
Physics of Music Notes, for the note frequencies
Luman Magnum, for the info on how to apply a sound envelope
Crusher-P, for creating the beautiful and addictive song ECHO
jazzermazzer99, for transcribing a piano adaptation (and making an easy version) of ECHO
Yes, bring back the Power Rangers references.
In this activity, we’d learn about morphological operations, specifically erosion and dilation, and apply them to the segmentation process. The real challenge for me here was familiarizing myself with the Image Processing Design kit by Dr. Harald Galda.
Erosion and Dilation
So remember how operations with sets work? Morphological operations rely on those. In this part of the activity, we focus on two simple operations: Erosion and Dilation.
Aside from our image, we need a structuring element to operate with. This can be almost any shape, as I’m about to show, and quantitatively dictates how an object changes.
I’ll get right down to showing how this works using the following base images.
Okay, so to put it simply, erosion removes pixels, while dilation adds them in. Just how many pixels are removed or added is determined by the structuring element.
Let’s start with dilation. We apply the structuring element to every pixel on the base image, noting the reference pixel of the element, and all the pixels not yet already part of the image are added in.
In erosion, we again apply the structuring element to every image on the base, again noting the reference pixel. If the applied structure element can fit into the base image, that pixel is kept. Got it? No? Maybe demonstrating it would be better.
And that’s how it works for a simple 2x2 square element. But a structuring element can be almost any shape, right? So here are more examples.
For a 2x1 vertical line:
For a 1x2 horizontal line:
For a 3x3 cross shape:
For a 2-pixel diagonal line:
Segmentation
So how can this help in segmentation? Well, in cleaning mostly. We can use morphological operations to clean up images for an easier time in segmentation.
We used the IPD kit by Dr. Harald Galda to perform morphological operations on an initial segmentation on the following image of “cells.”
That’s not all it can do, though. This kit allows us to identify and index continuous “blobs” in the image. We can then proceed to do things like measure the area of these blobs and use the histogram of these areas to find an acceptable range for an average “cell.”
So with this new range we can do things like filter out abnormal “cells.” Take this image for example. We were asked to highlight the “cancer cells,” using image segmentation.
Another life saved by Mighty Morphological Operations.
Like I’d said in the introduction, I struggled a bit when it came to using the functions of the IPD kit. I still managed to produce the required output, but again with not much exploration. However, I like to think my presentation deserves a bit of credit, so I’m giving myself 10 out of 10.
I’ve pretty much stopped attaching parts of my code, but that’s mostly because, besides the new functions of the IPD kit, it would mostly appear similar to at least one of the snippets I’d previously uploaded.
Time for some Mighty Morphin’ Acknowledgements:
Dr. Harald Galda, the legend himself who provided the IPD kit
Dr. Maricor Soriano, for guiding us in using the IPD kit to save the life of the cancerous paper
Joshua Abuel, my classmate who helped me realize the power of the IPD kit
This activity focused on using colour values to pick out objects in an image. I’ve always wondered how image segmentation worked with coloured objects, and this activity gave me a lot of insight in the process. I quite enjoyed bringing out the inner demon of my chosen image. You’ll see what I mean soon enough.
Okay, to begin, we had to familiarize myself with normalized chromaticity coordinates (NCC), which were used to represent the colour values in our images.
Basically, we can take the colour of a single pixel and represent it in values of Red, Green, and Blue. If we choose our coordinate system such that the sum of all three values representing red, green, and blue must be equal to 1, we can instead represent this colour using only two coordinates; chosen here to be r and g.
With that laid out, let me introduce you all to my test image: Happy Rabbit.
It’s so cute.
So we had two methods of segmentation to cover: parametric, and non-parametric.
Parametric Segmentation
This method relies on a Probability Distribution Function (PDF) generated from a selected subregion within the region of interest (ROI) we want highlighted. Here’s what I got from a particularly orange part of Happy Rabbit.
So this should give an idea of whether or not a pixel can be considered as part of your ROI depending on its r and g values.
With this, we found the probability for each pixel in the Happy Rabbit image to be part of the ROI. The result is... frightening.
The demon is rising. So to help Demon Rabbit enter our world, we set a probability threshold for which points should be considered in our ROI. Naturally, decreasing the threshold allows more points to be considered, but these points might not be part of Happy Rabbit, but instead part of the soft fabric it rests its happy body on. It’s up to us to find a good threshold.
Non-Parametric Segmentation
This method relies on histogram backprojection. Basically, we take the NCC histogram for our subregion, and use this to pick out pixels that should be included in our ROI. It’s kind of like having a a phone book to look up numbers you wanna call.
And with this histogram, we got Demon Rabbit v2.
No thresholds here. The segmentation is entirely dependent on the subregion chosen, and because of that, it can include a wider variety of colours in the segmentation, shown here with most of the shadow and face now being included, in contrast to our previous Demon Rabbit.
Now wasn’t that fun, summoning forth the Demon Rabbit from the abyss? I think it was. And hey, I also accomplished the requirements for this activity. Not bad, if I do say so myself. I’ll give myself a 10 out of 10 here for that. I didn’t really explore far beyond what was asked, so I don’t think I deserve any bonus points in this one. All in all, I had a good time.
Long live Happy Rabbit.
Here are my Acknowledgements and References:
Wondercat, the creator of Happy Rabbit
Anthony Fox, my classmate and friend who explained to me about histogram backprojection
I was excited to be working with the Fourier transform again, so I enjoyed this activity. The difficult part for me was figuring out how to do the masking thing. I was caught in between trying to automating the process of making masks or just drawing them myself in Photoshop, which I used in this activity with scilab 5.5.2.
Anamorphism in Fourier Space
Fourier space demonstrates anamorphism, which is a property that makes it so something wide in one axis in regular space becomes thin in Fourier space, and vice versa. To observe this easily, we used a pair of rectangular apertures.
They look squished, as expected. But what happens if I try it with two dots? In the last activity, we showed the FT of two dots on the x-axis is a corrugated roof, based on a sinusoid.
Looks like spreading the dots further apart increases the frequency of the sinusoid. This makes sense, since these dots would represent the frequencies of the sinusoid, and spreading them further apart would have them correspond to a greater displacement from the origin, hence a greater frequency.
Rotation Property
“You spin me right round, baby, right round. Like a record, baby, right round, round, round.”
-Dead or Alive, 1984
Using these convenient sinusoids, we can demonstrate the rotation property of the Fourier transform.
Look at them go. The transforms essentially follow the rotation of the original signals.
While we’re rotating sinusoids, we can show how the FT handles superposition.
Now that’s interesting. Here it is again with one of the sinusoids having a different frequency.
When the signals are added, their transforms seem to be mirrored somehow. Let’s take this a step further and add some rotated sinusoids and see what happens.
Crazy, but still following some form of symmetry.
Convolution 2: Electric Boogaloo
So convolution was fun in the last activity. Time to do more.
First, we get that our slightly overused two-dot pattern produces a sinusoid, but what happens if we replace the dots with circles? How about Gaussians?
As you can see, the result shows the corrugated roof, but it’s constrained to the pattern you’d expect to see when you take the FT of a single circle or Gaussian, as shown in the previous activity.
This actually shows a part of what happens in convolution.
The two-dot pattern can be considered as two dirac delta functions, having a value of 1 at two points, and zero everywhere else. When you convolve two functions, you multiply them in Fourier space. Look at the FT’s. They’re the result of a multiplication between a sinusoid and the FT of the circle/Gaussian! This means that the original images we used are actually the result of a convolution of the two-spot pattern and a corresponding single circle or Gaussian pattern.
Now let’s prove this using the conv2 function in scilab. I’ll be convolving two dirac delta patterns with arbitrary patterns.
Aha, so I understood it right. Nice.
Mask Charge!!
With the words “transform” and “mask” being brought up, there’s no way I couldn’t reference this.
Last activity, we used convolution to do edge detection. This time, we’ll use it to do filtering.
Filtering is an important part of image processing, as it can be used to remove or enhance patterns. Here, we’re gonna clean up a picture of a fingerprint, just like how the Power Rangers cleaned Angel Grove of the evil of Rita Repulsa’s monsters.
Like the Power Rangers, we need masks to do the job. Our masks, however, don’t hide our faces while we perform martial arts moves, as much as I want them to. Instead, they hide unwanted frequencies in Fourier space.
In this case, I use a filtering mask in Fourier space to eliminate those nasty frequencies I don’t want showing up in my image so that I can get a clearer look at this fingerprint, generously contributed by my classmate, Roland.
After the filtering, it was much easier to binarize the image, so we can see all those lines and whorls in their curvaceous glory.
Mask Change I
Since we literally started small with fingerprints, let’s think big. Like really big. Like moon big.
Look at this picture of the moon.
How are the Power Rangers supposed to find Rita Repulsa’s moon palace with all these horizontal lines messing up the view? Well, we can help them out with our filtering powers!
I knew which frequencies to filter out by applying that anamorphic property we learned about earlier. And with this, we can get back to our fight against evil finishing this activity.
Mask Change II
Here’s something that I’d find challenging to sneak a Power Rangers reference into.
We can use filtering masks to remove the canvas from this image of a painting. This can be used to study the brush strokes or the pigmentation used by the artist. Not only that, but by studying the canvas pattern itself, we can determine what kind of canvas the painting was made on, which can help when dealing with art fraud and the like. Then we can toss them into a space dumpster where they will be trapped for 10,000 years. Ha, I did it.
While we’re disposing of evil people that seek to profit from lies and deceit, we should take another look at that filtering mask I made using Photoshop. I’ll admit it, I couldn’t automate this one.
The resulting pattern looks kinda like the canvas, which shouldn’t be surprising, since that’s the very pattern we were filtering out.
I know this time I didn’t do a lot of exploring out of the scope of the activity, but I was pressed for time. Maybe if I were a Power Ranger...
Anyway, I hope I did a good job explaining what the activity was about. I didn’t put any codes in this time, as they were similar if not identical to those used in the last one. I’ll give myself a 9 out of 10 for this one. I satisfied the requirements, but I think I could’ve done more. Again, maybe if I were a Power Ranger...
I’d like to acknowledge these people as my Power Rangers, arriving just in time to save me from not finishing my codes:
Joshua Abuel, for helping me eliminate the distracting noise and enhance the power of the fingerprint pattern
Louie Rubio, for guiding me in figuring out which frequencies to block in order to scan the moon for Rita Repulsa’s moon palace
Roland Romero, for lending me a hand, or finger, in the fight against evil
Dr. Maricor Soriano, for allowing me to post in this format in which I can reference the Power Rangers
I really like the Power Rangers. I hope the movie next year will be good.
Ah, the Fourier Transform. This is a topic that’s come up often, as it’s used a lot in signal processing. I knew scilab had a function for the Fourier Transform, so I figured this activity wouldn’t be too bad. The tricky part was dealing with the image conversion to matrices, and vice-versa. I used scilab 5.5.2 and Paint for this activity.
Discrete Fourier Transform
So what is the Fourier Transform? Well, it’s a linear transform that converts a signal into frequency or Fourier space. For example, applying the Fourier Transform to a signal dependent on time t, it’ll express the signal in terms of frequency 1/t. This holds true for space, too, as occurs in this activity.
The Discrete Fourier Transform or DFT is used when we have a discrete signal, like a matrix of pixel values that make up a bitmap image.
Here’s the code I used to apply the DFT to any monochrome bitmap images I had. FFT here stands for Fast Fourier Transform.
//reads the image
I = imread('5a-img.bmp');
//takes fft of the 2-d image
FI = fft2(double(I));
//plots fft of image
imwrite(uint8(255*abs(FI)/max(abs(FI))), "5a-img-fft.bmp");
//abs computes the modulo
//fft matrix is normalized by dividing by maximum
//255 adjusts values for plotting
//plots shifted fft of image
imwrite(uint8(fftshift(255*abs(FI)/max(abs(FI)))), "5a-img-fftshift.bmp");
//fftshift restores the quadrants to their original places
//plots the fft of fft of image
FI2 = fft2(FI);
imwrite(uint8(255*abs(FI2)/max(abs(FI2))), "5a-img-fftwice.bmp");
Now I have a bunch of shapes to try this on thanks to the last activities, but the most basic one to try it on is the circle.
Take a look at that, a nice pattern for sure. Note the difference between the unshifted and shifted FFT patterns. This shows what the fftshift function does; it restores the quadrants of the image to their original places after they were moved around by the transform. Also, taking the FFT of the FFT seems to restore the image to it’s original state... or does it?
Will you look at that, it actually rotates the image upside down. It just isn’t noticeable on the radially symmetrical circle. Again, this has to do with the shifting of the quadrants of a matrix when the fft2 function is applied to it.
Now how about we apply it to a corrugated roof?
Huh, that doesn’t seem right, does it? The corrugated roof is based on a sinusoid function, which should have only two frequencies. That means the shifted FFT image should have only two white points, not... five.
What’s going on here? Well it turns out that saving a matrix as a bitmap image eliminates all negative values, turning them into zero. That’s the reason why the corrugated roof image looks funny, and why more than one point shows up in the FFT image. So before saving the image, I can introduce a bias to the sinusoid so that its entire shape is saved. That should solve the problem... right?
Kinda. There’s a third point right there at the origin. That’s what happens when you take the FT of a constant, like the bias I just added. Alas, we are imperfect beings. At least this time I know the peak at the origin can be discarded, giving a result that coincides with theoretical expectations.
Here is the DFT applied to other patterns.
It’s worth mentioning that the FT of a Gaussian bell curve is another Gaussian bell curve, lacking in what is called the Airy pattern most notable in the FT of a circle and square. This admittedly cool property is why the Gaussian filter is popular.
The FFT of mashpotatos dog is underwhelming, isn’t it? It’s similar to what happens if you take the FT of a circle relatively large in comparison to its background, which is basically what mashpotatos dog is (aside from being soft and warm like mashpotatos).
Convolution and Edge Detection
Keeping it simple here, when two signals are convolved, it results in a signal similar in appearance to both. The process of convolution used in this activity involves multiplying the images in Fourier space.
As an example, convolution can be used to simulate how aperture size affects image quality. Here’s the code I used.
//reads the images
R1 = imread('5b-cir1.bmp'); //0.3 r circle aperture
R2 = imread('5b-cir2.bmp'); //0.5 r circle aperture
R3 = imread('5b-cir3.bmp'); //0.8 r circle aperture
a = imread('5b-vip.bmp');
Here, a larger aperture results in a sharper image, which simulates the effect of diffraction in a camera. It doesn’t mean that larger is always better, though. If the aperture gets too big, aberrations can occur. Learned that one from Physics 165.
So what else can we do with convolution? Well, edge detection is a pretty useful application, as shown in earlier activities, and it can be accomplished by convolving an image with an edge pattern, as seen below.
//reads the image
a = double(imread('5b-vip.bmp'));
//plots the convolved images
imwrite(mat2gray(Conv1), "5d-Edge1.bmp");
imwrite(mat2gray(Conv2), "5d-Edge2.bmp");
imwrite(mat2gray(Conv3), "5d-Edge3.bmp");
That’s pretty amazing. Naturally, I had to try it with mashpotatos dog and the spot pattern.
Oh dear.
Correlation and Template Matching
Correlation measures the similarity between two signals, or images in this case, and is something I’m really familiar with at this point. In a process similar to that of convolution, one image is multiplied to the complex conjugate of the other in Fourier space.
Since correlation measures similarity, we can use it to find areas in which one image is similar to another. This method is called Template Matching, and here’s an example in which I find the letter ‘A’ in a familiar saying.
//reads the images
R = imread('5c-rain.bmp');
A = imread('5c-a.bmp');
//takes fft2 of both images
Fr = fft2(double(R));
Fa = fft2(double(A));
//correlates the images by multiplying one by the conjugate of the other
FCorr = Fa.*conj(Fr);
//plots the results of the correlation
Corr = fftshift(fft2(FCorr));
imwrite(uint8(255*abs(Corr)/max(abs(Corr))), '5c-Corr.bmp');
The “peaks” of correlation are found at the whitest spots, and it should be no surprise that they correspond to the locations of the letter ‘A’ in the original image.
Now that I know it works, how about I apply it to a game where the player’s supposed to find an object in an image, like in Mystery Case Files! Would that be considered cheating? Time to find out.
This particular screenshot is from Ravenhearst Mansion. Ah, fond memories of playing this game... Anyway, let’s look at that list. Hey, a scorpion! I know where the scorpion is!
Guess it can’t be cheating if I have to find the object to apply the process. Okay, here’s what happens when I apply Template Matching to find the scorpion.
Don’t see it? I don’t blame you. Here, let me point out the peak. It should be over the scorpion in the original screenshot.
Not very effective here, is it? Well, that’s because Template Matching compares pixel values, and if there are a lot of pixels that end up matching values, correlation yields a high value. When converting the screenshot into a monochrome bitmap, the bright-coloured sky and the bright-coloured wall end up with similar values, hence the high correlation values. Unfortunately, this is an inherent limitation in the Template Matching method, so I can’t just do away with it.
At least I know I can’t cheat in Mystery Case Files using correlation.
And yes, I tried it with mashpotatos dog too. Naturally.
Science has gone too far.
This was a fun activity, though writing the blog took time. I owe a lot to our professor Ma’am Soriano and my classmate Roland for help in troubleshooting my code.
I’m giving myself a 12 out of 12 on the count of mashpotatos dog. Seriously though, I think I put a lot more into this activity than the previous ones. I believe I did a good job explaining what was happening in the activity through the text and figures, and I’d applied what I learned in previous activities as well.
I regret not doing edge detection with the mashpotatos dog sooner.
And here are my acknowledgements and references:
Dr. Maricor Soriano, for explaining the process of FT, convolution, and correlation in regards to image processing, and for discussing about the limitations of the Template Matching method
Dr. Percival Almoro, for teaching me about the optics stuff in cameras and other imaging devices
Scilab Online Help, for references on syntax and functions
Roland Romero, for helping me with plotting the result of the DFT
This one’s a bit heavier on the programming side than the others before. I’ll admit I was a little intimidated when I started reading the instructions, especially when Green’s Theorem was mentioned. But when I actually got started, it turned out to not be as bad as I thought. I still struggled a bit with the syntax and all of scilab, but in the end, I pushed through. This activity required the SIVP module in scilab 5.5.2, and ImageJ.
Green’s Theorem
So this activity was all about areas. And a way of computing for the area of a two-dimensional geometric shape is through Green’s Theorem. Here’s the form of the Theorem used in the construction of the code.
This form is necessary, as the area can be expressed as an integral, and integrals can only be solved for numerically through methods such as summation.
Alright, instead of waiting any further, I’ll show the code I used to calculate the area of a 2-D shape now. This code only works with monochromatic bitmap images, as is the limitation of the imread function used here.
//reading image
Img = imread("image.bmp");
E = edge(Img, 'prewitt'); //finding edge of shape
//sorting edge coordinates
[xe,ye] = find(E);
xef = xe - xcent; //fixing with respect to centroid
yef = ye - ycent;
r = sqrt(xef.^2 + yef.^2); //finding distance r
th = atan(yef, xef); //finding angle theta
coord = [];
for i = 1:length(xef); //compiling a list for sorting
elem = [xef(i) yef(i) r(i) th(i)];
coord = cat(1,coord,elem);
end
[sort,j] = gsort(coord(:,4),'lr','i'); //j = index wrt theta
coord = coord(j,:); //sorting wrt j
//computing area using Green's theorem
summ = 0
for i = 1:length(coord)/4;
if i == length(coord)/4;
summ = summ + (coord(i,1)*coord(1,2)-coord(i,2)*coord(1,1));
else;
summ = summ + (coord(i,1)*coord(i+1,2)-coord(i,2)*coord(i+1,1));
end
end
Area = 0.5*summ;
The comments should help in explaining how the code works, but I’ll give a quick rundown too.
After reading the image, edge is used to find points on the edge of the shapes. The second argument in the edge function interested me, since when I used the default values to find the edge of a square, it only gave the four corners. This argument turned out to be the edge detection method, with the default being sobel. Mich, a classmate and labmate, told me to try using prewitt, and that gave me something that looked more like a square’s edge.
Still, I wanted to see which of the edge detection methods worked best. To do this, I tested prewitt, canny, and the default sobel. There are two more methods available to scilab: fftderiv and LoG, but fftderiv wouldn’t work, and I didn’t want to use LoG just yet.
The centroid is found by taking the mean of the x and y positions of all the points on the shape. This is set as the origin when computing for the angle θ of a point on the shape’s edge. This angle is used to sort the edge points before summation. To be honest, figuring out this sorting part is what took up most of my time on this project.
Next, Green’s Theorem is applied, and we have the area. Yay.
Area of a Square and a Circle
For this activity, I tested the code first on two simple shapes: a circle and a square. I did the test for two resolutions, just to get an idea of how much the resolution affected accuracy.
Now you see the four points sobel returned for the edge of the square. Looking at the circle, sobel seems to have produced a series of disjointed points, for the most part. Prewitt edge detection gave a thicker edge, meaning more points. Finally, canny gave a thin edge, note the lack of the corner pixels for the square edge though. How did this affect the accuracy? Well,the results were pretty interesting.
For the square, prewitt and canny beat sobel, and it’s not hard to see why, with sobel only giving four points for the edge. It looks like the two are matched though, not being entirely perfect themselves.
And for the circle, what’s this?! Sobel wins! Prewitt comes in second, and canny pulls a third. That’s unexpected, but I guess it has to do with how the Green’s Theorem allows us to estimate the area. The world is full of surprises.
Area of... the Sunken Garden?
Alright, the code can compute areas for simple shapes, but what about the area of a place? How about a landmark like the sunken garden in U.P. Diliman? Let’s put it to the test, with the three edge detection methods. Here I used Google Maps, with the measurement feature in My Maps to get my theoretical area. The scale Google Maps provided at the zoom I used was 2.8 pixels/m.
Looks like all the expected traits present in the simple shape edge detection show up here as well. That’s good. Yay for consistency.
And prewitt wins this round! Of course, 8% isn’t exactly as good as the previous results with the simple shapes, but I think this can be improved on by working with larger resolutions, as the test with simple shape shows. The resolution of the image I worked with was 1076p x 705p, which is as big as I can get the whole sunken garden in Google Maps to appear on my screen.
Measurements in ImageJ
I’ll admit, I’m a Adobe Photoshop or Paint Tool Sai guy. I’m not really familiar with ImageJ, so I didn’t expect this feature. It turns out you can use ImageJ to analyze the size of objects present in the image. Granted, ImageJ was created with the idea of analyzing biological samples in mind, but hey, I did say I wasn’t familiar with it, right?
Basically, you can set the ratio of pixels to whatever unit you choose, and analyze the length of any segment you draw, or the area of any shape. That’s pretty cool, so I wanted to try it on something standardized, like a poker-size playing card.
This is the back of a card in the Bicycle Zombified deck, featuring art by Billy Tackett. I think it fits in the theme of biological specimens. Anyway, poker-size cards have a standard length of 88.9 mm, and a standard width of 63.5 mm, giving a rectangular area of 5,645.15 mm^2.
I set the length of the card into ImageJ after drawing a line segment along the card’s length. It returned a scale of 11.694 pixels/mm. Then I drew a rectangle over the card to get its area. I chose to ignore the curved bits at the corners this time, so I could get a clean rectangle. The area ImageJ gave me was 5,679.585 mm^2. This gives an error of 0.61%. Not bad at all.
Okay, that was fun. I owe a lot of the ideas behind the sorting part of the code to Roland, another classmate and friend, and once that was done, the rest just kinda flowed.
Alright, time to evaluate myself. Okay, the output is all there, and the figures and tables have captions at last! Good job, me.
The tables and images are all transparent now, and that is a plus in my own book. Finally put all the Photoshop know-how I’ve amassed across the years into making my tables transparent.
The ratios from Activity 2 and plotting from Activity 3 came in handy for this too. I’m going with 11 out of 10 this time. Earn those bonus points, me.
Finally, for my acknowledgements and references:
Scilab Online Help, for references on syntax and functions
Google Maps, for the map of the U.P. Diliman Sunken Garden
Michelle Cirunay, for enlightening me on the different edge detection methods
Roland Romero, for helping me solve the riddle of sorting edge coordinates
Okay, I admit to putting this activity off until the weekend after we came back from the Samahang Pisika ng Pilipinas conference in my hometown, Iloilo City, but only because I wasn’t feeling very well. Eventually, I sat down and got to it, thanks to the help of Albert, a classmate and friend. I used the Kingsoft Spreadsheets software, along with Paint and Photoshop for this activity.
The original image
The goal of this activity was to recreate a plot from an old journal. The plot had to either be hand drawn or created with an old XY plotter, which means we had to dig through the College of Science library for an old publication. I didn’t want to have to do that before the SPP conference, so I was planning on waiting until we came back and looking through the older theses in our lab’s archive.
Fortunately, Albert, who was a labmate as well as a classmate, offered to help find me an old plot while looking for his own. Thus, from the Canadian Journal of Physics Volume 29 published in January 1951, comes this plot.
It’s from a paper titled “The Scattering Lengths of the Deuteron” by D. Hurst, et. al. A deuteron is an isotope of hydrogen, consisting of a proton and a neutron, and would be the nucleus of the atom called deuterium. What makes a deuteron special is its stability, as a free neutron would decay, but a deuteron keeps its neutron through its binding energy.
The plot itself was drawn using an old XY plotter, and corrected for rotation and warping using Photoshop.
Recreating the plot
This is where most of the elbow grease went into. Our professor provided the outline of the method: open the plot in Paint, use the X and Y pixel values to determine the ratio of pixels to ticks on the plot, get points from the image and plot them, and compare the result to the original image.
First I went about finding that pixel to tick mark ratio using Paint. I recorded the locations of the ticks then took the average distance between them to serve as my ratio. I got 74.5 pixels per 0.1 units for the X-axis, and about 74.4 pixels per unit for the Y-axis. I used Kingsoft Spreadsheets because I don’t own a copy of MS Office myself, nor did I want to use Google Sheets for fear of connection issues taking a toll on my admittedly short patience.
That wasn’t too bad. Next came the actual points. Looking at the plot, there are six points where measurements were taken, corresponding to scattering angles. I took their pixel locations, not too bad since the points are aligned on the Y-axis, and also recorded the Y-intercepts for the three lines formed, just to ensure I got lines that reached the Y-axis.
As you can see above, I did corrections for my Y pixel values since Paint counts starting from the top of the canvas. I also corrected for the margins using the location of the origin. Only then did I finally make the plot. You can see the resulting line graph, but I used the scatterplot function in Spreadsheets, with connecting lines, since the values I had were technically ordered pairs of X and Y. I manually set my values for the axes, but the ticks were set using the ratios.
Comparing the two images
Again, I used Spreadsheets, which doesn’t have the function that lets you use an image as a background to a graph, as far as I know, so I copied both images onto Photoshop and set the layer containing the original plot to Multiply, eliminating the white of the paper and allowing a comparison to be made.
Overall, I’d say they line up quite well, don’t you?
Well, that wasn’t nearly as bad as I thought it’d be! That’s probably because I used only seven points for each line, not to mention the fact that they’re lines. Three lines, yes, but lines, nonetheless. I’d be shaming myself if I wasn’t able to recreate something like this. Still, my major concern was finding an actual plot to use, and with Albert’s help, that wasn’t a problem. I will forever be grateful to you, Albert.
Now it’s time for some self-evaluation. Okay, my output is there, and I did show the overlay comparison, although that was with the help of Photoshop. Speaking of Photoshop, I was able to use pixel ratios to fix the warping of the original plot, similar to what we learned in our first meeting of this class. That’s something, right? And I was able to recreate the plot well, even though it just consisted of three lines.
But, oops, no captions again. I’m starting to think I won’t use them at all, since Tumblr doesn’t really have a good format for those, unless I switch to using LaTeX, which I probably won’t be doing anytime soon.
Overall, I’d give myself a 9 out of 10. I do feel like I didn’t put in enough effort as I could have, with me not finding my own image and all, and it being an easy one to match.
Here come the acknowledgements and references,
D. Hurst, et. al., “The Scattering Lengths of the Deuteron,” Canadian Journal of Physics Volume 29 January 1951
HyperPhysics, for telling me what a deuteron is
My classmate and friend, Albert Yumol, for finding a plot for me to use
I gotta admit, I was pretty excited to do this activity, since I didn’t have a lot of experience with scilab in general, and I wanted to learn how to at least plot using it. I also chose to do this first because I wasn’t ready to get up and look for a journal just yet. For now, I wanted to sit and code. For reference, I coded this on scilab 5.5.2.
Circular Aperture
Alright, so our professor gave us an example that plotted a circular aperture. I used this as the basis for pretty much every other code I wrote here. The code and the plot produced is as follows:
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
r= sqrt(X.^2 + Y.^2); //note element-per-element squaring of X and Y
A = zeros (nx,ny);
A (find(r<0.7) ) = 1;
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Square Aperture
“Alright, cool. It’s easy to understand,” I told myself as I started on the square aperture. I immediately thought to use an absolute value function to define the shape. A good idea, but in my zeal, I didn’t notice I was using two absolute value functions instead of using an and operator, and long story short, I made a big cross shape. At least then I had an idea for the cross shape later on. Anyway, here’s my code and the plot it produced:
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
rx= abs(X); //absolute value function
ry= abs(Y);
A = zeros (nx,ny);
A (find(rx<0.6 & ry<0.6) ) = 1; //finds the perimeter of the square
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Corrugated Roof
Progress. Next was the corrugated roof. It needed a sine function, but I wanted to keep the plot to be from -1 to 1. Now if you know the sine function, you know it has a period of 2π. If I plotted the sine function from -1 to 1, it’d be pretty boring, wouldn’t it? The answer is yes.
Okay, so the solution I came up with was to map the sine function over something like -6π to +6π onto the plot we had now. So now the axes match, and my somewhat unnecessary problem is solved. Just a warning, it’s a little tiring on the eyes. Here’s what I did:
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
x1 = linspace(-7*%pi,7*%pi,nx); //range for the sinusoid
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
rx= cos(x1); //sinusoid function
A = zeros (nx,ny);
for j = 1:nx //inserts the sinusoid function into the array
A (j,:) = rx(j);
end
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Grating
Okay, for the grating, I wanted to take my sine result and apply an if condition that replaces all positive values with 1, and all negatives with 0. It worked and I was happy with my achievement. Especially after thinking very hard about the for loops.
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
x1 = linspace(-7*%pi,7*%pi,nx); //range for the sinusoid
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
rx= cos(x1); //sinusoid function
for i = 1:nx; //converts the sinusoid into a grate function
if rx(i)<0
rx(i)=0
else
if rx(i)>0
rx(i)=1
end
end
end
A = zeros (nx,ny);
for j = 1:nx //inserts the sinusoid function into the array
A (j,:) = rx(j);
end
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Annulus
Up next was the annulus, a scientific way of saying a ring shape. Not much explaining here, I just added another line of code defining a smaller circle. Ta-dah.
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
r= sqrt(X.^2 + Y.^2); //equation of a circle
A = zeros (nx,ny);
A (find(r<0.7) ) = 1; //forms the outer circumference
A (find(r<0.4) ) = 0; //forms the inner circumference
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Circular Aperture with Graded Transparency
This next one, I skipped at first. But I eventually came back and did it. The graded transparency is added to the circular aperture using a Gaussian function. This function was defined first, then applied to that A array I’ve been using. I then cut the rest of it off with the circular aperture from the sample code. The result looks pretty cool.
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
function y = gauss(x,s) //defines the gaussian function to be used
y = exp(-x^2/(2*s^2))
endfunction
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
r= sqrt(X.^2 + Y.^2); //note element-per-element squaring of X and Y
A = zeros (nx,ny);
A = gauss(r,8); //applies the gaussian to the A array
A (find(r>0.7) ) = 0; //adds the circular limit
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Ellipse
I’m not gonna lie, the ellipse was the hardest for me. Not because of anything with the equation, mind you. I know the equation of an ellipse. I could have sworn it was just refusing to work, no matter what I tweaked. I even asked my classmate and friend, Josh, about it and we were both stumped for a while... Then we saw it.
I had put the periods inside the parentheses beside X and Y like in the example code. I forgot to move them out. Correcting that, it produced an ellipse and I could finally breathe a sigh of relief.
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
x1 = linspace(-5,5,nx); //defines the range
y1 = linspace(-5,5,ny);
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
r= sqrt((X/0.7).^2 + (Y/0.5).^2); //equation of an ellipse
A = zeros (nx,ny);
A (find(r<1) ) = 1;
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Cross
And last but not least is the cross. No problems here, I applied what I’d learned earlier from making the square to get the horizontal and vertical parts of the cross. Then I defined a perimeter using the or operator that clipped the cross off where I wanted it.
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
rx= abs(X); //absolute value function
ry= abs(Y);
A = zeros (nx,ny);
A (find(rx<0.3 ) ) = 1; //creates the horizontal part of the cross
A (find(ry<0.3 ) ) = 1; //creates the vertical part
A (find(rx>0.6 | ry>0.6) ) = 0; //limits the cross
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Exploration
Alright, then came the fun part. We got to combine a bunch of images to form something cool. I already knew what I wanted to form, so I took it step by step, masking the images on like they were layers in an image editing software. It was nice getting to combine the and and or operators to define boundaries, particularly in the cross sections within the iris/pupil. Haven’t played around with logic like that since the philosophy GE courses.
That’s really all I gotta say on the method. So behold, an eye, definitely not representing how awake I am, writing this at 2am. Rather, it’s what came to mind when I thought of the word “exploration.” To explore is to see, or witness, in a sense.
stacksize(13000000); //changes the stacksize to allow for more memory
nx = 1000; ny = 1000; //defines the number of elements along x and y
x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y); //creates two 2-D arrays of x and y coordinates
riris= sqrt(X.^2 + Y.^2); //origin-centered circle
rtop = sqrt(X.^2 + (Y+0.4).^2); //circle centered at y=+0.4
rbot = sqrt(X.^2 + (Y-0.4).^2); //circle centered at y=-0.4
rx= abs(X); //absolute value function
ry= abs(Y);
A = ones (nx,ny);
A (find(riris<0.45) ) = 0.5; //draws the "iris" of the eye
A (find(riris<0.35) ) = 0.3;
A (find(riris<0.3) ) = 0.1;
A (find(rx<0.075 & riris<0.3) ) = 0.15; //creates the x part of the cross
A (find(ry<0.075 & riris<0.3) ) = 0.15; //creates the y part of the cross
A (find(rx<0.05 & riris<0.3) ) = 0.5; //creates the x part of the inner cross
A (find(ry<0.05 & riris<0.3) ) = 0.5; //creates the y part of the inner cross
A (find((rx>0.2 | ry>0.2) & ((rx<0.075 & riris<0.3) | (ry<0.075 & riris<0.3))) ) = 0.15; //limits the inner cross
A (find(rtop>1.) ) = 0; //cuts off bottom part of the eye
A (find(rbot>1.) ) = 0; //cuts off top part of the eye
f = scf();
grayplot(x,y,A);
f.color_map = graycolormap(32);
Overall, I’d say I really did have fun making all these. I learned quite a bit, and I liked looking at these plots, especially when I figured out how to set the axes to have the same scale. Look at them, perfect squares... Anyway, this was a good activity that I could do at a nice pace, especially with the SPP conference happening around the time. I’m thankful for being able to do it.
Now comes the part where I evaluate myself. Hmm, all the required outputs are there. And I took great pains to ensure they were of uniform size, but oh man, my captions could probably be better. Also, I was able to use past lessons on both programming and logic in this so hey. I’d give myself a solid 10 out of 10, with initiative making up for the overall lack of captioning.
As far as my acknowledgements and references go,
Scilab Online Help, for references on the syntax and such
My classmate and good friend, Joshua Abuel, for helping me with the dreaded ellipse
