The Crab Decomposition

work in progress

Last January Doc Muriel discussed something about laminar fluids emitting blue light and that turbulent fluids would emit at a lower wavelength. He showed a picture of the Crab Nebula and the decomposition to blue. The result showed a blue ring at the center of the Nebula. It surprised me that such phenomena as laminar fluids in space could be seen in a simple nebula image. I never saw it again.

I'd love to have a 7-band image of the nebula; all I have so far is a snapshot from Stellarium. I tried to get the blue channel. Of course it's an RGB image and the blue channel also contributes to form other colors.

Here is an attempt to further filter the blue channel by removing pixels that overlap higher level red and green (>0.7 with max 1.0). It's a crude attempt but I'll build on this and find an actual decomposition of the nebula for comparison.

the Crab Nebula taken from Stellarium


blue channel

result of filtering out the other colors

Here are nice crab photos

x-rays from Crab Nebula

Here is the code.

A =imread('crab.png');

red = A(:,:,1);
green = A(:,:,2);
blue = A(:,:,3);

//red channel
R = A;
R(:,:,2) = 0;
R(:,:,3) = 0;

//green channel
G = A;
G(:,:,1) = 0;
G(:,:,3) = 0;

//blue channel
B = A;
B(:,:,1) = 0;
B(:,:,2) = 0;

scf(2); imshow(B)
scf(3); imshow(B(:,:,3));

//640-490nm
rc = R(:,:,1);
rc(find(rc>0.7)) = 0;

gc = G(:,:,2);
gc(find(gc>0.7)) = 0;

blue = B(:,:,3).*rc.*gc;
B(:,:,3) = blue;

scf(4);imshow(blue);
imwrite(blue, 'bluecrab.png');
imwrite(B(:,:,3), 'bluechannel.png');
imwrite(B, 'bluecrab3.png');

FOSSY Things: Histogram Equalization

Image contrast is changed by manipulating the PDF of the image. One can easily use popular image processing software like Photoshop and GIMP to achieve similar results. Inkscape, a vector-graphics software can also do the same thing. The software is a parallel of Illustrator. It is specialized for creating vectors but can also import raster graphics and do a variety of enhancements.

Here is a sample of some Inkscape tricks. The image used here was created with TuxPaint, a FOSS kiddie version of Photoshop and Paint.

Image created from TuxPaint with the Stamps extension


Here is a grayscale version of the image
Conversion was done with GIMP 2.6

You can't really manipulate the CDF curves like what we did in the previous post but it can do the same trick. With Levels, you simply change white and black points, like having a threshold for both black and white. I assume that Equalize performs a histogram equalization. For the curious, Inkscape effects are written in Python and codes can be viewed after installation.

The gray and white points can be edited with the Raster>Level feature in Inkscape 0.46


Equalize does the same trick but it doesn't have customizable settings

Using GIMP, one can manipulate the CDF by hand and see the PDF change live.

Colors>Curves trick from GIMP 2.6

Here are the final images.

Enhanced with Inkscape's Equalize


Enhanced with Inkscape's Level


Enhanced with GIMP's Curves

Inkscape, GIMP and TuxPaint are Free and Open Source and can be easily downloaded over the internet. Thanks to Leon of CPU and Last Year's Software Freedom Day for the OpenEduDisc. I stumbled upon TuxPaint on that disc in search for graphics software goodies.

Many thanks to Shugo Tokumaru and Lost in Found for making beautiful music to keep me company.

5. Enhancement by Histogram Manipulation

In the previous post on image types and formats, I used the histogram of a grayscale image to select the threshold for binary conversion. The histogram gives the distribution of pixels over various graytones or levels (PDF, Probability Distribution Function when normalized). Image enhancement can be done by manipulating the PDF of an image.

The original image with poor contrast and image enhanced by linear backprojection of the CDF

To increase the contrast of images, the pixels must be distributed over a wide range of gray levels. Histogram manipulation can be accomplished by using the CDF (Cumulative Distribution Function). The CDF is described by the following equation,

Definition of the Cumulative Distribution Function

A linear CDF would then do the trick of increasing PDF spread. Having CDF values plotted against itself would yield a linear plot. Each pixel must then take on its corresponding CDF value as its new gray level.

Notice that in the new PDF, pixels distribution are more spread especially at the 0-0.6 region. The old PDF is in 0-255 format since this is the default configuration of the imwrite() function. The new PDF is in 0-1.0 formats because it was taken from normalized values of the old CDF.

PDF and CDF of original image (left) and the enhance image (right)


linear CDF as reference for backprojection

Aside from a linear CDF, a simple matrix manipulation can create a nonlinear CDF. A logarithmic CDF (y = log(x)) was applied. CDF values do not change in the white end of the scale and more pixels are concentrated in the dark region. Consequently, the result is a low contrast, dark image.

Image enhanced by a logarithmic CDF


PDF and CDF of image with logarithmic enhancement

Since I used histplot() on the last activity, I was fixated on how to access the values from histplot(). Histplot uses dsearch() which in turn makes use of tabul(). Thanks to the excellent Scilab Help modules. A simple use of find() by Gilbert Gubatan accomplishes the same thing. The difference with the latter is that tabul() does not store values for empty classes. In the case of the featured image, the length of the array is equal to 244 instead of 256. It is necessary to consider this in backprojection. Intuitively, one would index the CDF values with 0:255, which will not work for a 244 long array. See the code below.

Special thanks to Carmen Lumban who took the photo in 2008.

I give myself a 10/10 for this activity.

Code

A =imread('e_gray.jpg');
Amax = max(A);
Amin = min(A);

//get PDF of original image
p = tabul(A,'i');
scf(0); subplot(211);
plot(p(:,1),p(:,2)/max(p(:,2)));
title('PDF');
xlabel('Grayscale levels');
ylabel('# of pixels');

