Tuesday, June 22, 2010

Scilab Basics... or not

According to the manual prepared by Dr. Maricor Soriano, Scilab is a high level programming language which is an acceptable substitute for Matlab but most importantly, it is FREE. Although I can't really compare its functions with Matlab (because I haven't used Matlab yet), I believe the programming language is pretty cool and neat primarily because it treats variables as arrays and matrices.

The first task was to create a centered circle. This exercise allowed us to basically test some of the features of Scilab and in doing so, we were able to create more complex images.
Here is the code on how to create a centered circle.

code 1:
nx = 200; ny = 200; //defines the number of elements along x and y
x = linspace(-100,100,nx); //defines the range
y = linspace(-100,100,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); //create an array of zeros
A (find(r<70)>//circles with radius less than 70 is given a value of 1
imshow (A, []); //the image is shown

Note that there are characters in green, they are comments detailing what that particular line of code does for the program. Figure 1 shows the image produced by running the code 1.

Figure 1: Synthetic Centered Circle. A white center circle indicates a value of 1 and the black boundary is 0.

From code 1, increasing the number of elements, nx and ny, would be very conventional because the greater the number of elements, the greater is the image's resolution but alas, the limitation of Scilab comes into play. Scilab has a limited total memory available, meaning, Scilab cannot handle any high resolution images. I tried setting the value of nx and ny to 1000 but a message popped out detailing an error about total memory availability.

With these basic informations about Scilab and its features, I was able to create more elaborate images. First up, a centered square aperture with its source code given by code 2.

code 2:
nx = 200; ny = 200;
x = linspace(-100,100,nx);
y = linspace(-100,100,ny);
[X,Y] = ndgrid((x/100),(y/100)); //scaling; the value would still be from -1 to 1 but with a bigger size

A = zeros(nx,ny);
A(find(abs(X)<0.5>//joint conditions, any point in the array that satisfies BOTH condition would be given a value of 1

imshow(A, []);

Figure 2: Centered Square Aperture.

New about code 2 is the joint condition, wherein both conditions must be satisfied in order for the action/consequence to take place. In this case, abs(X) must be less than 0.5 AND abs(Y) must also be less than 0.5 for that particular point to be given a value of 1. Note that X and Y are the coordinates/pixel location of a particular point in the array.

code 3:
nx = 200; ny = 200;
x = linspace(-100,100,nx);
y = linspace(-100,100,ny);
[X,Y] = ndgrid((x/100),(y/100));

func = (1+sin(2*%pi*Y.*5))/2; //produces a sinusoidal function with a period of 5
imshow(func, []);

Figure 3: Sinusoid along the x direction. Looks like a corrugated roof from this view (top view).

Next is code 3 and as figure 3 shows, the code produces a sinusoid along the x direction. For me, the important thing about code 3 is it allowed me to examine how a simple sinusoid function is used and projected in scilab. Arvin Mabilangan pointed out that it is necessary to shift and scale the sinusoidal function so as not to obtain negative values which could become troublesome when handled.

code 4:
nx = 200; ny = 200;
x = linspace(-100,100,nx);
y = linspace(-100,100,ny);
[X,Y] = ndgrid((x/100),(y/100));

func2 = (1+sin(2*%pi*Y*5))/2; //sinusoid function is scaled so as not to produce a negative value

func2(find(func2<0.5))>//if the value of func2 is less than 0.5, instead of using its true value, it is assigned a value of 0

func2(find(func2>0.5)) = 1; //if the value of func2 is greater than 0.5, instead of using its true value, it is assigned a value of 1.

imshow(func2, []);

Figure 4: Grating along the x direction.

In order to understand how code 4 produces figure 4, notice first the resemblance between figure 3 and figure 4. I noticed that figure 4 is basically the same as figure 3 without the gradual change in gray scale (indicating magnitude/value). Figure 3 was produced using a sinusoid function, thus it would be intuitive that Figure 4 may also be produced using a sinusoid function with a minor manipulation. As code 4 shows, a boundary value(0.5) was used to separate the white stripes (value of 1) with the dark ones (value of 0). If the sinusoid function produces a value greater than 0.5, instead of using its true value, it is replaced with 1 (white). And if the sinusoid function produces a value less than 0.5, that value is replaced with 0 (dark). With these 2 conditions governing the program, figure 4 is produced.

code 5:
nx = 200; ny = 200;
x = linspace(-100,100,nx);
y = linspace(-100,100,ny);
[X,Y] = ndgrid((x/100),(y/100));

r = sqrt(X.^2 + Y.^2);
circle = zeros(nx,ny);
circle(find(r>=40 & r<=70)) = 1; //only the points within and including the radius 40 and 70 are given a value of 1. Another case of joint condition

imshow(circle, []);

Figure 5: Annulus. Only the points within a certain radius are given a value of 1.

Figure 5 and code 5 has a lot of similarities with figure 1 and code 1, the only difference is the additional condition presented in code 5 to produce figure 5. In code 1, the only condition was, any point within the boundary of the circle with radius 70 is given a value of 1, in code 5, an additional condition (joint condition) is stated, the point should also be outside the boundary of the concentric circle with radius 40. With this additional condition, an annulus is produced.

code 6:
nx = 200; ny = 200;
x = linspace(-100,100,nx);
y = linspace(-100,100,ny);
[X,Y] = ndgrid((x/100),(y/100));

r = (X.^2 + Y.^2);
circle = sqrt(r);
aperture = exp((-1*(r)^2)/100); //Gaussian function

aperture(find(circle>0.6)) = 0; //any point outside the circle with radius 0.6 is given a value of 0

imshow(aperture, []);

Figure 6: Gaussian function used in code 6

Figure 7: Circular Aperture with graded transparency (Gaussian transparency)

Figure 7, produced by code 6, is basically figure 1, but instead of a clear, white hole as the center circle, a graded transparency (Gaussian function) is used in Figure 7. As code 6 indicates, the Gaussian function was created first, shown in figure 6, instead of the circular aperture. I did this because it was easier to do so. And then, much like what I did in the previous code, I indicated a boundary circle of radius 0.6, but this time, any point outside of this circle, I assigned a value of 0. The resulting image looks very much like a circular aperture with graded transparency(Figure 7).

This activity made me think a lot mainly because it has been a long time since I last did programming. It got me frustrated, I had to familiarize myself again with the jargons and logic involved in programming. But once I got the hang of it, I immediately enjoyed this activity. At the end of the day, I'll give my self a pat on the back by giving me a grade of 10, 5 for technical correctness (I completely understood the lesson, even though I got frustrated at first, and I was also able to produce all the necessary outputs), and another 5 for quality of presentation (given the limitation presented by scilab, I was still able to produce relatively high quality images complete with captions that clearly describes the image). 'twas a good day after all.




No comments:

Post a Comment