Introduction
In this article we will do a walkthrough of the smaller problems we had to solve in order to get to a good end result. The end result in itself is too specific to be of interest, but the different smaller problems might just be something you need somewhere down the road. I will only post code inline (and not as an article resource) because of this nature and various other reasons (you need a datasource, ...).
Note that this comes from a prototype. Chances are I will change things when I will finally integrate in the real application.
Background
I was put on an existing project were I was asked to integrate surface calculation of sunspots based on scanned drawings. Since I wasn't familiar with the existing project and also because I saw right away this wouldn't be straightforward I chose to write a prototype first.
A sunspot drawing is made by putting a template paper under a telescope and manually drawing the spots (with a pencil).
A drawing looks like this: http://sidc.oma.be/DATA/uset/archdrawings/2014/08/usd201408270840.jpg
The Problems
This project posed several problems.
- How to load the image from disc?
- How to detect the sunpots?
- How to calculate their surface?
- How to get out unwanted features?
Each of these problems had smaller problems to deal with as well and we'll dive deeper into the rabbit hole to discover (and solve) those as well.
The Sub-Problems
How to load the image from disc?
An easy problem to start with.
You start by creating a System.Drawing.Bitmap object.
System.Drawing.Bitmap bmp = new System.Drawing.Bitmap("[path to file]");
And I passed it to this function:
private void SetModifiedImage(System.Drawing.Bitmap bmp, bool newselected){
bmp = CreateNonIndexedImage(bmp);
BitmapImage bmpi = new BitmapImage();
bmpi.BeginInit();
bmpi.StreamSource = new System.IO.MemoryStream((byte[])(new System.Drawing.ImageConverter()).ConvertTo(bmp, typeof(byte[])));
bmpi.EndInit();
image_selected.Source = bmpi;
}
This is nothing exotic, except for the CreateNonIndexedImage function. This was one of the problems I encountered during a filter I had to apply which was in turn a problem to be solved with "how to detect sunspots".
To make sure the problem was solved right away I did so on loading the image, hence the function. I will dive deeper into this later on.
The image_selected object is a System.Windows.Controls.Image object from WPF. Similar to a picturebox in Winforms.
This is how the original looks like:
How to detect the sunspots?
Now for the serious work. To detect the features on the disc I basically thought of 2 options.
- feature detection (edges, circles, ...) and calculate the surface of that.
- counting "dirty" pixels.
I was very doubtful the first method would work very well, since many of the features are NOT sunpots and you would end up in (manually) grouping features since we needed the surface of an entire group, not a single sunspot. With selecting option 2 I immediately went down a road of other problems, but I saw the end of a tunnel, whereas with option 1 I had no clue where I would end up.
Quote:
Advice: Never start anything if you're not sure where it will end up.
Counting dirty pixels. Since this image is a "grayscale" it is kind of hard to determine which pixels are "dirty" (and thus a feature) and which ones aren't, but we can assist the computer by doing something I called "Contrast Fix" in my application. Basically you define a (configurable) treshold and loop through the image assigning a pixel a complete black or complete white color based on whether the initial value was below or above this treshold. There are two options to do this:
- GetPixel / SetPixel
- LockBits/UnLockBits
At first I used option 1 inside a designated area I indicated on the image, but for other reasons I actually had to do the Contrast Fix on the entire image and here option 1 is way too slow.
Here is the code for the Contrast Fix:
(note I used lightgray instead of complete white as color 
)
private void HighLightDirtyPixelsLockBits(){
string fullpath = System.Configuration.ConfigurationManager.AppSettings["workingfolder"] + "\\" + listbox_files.SelectedItem;
System.Drawing.Bitmap bmp = (System.Drawing.Bitmap)System.Drawing.Bitmap.FromFile(fullpath);
System.Drawing.Bitmap tmpbmp = bmp = CreateNonIndexedImage(bmp);
System.Drawing.Rectangle tmprect = new System.Drawing.Rectangle(0, 0, bmp.Width, bmp.Height);
System.Drawing.Imaging.BitmapData bmpdata = tmpbmp.LockBits(tmprect, System.Drawing.Imaging.ImageLockMode.ReadWrite, tmpbmp.PixelFormat);
IntPtr pointer = bmpdata.Scan0;
int nrofbytes = Math.Abs(bmpdata.Stride) * bmpdata.Height; //stride = scan width in bytes (not pixels)
byte [] bytes = new byte[nrofbytes];
System.Runtime.InteropServices.Marshal.Copy(pointer, bytes, 0, nrofbytes);
System.Drawing.Color color_black = System.Drawing.Color.Black;
System.Drawing.Color color_lightgray = System.Drawing.Color.LightGray;
long count = 0;
int r,g,b;
for(int i = 0; i < bytes.Length; i += 4){
b = (int)bytes[i];
g = (int)bytes[i+1];
r = (int)bytes[i+2];
if(b > treshold && g > treshold && r > treshold){
bytes[i] = color_black.B;
bytes[i+1] = color_black.G;
bytes[i+2] = color_black.R;
bytes[i+3] = color_black.A;
} //end if
else{
count++;
bytes[i] = color_lightgray.B;
bytes[i+1] = color_lightgray.G;
bytes[i+2] = color_lightgray.R;
bytes[i+3] = color_lightgray.A;
} //end else
} //end for
System.Runtime.InteropServices.Marshal.Copy(bytes, 0, pointer, nrofbytes);
tmpbmp.UnlockBits(bmpdata);
bmp_contrastfix = (System.Drawing.Bitmap)tmpbmp.Clone();
SetModifiedImage(tmpbmp, false);
}
And here is where we ran into the problem that the CreateNonIndexedImage function, mentioned above ,solves. The problem was that the HighLightDirtyPixelsLockBits crashed on some images. The reason for this crash was the PixelFormat. Some images have a Format24bppRgb format while some have a Format8bppIndexed format. The CreateNonIndexedImage function makes sure that all images have the same PixelFormat in memory. take note that I use i += 4 in the for loop, because of the ARGB (32 bit) nature of the memory bitmap. If you have another format you might need to change that. eg a 24 bit image (RGB) would only need i+= 3.
It is simple enough:
private System.Drawing.Bitmap CreateNonIndexedImage(System.Drawing.Bitmap bmp){
System.Drawing.Bitmap newbmp = new System.Drawing.Bitmap(bmp.Width, bmp.Height, System.Drawing.Imaging.PixelFormat.Format32bppRgb);
System.Drawing.Graphics g = System.Drawing.Graphics.FromImage(newbmp);
g.DrawImage(bmp, 0, 0);
return newbmp;
}
The HighLightDirtyPixelsLockBits function puts the image in a 2 color mode in a fraction of a second, while the GetPixel/SetPixel method would take several seconds to complete. in fact it is so fast that I used a slider to change the treshold if the result wasn't good enough and I change the contrast while changing the slider's value.
The user still has to indicate manually where the sunspots are, currently with drawing a rectangle on the image (mousedown, mousemove, mouseup events on the Image control) and with the GetPixel/SetPixel method (looping through the drawn rectangle's coordinates instead of the entire image). With that we can highlight the "dirty" pixels belonging to the sunspots.
Drawing that rectangle proved to be a pain in the *ss. The reason for this pain is that you need to translate the mouse coordinates to coordinates on the image, else your rectangle does not end up where you expect it. I'm not sure I solved that the best way, but it works.
The MouseDown event handler:
private void image_selected_MouseLeftButtonDown(object sender, System.Windows.Input.MouseButtonEventArgs e) {
if(drawstate == DRAWSTATE.RECTANGLE){
gaugepoint = image_selected.PointToScreen(new System.Windows.Point(0,0));
double x = e.GetPosition((System.Windows.Controls.Image)sender).X + gaugepoint.X + correction.X;
double y = e.GetPosition((System.Windows.Controls.Image)sender).Y + gaugepoint.Y + correction.Y;
start = new System.Drawing.Point((int)System.Math.Floor(x), (int)System.Math.Floor(y));
Rect r = new Rect(start.X + sv.HorizontalOffset, start.Y + sv.VerticalOffset, 0, 0);
System.Windows.Shapes.Rectangle newrect = new System.Windows.Shapes.Rectangle();
newrect.Height = r.Height;
newrect.Width = r.Width;
newrect.Stroke = new SolidColorBrush(System.Windows.Media.Colors.Red);
newrect.MouseLeftButtonUp += newrect_MouseLeftButtonUp;
Canvas.SetLeft(newrect, r.X);
Canvas.SetTop(newrect, r.Y);
c.Children.Add(newrect);
} }
as you can see the mouse coordinates are translated with a correction factor and a Rectangle shape is added to the canvas object that also holds the Image object. The sv variable is the correction for the Scrollbar position.
I will discuss the DRAWSTATE if statement later on in this article.
The MouseMove handler:
private void image_selected_MouseMove(object sender, System.Windows.Input.MouseEventArgs e) {
if(drawstate == DRAWSTATE.RECTANGLE){
if(e.LeftButton == System.Windows.Input.MouseButtonState.Pressed){
double x = e.GetPosition((System.Windows.Controls.Image)sender).X + gaugepoint.X + correction.X;
double y = e.GetPosition((System.Windows.Controls.Image)sender).Y + gaugepoint.Y +correction.Y;
end = new System.Drawing.Point((int)System.Math.Floor(x), (int)System.Math.Floor(y));
Rect r = new Rect(start.X, start.Y, (int)Math.Abs(end.X-start.X), (int)Math.Abs(end.Y-start.Y));
((System.Windows.Shapes.Rectangle)c.Children[c.Children.Count-1]).Width = r.Width;
((System.Windows.Shapes.Rectangle)c.Children[c.Children.Count-1]).Height = r.Height;
} } }
Which only updates the rectangle's size.
and finally the MouseUp event handler:
private void image_selected_MouseLeftButtonUp(object sender, System.Windows.Input.MouseButtonEventArgs e) {
if(drawstate == DRAWSTATE.RECTANGLE){
double x = e.GetPosition((System.Windows.Controls.Image)sender).X + gaugepoint.X + correction.X;
double y = e.GetPosition((System.Windows.Controls.Image)sender).Y + gaugepoint.Y + correction.Y;
end = new System.Drawing.Point((int)System.Math.Floor(x), (int)System.Math.Floor(y));
Rect r = new Rect(start.X + sv.HorizontalOffset, start.Y + sv.VerticalOffset, (int)Math.Abs(end.X-start.X), (int)Math.Abs(end.Y-start.Y));
((System.Windows.Shapes.Rectangle)c.Children[c.Children.Count-1]).Stroke = new SolidColorBrush(System.Windows.Media.Colors.LightGreen);
((System.Windows.Shapes.Rectangle)c.Children[c.Children.Count-1]).Width = r.Width;
((System.Windows.Shapes.Rectangle)c.Children[c.Children.Count-1]).Height = r.Height;
System.Windows.Point one = new System.Windows.Point(start.X - gaugepoint.X, start.Y-gaugepoint.Y-correction.Y);
System.Windows.Point two = new System.Windows.Point(end.X - gaugepoint.X, end.Y-gaugepoint.Y-correction.Y);
bmprect = new Rect(one, two);
BusinessObjects.Shapes.Square square = new BusinessObjects.Shapes.Square();
square.Id = c.Children.Count-1;
square.Start = new System.Windows.Point(bmprect.X, bmprect.Y);
square.Size = new System.Windows.Point(bmprect.Width, bmprect.Height);
squares.Add(square);
HighlightDirtyPixels(bmprect, square);
} else{
if(drawstate == DRAWSTATE.POLYLINE){
double x = e.GetPosition((System.Windows.Controls.Image)sender).X + gaugepoint.X;
double y = e.GetPosition((System.Windows.Controls.Image)sender).Y + gaugepoint.Y; if(tmppolygon == null){
tmppolygon = new System.Windows.Shapes.Polygon();
c.Children.Add(tmppolygon);
} tmppolygon.Points.Add(new System.Windows.Point(x, y));
tmppolygon.Stroke = new SolidColorBrush(System.Windows.Media.Colors.Red);
} } }
Note that I put the rectangle in green here (it started as red), that's just a detail for the user, so he knows his rectangle is "ready". (and also for debug purposes, see next paragraph)
One problem here is that the MouseUp event did not fire most of the time and I ended up with multiple (red) rectangles which didn't do things as expected. The reason for this is because we update the drawn rectangle during mousemove, when you release the mousebutton you're on top of the rectangle, NOT on top of the image_selected control. Hence I added:
private void newrect_MouseLeftButtonUp(object sender, System.Windows.Input.MouseButtonEventArgs e) {
image_selected_MouseLeftButtonUp(image_selected, e);
}
Also you'll notice the there is a drawstate polyline, which is to solve a problem I will mention later on.
And here's the function that uses the GetPixel/SetPixel method:
private void HighlightDirtyPixels(Rect r, BusinessObjects.Shapes.Square square){
long count = 0;
long pixelindex = 0;
string msg = "";
System.Drawing.Bitmap bmp = new System.Drawing.Bitmap(bmp_contrastfix);
try{
System.Drawing.Color color_black = System.Drawing.Color.Black;
square.Pixels = new List<BusinessObjects.Shapes.Pixel>();
r.X -= correction.X;
r.Y += 2;
BusinessObjects.Objects.SunspotArea ssa = new BusinessObjects.Objects.SunspotArea();
ssa.Id = square.Id;
for (int x = (int)r.X; x < (int)r.X + r.Width; x++){
for (int y = (int)r.Y; y < (int)r.Y + r.Height; y++){
BusinessObjects.Shapes.Pixel p = new BusinessObjects.Shapes.Pixel();
p.Id = pixelindex;
p.Index = new System.Windows.Point(x, y);
double x_from_center = x - (int)sun.CenterPixel.X;
double y_from_center = y - (int)sun.CenterPixel.Y;
p.DistanceToCenter = Math.Sqrt(Math.Pow(x_from_center, 2) + Math.Pow(y_from_center, 2));
bool insidecircle = p.DistanceToCenter <= sun.RadiusInPixels;
System.Drawing.Color c = bmp.GetPixel(x, y);
bool insidepolygon = PointInPolygon(tmppolygon, new System.Windows.Point(x, y));
bool draworange = insidepolygon || drawstate == DRAWSTATE.RECTANGLE;
if( (c.R > color_black.R && c.G > color_black.G && c.B > color_black.B) && insidecircle && draworange){
bmp.SetPixel(x, y, System.Drawing.Color.Orange);
p.PixelColor = System.Drawing.Color.Orange;
p.DeprojectedSurfaceInMeters = sun.SquareMeterPerPixel * p.DeprojectedSurfaceInPixels;
ssa.DeprojectedSurfaceInMeters += p.DeprojectedSurfaceInMeters;
count++;
msg = "orange\t" + "\tx:" + x + "\ty:" + y;
} else{
if(drawstate == DRAWSTATE.POLYLINE && !insidepolygon){
bmp.SetPixel(x, y, System.Drawing.Color.Black);
p.PixelColor = System.Drawing.Color.Black;
msg = "black\t" + "\tx:" + x + "\ty:" + y;
} if(drawstate == DRAWSTATE.POLYLINE && insidepolygon){
bmp.SetPixel(x, y, System.Drawing.Color.White);
p.PixelColor = System.Drawing.Color.White;
msg = "white\t" + "\tx:" + x + "\ty:" + y;
} } pixelindex++;
square.Pixels.Add(p);
} } ssa.SurfaceInPixels = count;
ssa.SurfaceInMeters = ssa.SurfaceInPixels * sun.SquareMeterPerPixel;
ssa.DeprojectedSurfaceInPixels = ssa.CalculatePixelSurface(square, sun, BusinessObjects.Objects.SunspotArea.UNIT.Pixel);
ssa.DeprojectedSurfaceInMeters = ssa.CalculatePixelSurface(square, sun, BusinessObjects.Objects.SunspotArea.UNIT.Meter);
ssa.DeprojectedSurfaceMsh = ssa.CalculatePixelSurface(square, sun, BusinessObjects.Objects.SunspotArea.UNIT.Msh);
sunspotareas.Add(ssa);
} catch(Exception ex){
WriteInfo(msg + "\tError: " + ex.Message);
} bmp_contrastfix = new System.Drawing.Bitmap(bmp);
SetModifiedImage(bmp, false);
}
This asks for some explanation:
The bmp_contrastfix variable is the object that holds the image after the contrast fix. There are also different objects used here:
- sun (BusinessObjects.Objects.Sun)
- ssa (BusinessObjects.Objects.SunspotArea)
- p (BusinessObjects.Shapes.Pixel)
They are used for calculating different properties along the way, the class definitions are in the References section below.
So let's take a moment to reflect. We loaded the image in memory, we fixed it's pixelformat in order to use the lockbits/unlockbits method for a Contrast Fix and then we can draw a rectangle on the image to use the GetPixel/SetPixel method inside that smaller area to highlight the sunspot (group). We fixed several problems along the way:
- determined what defines a "dirty" pixel
- Cleaned up the issue with the PixelFormat.
- Fixed a very annoying and subtle bug where the mouse was over the rectangle during the mouseup event instead of over the image.
How to calculate their surface?
All the previous work was done in order to do exactly this. Calculating the area of a sunspot(group). So let's start counting the "dirty" pixels inside a drawn rectangle and see how many pixels we have inside the circle (= the solar surface) and we have it! Yeah!
Wait, ... circle? euhm, where is that circle exactly ? Yes, you can see it clearly on the image, but your program doesn't know if the circle's edge is part of a sunspot or something else...
Welcome, to
Problem 1 - Detect the circle that defines the Sun's limb
This is surprisingly easy, . . . if you know how. I tried 2 libraries AForge.Net and Emgu (a wrapper around OpenCV). the Emgu method was clearly better, though I had to provide some manual methods to move and resize the detected circle a bit to fit the drawing "exactly".
This is the detection method:
private void btn_shapedetectionemgu_Click(object sender, RoutedEventArgs e) {
Emgu.CV.Image<Emgu.CV.Structure.Bgr, byte> img;
if(bmp_contrastfix != null){
img = new Image<Bgr,byte>(bmp_contrastfix);
} else{
img = new Image<Bgr,byte>(bmp_ori);
}
Emgu.CV.Image<Emgu.CV.Structure.Gray, byte> img_gray = img.Convert<Emgu.CV.Structure.Gray, byte>().PyrDown().PyrUp();
Gray cannytreshold = new Gray(int.Parse(txtbox_cannytresholdintensity.Text));
Gray circleaccumulatortreshold = new Gray(int.Parse(txtbox_accumulatortresholdintensity.Text));
CircleF [] circles = img_gray.HoughCircles(cannytreshold, circleaccumulatortreshold, double.Parse(txtbox_resolution.Text),double.Parse(txtbox_minimumdistance.Text),int.Parse(txtbox_minimumradius.Text),int.Parse(txtbox_maximumradius.Text))[0];
for(int i = 0; i < circles.Length; i++){
System.Windows.Point center = new System.Windows.Point(circles[i].Center.X, circles[i].Center.Y);
Cross2DF cross = new Cross2DF(new System.Drawing.PointF((float)center.X, (float)center.Y), 10, 10);
if(IsCircleInDrawing(center, circles[i].Radius, img.Width, img.Height)){
img.Draw(circles[i], new Bgr(System.Drawing.Color.Magenta), 2);
img.Draw(cross, new Bgr(System.Drawing.Color.Magenta), 2);
BusinessObjects.Shapes.Circle c = new BusinessObjects.Shapes.Circle();
c.Center = new System.Windows.Point(circles[i].Center.X, circles[i].Center.Y);
c.Radius = (int)circles[i].Radius;
c.Id = i;
detectedcircles.Add(c);
sun = new BusinessObjects.Objects.Sun(radiusinmeter, c.Radius);
sun.CenterPixel = new System.Windows.Point(c.Center.X, c.Center.Y);
} } SetModifiedImage(img.ToBitmap(), false);
}
Most of it are Emgu specific objects and I refer to their page for more info: http://www.emgu.com/wiki/index.php/Main_Page and many samples can be found on the internet.
With the HoughCircles method you can find circles in an image. ALL circles, but you can play with the properties a bit to clear most if it out. (eg the circles minimum size, ...) It also finds circles that are partly out of the image so I use a IsCircleInDrawing method which looks like this:
private bool IsCircleInDrawing(System.Windows.Point center, float radius, int bmp_width, int bmp_height){
bool indrawing = false;
if( center.X-radius > 0 && center.X+radius < bmp_width &&
center.Y-radius > 0 && center.Y+radius < bmp_height){
indrawing = true;
} return indrawing;
}
Which only takes circles that are completely inside the image boundaries.
And effectively, it almost always finds the correct circle. In the rare case it doesn't you can change the contrast fix and try to detect again.
This is how the contrast fix and shape detection looks like:
With that you have the center point of the Sun as well as it's radius. You need that in order calculate the surface. Here's where the fun part starts. The drawing is a circle and thus the logical thought is to count how much percent of the sunspot area is covering the disc. Since you know the diameter of the sun (http://en.wikipedia.org/wiki/Sun) and the radius in pixels you can calculate the square meters by adding up all the dirty pixels, right?
WRONG!
Problem 2 - The Sun is a sphere
The surface you calculate there is a projected surface. The Sun is a sphere not a disc and so the surface "increases" with increasing distance from the center of the disc. This part is not easily found on the internet.
Here is the code part:
public double CalculatePixelSurface(BusinessObjects.Shapes.Square square, BusinessObjects.Objects.Sun sun, UNIT unit){
double surface = 0;
double X = 0;
double Y = 0;
for(int i = 0; i < square.Pixels.Count; i++){
if(square.Pixels[i].PixelColor == System.Drawing.Color.Orange){
X = (square.Pixels[i].Index.X-sun.CenterPixel.X) / sun.RadiusInPixels;
Y = (square.Pixels[i].Index.Y-sun.CenterPixel.Y) / sun.RadiusInPixels;
double cosalpha = Math.Sqrt(1-Math.Pow(X, 2)-Math.Pow(Y, 2)); cosalpha = (cosalpha > 16 ? 16 : cosalpha); surface += ((1/Math.Pow(sun.RadiusInPixels, 2)) * (1/(2*Math.PI*cosalpha)));
} } X = ((square.Start.X+(square.Size.X)/2)-sun.CenterPixel.X) / sun.RadiusInPixels;
Y = ((square.Start.Y+(square.Size.Y)/2)-sun.CenterPixel.Y) / sun.RadiusInPixels;
double avgcos = Math.Sqrt(1-Math.Pow(X, 2)-Math.Pow(Y, 2));
Angle = (180/Math.PI)*Math.Acos(avgcos);
switch(unit){
case UNIT.Fraction: break;
case UNIT.Meter: surface *= sun.SquareMeterPerPixel; break;
case UNIT.Msd: break;
case UNIT.Msh: surface *= 1000000; break;
case UNIT.Pixel: break;
}
return surface;
}
Here we loop through a Square object that holds the Sunspot and calculates, for each pixel that is "dirty" it's de-projected surface. Once you have the formula, the method is not too hard. Getting to that formula however is not that easy. Furthermore I take the center of the square and calculate it's angle against the center to check if it is more or less correct with what I selected. eg. If I take a square near the center and I end up with 78 degrees for example, that is wrong. Mostly the area is expressed as a fraction of the disc in millionth of a solar hemisphere (Msh).
The mathematics.
This is how the raw formula looks like:
For a non mathematician (like me) this is useless, so I talked to a colleague who explained things to me.
In a more coding manner this looks like:
(1/N^2) * (1/(2*PI*cos(alpha)) = (1/N^2) * (1/(2*PI*SQRT(1-x^2-y^2))
How to get there?
First thing you need to do is scale the sphere to 1,1. (That means that a radius of 450 pixels eg, = 1).
N = # pixels of the radius.
We need this scaling because it makes things easier in finding the formula. eg. if the surface of a circle = Pi.r² you can say that for a radius 1 the surface is just Pi (r² = 1² = 1).
The surface (on the circle) of a sunspot area is the sum of all pixels. To get this same area on the scaled circle you get:
sum of pixels / N².
and as a fraction of the disc it becomes:
sum of pixels / (Pi*N²)
The previous formulas are simple, but still in 2D. We are not working in 2D, but 3D and the formula for the surface of a sphere is:
4*Pi*r²
Great. We have r = 1 so that's reduced to 4*Pi, but also note we're not working on a full sphere, just half of it. So the circle we see in the drawing has a (scaled) surface of:
2*Pi
How to get the surface of just a part of our sphere?
For this you need to know how far the sunspot is located relative to the center. In human terms, the pixel in the center of the disc (the circle) has a surface of 1 pixel (That is, the pixel has an equal "surface value" in 2D and 3D) The more your pixel is moved away from the center, the more it is "projected" and thus, the more you need to correct. In order to get this "de-projected" value you'll need to know the angle (alpha) of the line pixel/solarcenter with the Z-axis. That is correct, the Z-axis comes out of (of goes into) the page. We have no information yet about anything concerning this third axis (and the 3D properties of the drawing).
It would look something like this:
(Note that the axes have turned compared to the previous (2D) drawing!)
The red cros is the point we want to calculate. The line C (full line) is from the center to the pixel located on the sphere. A is the line projecting the pixel in the XY plane. B is the distance from that projected pixel towards the center. alpha is the angle between the Z-axis and the line C.
And if you remember from school:
cos(alpha) = B / C
Where B equals the adjacent side of the right-angled triangle and C the hyptenuse. (http://en.wikipedia.org/wiki/Trigonometric_functions)
And C we know, this is the radius of the sphere. Scaled this gives 1. Thus we only need to find the value for B.
In our case the B in the formula is the line A. Which is just the value on the Z-axis. Since we only have x,y information we need to find the z-value.
In a sphere counts:
X² + Y² + Z² = r²
Since the radius is 1 we have:
X² + Y² + Z² = 1
=>
Z² = 1 - X² - Y²
=>
Z = SQRT(1 - X² - Y²)
So cosine of alpha becomes:
cos(alpha) = SQRT(1 - X² - Y²) / 1
So here the fraction of the surface of a sunspot(group) for sphere becomes:
sumofpixels/N² * 1/2*Pi*cos(alpha)
=>
sumofpixels/N² * 1/2*Pi*SQRT(1-X²-Y²)
so what we do in code is for each pixel:
1/N² * 1/2*Pi*SQRT(1-X²-Y²)
and add those (the surface value of each pixel) up for a total.
In code this becomes (CalculatePixelSurface function):
for(int i = 0; i < square.Pixels.Count; i++){
if(square.Pixels[i].PixelColor == System.Drawing.Color.Orange){
X = (square.Pixels[i].Index.X-sun.CenterPixel.X) / sun.RadiusInPixels;
Y = (square.Pixels[i].Index.Y-sun.CenterPixel.Y) / sun.RadiusInPixels;
double cosalpha = Math.Sqrt(1-Math.Pow(X, 2)-Math.Pow(Y, 2)); cosalpha = (cosalpha > 16 ? 16 : cosalpha); surface += ((1/Math.Pow(sun.RadiusInPixels, 2)) * (1/(2*Math.PI*cosalpha)));
} }
Note the check against the cosine value (cosalpha = (cosalpha > 16 ? 16 : cosalpha) that's because the Area AtDiskCenter will go to infinity as the pixel goes toward the eastern or western limb. In order to avoid this, we use the Paul Higgins's approach to limit the correction factor (i.e., the terms of the sum) to 16, which corresponds to a LOS angle of ' 86.
Update 14/03/2016: After some comparisons we noticed an error in the results, leading back to the line:
cosalpha = (cosalpha > 16 ? 16 : cosalpha);
We changed it to
double tresholdradians = (2*Math.PI/360.0)*86.0;
cosalpha = (cosalpha < Math.Cos(tresholdradians) ? Math.Cos(tresholdradians) : cosalpha);
which solved the problem.
How to get out unwanted features?
1. Point In Polygon Algorithm
By using a rectangle you're somewhat limited in your selection method. What if the sunspot group was on the edge of the circle (you select the edge with it and those are "dirty" pixels as well). You also have unwanted features like lines, filaments, a cross, ...) that sometimes creep in the selection.
We can get out anything outside the circle (part of HightLightDirtyPixels method):
p.Index = new System.Windows.Point(x, y);
double x_from_center = x - (int)sun.CenterPixel.X;
double y_from_center = y - (int)sun.CenterPixel.Y;
p.DistanceToCenter = Math.Sqrt(Math.Pow(x_from_center, 2) + Math.Pow(y_from_center, 2));
bool insidecircle = p.DistanceToCenter <= sun.RadiusInPixels;
basically we calculate for each pixel in the selected rectangle if that pixel is within the radius of the sun.
That still leaves us with the other unwanted features and therefore we did research on using a polygon instead of the rectangle. The user can manually define a polygon and calculate the dirty pixels inside that polygon. It is easy enough detecting a boundary rectangle to keep the loop short enough. In the prototype I developed I used a button to switch between polygon and rectangle. here's that code:
private void btn_usepolyline_Click(object sender, RoutedEventArgs e) {
Button tmpbtn = (Button)sender;
if( tmpbtn.Content.ToString() == "Polyline Off"){
tmpbtn.Content = "Polyline On";
drawstate = DRAWSTATE.POLYLINE;
} else{
if( tmpbtn.Content.ToString() == "Polyline On"){
tmpbtn.Content = "Polyline Off";
var x = (from p in tmppolygon.Points select p.X).Min();
var y = (from p in tmppolygon.Points select p.Y).Min();
var width = (from p in tmppolygon.Points select p.X).Max() - x;
var height = (from p in tmppolygon.Points select p.Y).Max() - y;
Rect r = new Rect(x + sv.HorizontalOffset, y + sv.VerticalOffset, width, height);
r.X -= 20;
r.Width = += 20; System.Windows.Shapes.Rectangle newrect = new System.Windows.Shapes.Rectangle();
newrect.Height = r.Height;
newrect.Width = r.Width;
newrect.Stroke = new SolidColorBrush(System.Windows.Media.Colors.Blue);
newrect.MouseLeftButtonUp += newrect_MouseLeftButtonUp;
Canvas.SetLeft(newrect, r.X);
Canvas.SetTop(newrect, r.Y);
c.Children.Add(newrect);
BusinessObjects.Shapes.Square square = new BusinessObjects.Shapes.Square();
square.Id = c.Children.Count-1;
square.Start = new System.Windows.Point(r.X, r.Y);
square.Size = new System.Windows.Point(r.X+r.Width, r.Y+r.Height);
squares.Add(square);
HighlightDirtyPixels(r, square);
drawstate = DRAWSTATE.RECTANGLE;
tmppolygon = null;
} } }
Note that we make the bounding rectangle slightly larger on the left. This is due to the fact that the algorithm always would skip a few pixels on the left side of the polygon (see image below) . Making the rectangle a little larger fixes this problem.
So the next difficulty is to determine if a pixel is inside a polygon or not. here's an interesting article about it: http://alienryderflex.com/polygon/
The code, once you find it, is easy enough:
private bool PointInPolygon(System.Windows.Shapes.Polygon polygon, System.Windows.Point point){
int sides = polygon.Points.Count;
bool pinp = false;
int j = sides-1;
double [] point_x = (from p in polygon.Points select p.X).ToArray();
double [] point_y = (from p in polygon.Points select p.Y).ToArray();
for(int i = 0; i < sides; i++){
if( (point_y[i] < point.Y && point_y[j] >= point.Y || point_y[j] < point.Y && point_y[i] >= point.Y) &&
(point_x[i] <= point.X || point_x[j] <= point.X)){
pinp ^= (point_x[i]+(point.Y-point_y[i])/(point_y[j]-point_y[i])*(point_x[j]-point_x[i]) < point.X);
} j = i;
} return pinp;
}
int sides =polygon.Points.Count is not completely correct if your edges cross one another, but in our case the user should define a polygon AROUND something, so lines should not cross.
And here a final image where the sunspot was selected with a polygon:
(Note: the small black tip on the right of the polygon issue is solved by making the bounding rectangle slightly larger)
2. Image Mask Method
An alternative way of finding the pixels inside a polygon is by making a copy of the image and draw a filled polygon on it. (Graphics.FillPolygon method). By filling the polygon with a weird color value (let's say Aqua) you can loop through the same bounding rectangle and check if the color has an Aqua color (inside polygon) or another color (outside polygon). This is how it looks like...
private void HighlightDirtyPixelsPolygon(Rect r, BusinessObjects.Shapes.Square square){
long count = 0;
long pixelindex = 0;
string msg = "";
System.Drawing.Bitmap bmp = new System.Drawing.Bitmap(bmp_contrastfix);
System.Drawing.Bitmap bmp_mask = new System.Drawing.Bitmap(bmp);
if(tmppolygon != null){
System.Drawing.Graphics g = System.Drawing.Graphics.FromImage(bmp_mask);
System.Drawing.Point [] points = new System.Drawing.Point[tmppolygon.Points.Count];
for(int i = 0; i < tmppolygon.Points.Count; i++){
points[i].X = (int) tmppolygon.Points[i].X;
points[i].Y = (int) tmppolygon.Points[i].Y;
}
System.Drawing.Color c = System.Drawing.Color.Aqua;
g.FillPolygon(new System.Drawing.SolidBrush(c), points);
}
try{
System.Drawing.Color color_black = System.Drawing.Color.Black;
square.Pixels = new List<BusinessObjects.Shapes.Pixel>();
r.X -= correction.X;
r.Y += 2;
BusinessObjects.Objects.SunspotArea ssa = new BusinessObjects.Objects.SunspotArea();
ssa.Id = square.Id;
for (int x = (int)r.X; x < (int)r.X + r.Width; x++){
for (int y = (int)r.Y; y < (int)r.Y + r.Height; y++){
BusinessObjects.Shapes.Pixel p = new BusinessObjects.Shapes.Pixel();
p.Id = pixelindex;
p.Index = new System.Windows.Point(x, y);
double x_from_center = x - (int)sun.CenterPixel.X;
double y_from_center = y - (int)sun.CenterPixel.Y;
p.DistanceToCenter = Math.Sqrt(Math.Pow(x_from_center, 2) + Math.Pow(y_from_center, 2));
bool insidecircle = p.DistanceToCenter <= sun.RadiusInPixels;
System.Drawing.Color c = bmp.GetPixel(x, y);
bool insidepolygon = PointInPolygon2(bmp_mask, new System.Drawing.Point(x, y));
if( (c.R > color_black.R && c.G > color_black.G && c.B > color_black.B) && insidecircle && insidepolygon){
bmp.SetPixel(x, y, System.Drawing.Color.Orange);
p.PixelColor = System.Drawing.Color.Orange;
p.DeprojectedSurfaceInMeters = sun.SquareMeterPerPixel * p.DeprojectedSurfaceInPixels;
ssa.DeprojectedSurfaceInMeters += p.DeprojectedSurfaceInMeters;
count++;
msg = "orange\t" + "\tx:" + x + "\ty:" + y;
} else{
if(insidepolygon){
bmp.SetPixel(x, y, System.Drawing.Color.White);
p.PixelColor = System.Drawing.Color.White;
msg = "white\t" + "\tx:" + x + "\ty:" + y;
} } pixelindex++;
square.Pixels.Add(p);
} } ssa.SurfaceInPixels = count;
ssa.SurfaceInMeters = ssa.SurfaceInPixels * sun.SquareMeterPerPixel;
ssa.DeprojectedSurfaceInPixels = ssa.CalculatePixelSurface(square, sun, BusinessObjects.Objects.SunspotArea.UNIT.Pixel);
ssa.DeprojectedSurfaceInMeters = ssa.CalculatePixelSurface(square, sun, BusinessObjects.Objects.SunspotArea.UNIT.Meter);
ssa.DeprojectedSurfaceMsh = ssa.CalculatePixelSurface(square, sun, BusinessObjects.Objects.SunspotArea.UNIT.Msh);
sunspotareas.Add(ssa);
} catch(Exception ex){
WriteInfo(msg + "\tError: " + ex.Message);
} bmp_contrastfix = new System.Drawing.Bitmap(bmp);
SetModifiedImage(bmp, false);
}
Note here we take a copy of the bmp object into bmp_mask and use the Grapics object to draw the filled polygon into the image (compared to on top of the image). We can loop through the bounding rectangles coordinates and check the color of the mask to know if the pixel on the original is inside the polygon (drawn on the mask).
The PointInPolygon2 method is simplified in comparison to the PointInPolygon method:
private bool PointInPolygon2(System.Drawing.Bitmap bmp, System.Drawing.Point point){
bool pinp = false;
System.Drawing.Color c = bmp.GetPixel(point.X, point.Y);
if(c.A == 255 && c.R == 0 && c.G == 255 && c.B == 255){ pinp = true;
}
}
Both the "Point In Polygon Algorithm" as the "Image Mask Method" seem to work equally well, though the second method takes more memory.
Points of Interest
- When I started this project I had no experience in image manipulation, I learned this stuff on the go.
- When doing complex stuff like this always start out with a prototype until you're satisfied. Then, and only then, put it in the real application.
- You need to understand that for software, when we did the contrast fix, all "white parts" are considered "dirty". That includes anything outside the disc, the disc edge itself, all the sunspots and all the unwanted features.
- This still is a prototype. I needed a lot of help from various people (including here on codeproject) to get it working. However, bugs can still be in the code somewhere. (there is still one in the polyline method of selecting the sunspot groups fro example)
- One thing at a time. As you saw, some changes had it's repercussions on code written earlier. eg. The CreateNonIndexedImage function came only after I encountered errors during the "contrast fix" and was used on the image load (practically the first thing I wrote) For that reason make sure to tackle things one at a time. Even then, things can become pretty complex very quickly. Trying to do too many things at once is a road to disaster.
- What did I (we) learn during this process:
- The PixelFormat of an image can mess up the lockbits/unlockbits method. even if it doesn't mess it up (raises and exception) you still need to know this information for the for-loop.
- Mouse events work as expected, but not necessarely as you expected. (the example here was the mouseup event triggering on the square object instead of the image object.
- Rectangle "selection" is the easiest way, but not very accurate.
- Feature detection with Emgu. We managed to get out unwanted (detected) features, by changing the properties of the search and limit the results by checking circles completely inside the drawing.
- A lot of mathematics since we know the sun is a sphere (and not a circle). I very much enjoyed the elegance in the reasoning towards the final formula.
Interesting facts
- In our dataset there are about 23000 drawings.
- The oldest drawing dates from March 5th 1940.
- Drawings are still made today (not only from our institute).
History
- September 2014: First version.
References
Objects used
Here are the objects that I used throughout the code. I just added them for reference.
Class Sun
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace BusinessObjects.Objects {
public class Sun {
public enum UNIT{
Meter,
Pixels
};
private int radiusinmeter;
private int radiusinpixels;
private System.Windows.Point centerpixel;
public Sun(int radius_in_meter, int radius_in_pixels){
radiusinmeter = radius_in_meter;
radiusinpixels = radius_in_pixels;
}
public System.Windows.Point CenterPixel{
set{
centerpixel = value;
}
get{
return centerpixel;
}
}
public int RadiusInMeter{
set{
radiusinmeter = value;
} get{
return radiusinmeter;
} }
public int RadiusInPixels{
set{
radiusinpixels = value;
} get{
return radiusinpixels;
} }
public double MeterPerPixel{
get{
return (double)radiusinmeter / (double)radiusinpixels;
} }
public double SquareMeterPerPixel{
get{
double surfacep = TotalSolarSurface(UNIT.Pixels);
double surfacem = TotalSolarSurface(UNIT.Meter);
return surfacem / surfacep;
} }
public double TotalSolarSurface(UNIT unit){
double result = -1;
switch(unit){
case UNIT.Meter: result = 4 * Math.PI * Math.Pow(radiusinmeter, 2); break;
case UNIT.Pixels: result = 4 * Math.PI * Math.Pow(radiusinpixels, 2); break;
} return result;
}
public double TotalSolarVolume(UNIT unit){
double result = -1;
switch(unit){
case UNIT.Meter: result = (4/3) * Math.PI * Math.Pow(radiusinmeter, 3); break;
case UNIT.Pixels: result = (4/3) * Math.PI * Math.Pow(radiusinpixels, 3); break;
} return result;
}
public override string ToString() {
StringBuilder builder = new StringBuilder("");
builder.Append("Meter / Pixel: ").Append(this.MeterPerPixel).Append(Environment.NewLine);
builder.Append("Square Meter / Pixel: ").Append(this.SquareMeterPerPixel).Append(Environment.NewLine);
builder.Append("Center Pixel: ").Append(this.CenterPixel.ToString()).Append(Environment.NewLine);
builder.Append(Environment.NewLine);
builder.Append("Radius (meters): ").Append(this.RadiusInMeter).Append(Environment.NewLine);
builder.Append("Surface (meters): ").Append(this.TotalSolarSurface(BusinessObjects.Objects.Sun.UNIT.Meter)).Append(Environment.NewLine);
builder.Append("Volume (meters): ").Append(this.TotalSolarVolume(BusinessObjects.Objects.Sun.UNIT.Meter)).Append(Environment.NewLine);
builder.Append(Environment.NewLine);
builder.Append("Radius (pixels): ").Append(this.RadiusInPixels).Append(Environment.NewLine);
builder.Append("Surface (pixels): ").Append(this.TotalSolarSurface(BusinessObjects.Objects.Sun.UNIT.Pixels)).Append(Environment.NewLine);
builder.Append("Volume (pixels): ").Append(this.TotalSolarVolume(BusinessObjects.Objects.Sun.UNIT.Pixels)).Append(Environment.NewLine);
return builder.ToString();
}
}
}
Class SunspotArea
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace BusinessObjects.Objects {
public class SunspotArea {
public enum UNIT{
Fraction, Meter,
Msd, Msh, Pixel
};
public int Id{ set; get; }
public double SurfaceInPixels{ set; get; }
public double SurfaceInMeters{ set; get; }
public double DeprojectedSurfaceInPixels { set; get; }
public double DeprojectedSurfaceInMeters { set; get; }
public double DeprojectedSurfaceMsh { set; get; }
public double Angle { set; get; }
public override string ToString() {
StringBuilder builder = new StringBuilder("");
builder.Append("Id: ").Append(this.Id).Append(Environment.NewLine);
builder.Append("Surface (Pixels): ").Append(this.SurfaceInPixels).Append(Environment.NewLine);
builder.Append("Surface (Meters): ").Append(this.SurfaceInMeters).Append(Environment.NewLine);
builder.Append("de-Projected Surface (Pixels - fraction): ").Append(this.DeprojectedSurfaceInPixels).Append(Environment.NewLine);
builder.Append("de-Projected Surface in Meters: ").Append(this.DeprojectedSurfaceInMeters).Append(Environment.NewLine);
builder.Append("de-Projected Surface in Msh: ").Append(this.DeprojectedSurfaceMsh).Append(Environment.NewLine);
builder.Append("Angle (degrees): ").Append(this.Angle).Append(Environment.NewLine);
return builder.ToString();
}
public double CalculatePixelSurface(BusinessObjects.Shapes.Square square, BusinessObjects.Objects.Sun sun, UNIT unit){
double surface = 0;
double X = 0;
double Y = 0;
for(int i = 0; i < square.Pixels.Count; i++){
if(square.Pixels[i].PixelColor == System.Drawing.Color.Orange){
X = (square.Pixels[i].Index.X-sun.CenterPixel.X) / sun.RadiusInPixels;
Y = (square.Pixels[i].Index.Y-sun.CenterPixel.Y) / sun.RadiusInPixels;
double cosalpha = Math.Sqrt(1-Math.Pow(X, 2)-Math.Pow(Y, 2)); double tresholdradians = (2*Math.PI/360.0)*86.0;
cosalpha = (cosalpha < Math.Cos(tresholdradians) ? Math.Cos(tresholdradians) : cosalpha);
surface += ((1/Math.Pow(sun.RadiusInPixels, 2)) * (1/(2*Math.PI*cosalpha)));
} } X = ((square.Start.X+(square.Size.X)/2)-sun.CenterPixel.X) / sun.RadiusInPixels;
Y = ((square.Start.Y+(square.Size.Y)/2)-sun.CenterPixel.Y) / sun.RadiusInPixels;
double avgcos = Math.Sqrt(1-Math.Pow(X, 2)-Math.Pow(Y, 2));
Angle = (180/Math.PI)*Math.Acos(avgcos);
switch(unit){
case UNIT.Fraction: break;
case UNIT.Meter: surface *= sun.SquareMeterPerPixel; break;
case UNIT.Msd: break;
case UNIT.Msh: surface *= 1000000; break;
case UNIT.Pixel: break;
}
return surface;
}
}
}
Class Circle
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace BusinessObjects.Shapes {
public class Circle {
public int Id { set; get; }
public System.Windows.Point Center { set; get; }
public int Radius { set; get; }
}
}
Class Pixel
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace BusinessObjects.Shapes {
public class Pixel {
public long Id { set; get; }
public System.Windows.Point Index { set; get; }
public System.Drawing.Color PixelColor { set; get; }
public double DistanceToCenter { set; get; }
public double DeprojectedSurfaceInMeters { set; get; }
public double DeprojectedSurfaceInPixels { set; get; }
}
}
Class Square
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace BusinessObjects.Shapes {
public class Square {
public int Id { set; get; }
public System.Windows.Point Start { set; get; }
public System.Windows.Point Size { set; get; }
public List<Pixel> Pixels { set; get; }
}
}
Conclusion
This was the first time I had to create such a complex prototype, and it is thanks to experience (and some help) I pulled it of. The best decision was to not start coding in the existing application. Might be obvious for more experienced programmers, but for beginners this is sound advice.