//get CDF of original image
c0 = cumsum(p(:,2)); //2nd col of p
c = c0/max(c0);
scf(0); subplot(212);
plot(p(:,1),c);
title('CDF')
xlabel('Grayscale levels');
ylabel('# of pixels');

//plot uniformly distributed CDF
x = c;
scf(1);
plot(x,x);
title('linear CDF');
xlabel('Grayscale levels');
ylabel('# of pixels');


//backplot uniform CDF to image

B = []; //enhanced image matrix
i = 0;
j = 0;
k= 0;

for i = 1:size(A,1)
for j = 1:size(A,2)
k = find(p(:,1)==A(i,j));
B(i,j) = c(k);
end
end

//get new PDF
p2 = tabul(B,'i');
scf(2); subplot(211);
plot(p2(:,1),p2(:,2)/max(p2(:,2)));
title('PDF of Enhanced Image');
xlabel('Grayscale levels');
ylabel('# of pixels');


//get new CDF
c02 = cumsum(p2(:,2)); //2nd col of p
c2 = c02/max(c02);
scf(2); subplot(212);
plot2d(p2(:,1),c2);
title('CDF of Enhanced Image');
xlabel('Grayscale levels');
ylabel('# of pixels');

C = []; //enhanced image, log
i = 0;
j = 0;
k= 0;
clog0 = abs(log(c));
clog = clog0/max(clog0);


for i = 1:size(A,1)
for j = 1:size(A,2)
k = find(p(:,1)==A(i,j));
C(i,j) = clog(k);
end
end

//get PDF of new image
p3 = tabul(C,'i');
scf(3); subplot(211);
plot(p3(:,1),p3(:,2)/max(p3(:,2)));
title('PDF of Enhanced Image: Logarithmic CDF');
xlabel('Grayscale levels');
ylabel('# of pixels');


//get new CDF
c03 = cumsum(p2(:,2)); //2nd col of p
c3 = c03/max(c03);
scf(3); subplot(212);
plot2d(p3(:,1),c3);
title('CDF of Enhanced Image: Logarithmc CDF');
xlabel('Grayscale levels');
ylabel('# of pixels');

imwrite(B, 'C:\Users\regaladys\Desktop\e_gray2.jpg');
imwrite(C, 'C:\Users\regaladys\Desktop\elog_gray2.jpg');

4. Area Estimation for Images with Defined Edges

With Green's Theorem one can easily calculate the area of images. This requires defined edges and Scilab with the SIP Toolbox accomplishes edge detection with follow().

Img = imread('circle.bmp');
[x,y] = follow(Img);
plot2d(x,y, rect = [0,0,200,200]);

l = length(x);
area = 0;
area_updated = 0;
for i=1:l-1
area_add = (x(i)*y(i+1)) - (x(i+1)* y(i));
area_updated = area_updated + area_add;
end

area = area_updated/2;

plot2d() displays the edges of the image defined in matrix [x,y]. Area is calculated in terms of pixels. The analytic value of the area is the number of pixels defining the image. With a white image, each matrix element is of the value 1. The area is simply given by sum(). Accuracy of the measurement is calculated by taking the difference of calculated area and the analytic value.

area_pixel= sum(Img);
dif = area - area_pixel;
area_acc = abs(dif*100/area_pixel);

Using synthetic images generated with Scilab, I calculated their areas with Green's Theorem.

Synthetic images with calculated areas of up to 97.8% and 97.3% accuracy

With the available satellite images of Wikimapia, I was able to take images of the subdivision where we live in Bacolod City.

Satellite image of St. Benilde Homes, Mansilingan City
It appears to be taken before 2006
photo courtesy of Wikimapia

Zooming in on Block 5, the Regalado Residence is indicated by the cross marking
photo courtesy of Wikimapia

I intend to calculate the total area of the subdivision and Block 5. To create solid bw images of these photos I used GIMP 2.6. I used the free select tool to create a mask of the block's shape and filled a white layer with black.

Black and white bitmap file of Block 5 for image detection and area estimation

With the map's scale, the conversion is approximately 5.25 m/pixel. This is calculated with the method described in Activity 1.

The area was calculated to 95.6% accuracy with the block measuring up to 11,103 square meters by area estimation and 11,620 square meters as the analytic value.

Here is the full code used to calculate the Block 5 area. Note that the bitmap file is a truecolor image though it appears to be in black and white. It must be first converted to a binary image before area estimation. Without this step, follow() function will not detect images as the truecolor file will have 3 channels for RGB, not compatible for the operation of follow().

//area estimation by green's function

Img = gray_imread('block5_w2.bmp');
histplot([0 : 0.05 : 1.0], Img)
Img2 = im2bw(Img, 0.5);

[x,y] = follow(Img2);
plot2d(x,y, rect = [0,0,200,200]);

l = length(x);
area = 0;
area_updated = 0;
for i=1:l-1
area_add = (x(i)*y(i+1)) - (x(i+1)* y(i));
area_updated = area_updated + area_add;
end

area = area_updated/2;
area_pixel= sum(Img);
dif = area - area_pixel;
area_acc = 100 - abs(dif*100/area_pixel);

Thanks to Ricardo Fabri, author of most of SIP, for keeping beautiful and simple docstrings on his macros. It took me a while to realize that the .bmp file must be converted to a binary image. Nevertheless, I have completed this exercise with more ease than the previous two. I feel a considerable improvement in my familiarity with Scilab and SIP.

For this activity I give myself a 10/10. I think this entry has achieved a simpler format in comparison to the previous discussions. Thanks to "Mother Earth" (mom) for reading this blog and telling me that she cannot understand anything. It has compelled me to write simpler and I hope to post more entries in the future.