检测形状数学 OpenCV 图像处理中级开发 Windows C#

太阳黑子表面计算（解决复杂问题的案例研究）

5.00/5 (11投票s)

2014年9月1日

CPOL

17分钟阅读

29345

通过将复杂问题分解为较小问题来解决的案例研究。

引言

在本文中，我们将逐步讲解为了获得一个好的最终结果而需要解决的各个小问题。最终结果本身过于具体，可能不那么有趣，但其中不同的小问题可能正是您将来可能需要的东西。由于这个性质和各种其他原因（您需要数据源，...），我将只在文中内嵌代码（而不是作为文章资源）。

请注意，这是来自一个原型。当我最终将其集成到实际应用程序中时，我可能会更改某些内容。

背景

我被分配到一个现有项目，被要求根据扫描图集成太阳黑子的表面积计算。由于我并不熟悉现有项目，而且我立刻发现这不会是一件简单的事情，所以我选择先编写一个原型。

太阳黑子的图纸是通过将一张模板纸放在望远镜下方，然后手动（用铅笔）绘制黑子来制作的。

图纸看起来是这样的： http://sidc.oma.be/DATA/uset/archdrawings/2014/08/usd201408270840.jpg

问题

这个项目带来了一些问题。

如何从磁盘加载图像？
如何检测太阳黑子？
如何计算它们的表面积？
如何去除不需要的特征？

这些问题中的每一个都有更小的子问题需要处理，我们将深入挖掘，发现（并解决）那些子问题。

子问题

如何从磁盘加载图像？

这是一个容易开始的问题。

您首先创建一个System.Drawing.Bitmap对象。

System.Drawing.Bitmap bmp = new System.Drawing.Bitmap("[path to file]");

然后我将其传递给此函数

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;
}

这没什么特别的，除了CreateNonIndexedImage函数。这是我在应用一个过滤器时遇到的问题之一，而这个过滤器又是需要用“如何检测太阳黑子”来解决的问题。

为了确保问题立即得到解决，我在加载图像时就解决了这个问题，因此有了这个函数。我稍后会更详细地介绍这一点。

image_selected 对象是WPF的System.Windows.Controls.Image对象。类似于Winforms中的picturebox。

原始图像看起来是这样的

Original

如何检测太阳黑子？

现在是真正的工作。为了检测圆盘上的特征，我基本上考虑了两种选择。

特征检测（边缘、圆形等）并计算其表面积。
计算“脏”像素的数量。

我非常怀疑第一种方法是否能很好地工作，因为许多特征并非太阳黑子，而您将不得不（手动）对特征进行分组，因为我们需要的是整个组的表面积，而不是单个太阳黑子的表面积。选择第二种方法立即引出了其他问题，但我看到了隧道的尽头，而对于第一种方法，我则完全不知道会走向何方。

引用

建议：在不确定结果的情况下，切勿开始任何事情。

计算脏像素。由于此图像是“灰度”的，因此很难确定哪些像素是“脏”的（从而成为特征），哪些不是，但我们可以通过执行我应用程序中称为“对比度修复”的操作来协助计算机。基本上，您定义一个（可配置的）阈值，然后遍历图像，根据初始值高于或低于此阈值，为像素赋予完全黑色或完全白色的颜色。有两种选择可以做到这一点：

GetPixel / SetPixel
LockBits/UnLockBits

起初，我在图像上的一个指定区域使用了第一种方法，但由于其他原因，我实际上不得不对整个图像进行对比度修复，而这时第一种方法太慢了。

这是对比度修复的代码
（请注意，我使用了浅灰色而不是完全白色作为颜色;-））

private void HighLightDirtyPixelsLockBits(){
    //http://msdn.microsoft.com/en-us/library/system.drawing.bitmap.lockbits(v=vs.110).aspx 
    //http://stackoverflow.com/questions/6251599/how-to-convert-pixel-formats-from-32bpprgb-to-16bpp-grayscale-in-c-sharp
    //http://stackoverflow.com/questions/22257878/set-white-and-black-pixel-depends-on-brightness
    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);
}

在这里，我们遇到了前面提到的CreateNonIndexedImage函数无法解决的问题。问题是HighLightDirtyPixelsLockBits在某些图像上崩溃了。崩溃的原因是PixelFormat。某些图像具有Format24bppRgb格式，而某些图像具有Format8bppIndexed格式。 CreateNonIndexedImage函数确保所有图像在内存中具有相同的PixelFormat。请注意，在for循环中我使用了i += 4，因为位图内存是ARGB（32位）的。如果您有其他格式，则可能需要更改此设置。例如，24位图像（RGB）只需要i+= 3。

这很简单

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;            
}

HighLightDirtyPixelsLockBits函数在一瞬间将图像设置为2色模式，而GetPixel/SetPixel方法则需要几秒钟才能完成。事实上，它非常快，我使用滑块来更改阈值，如果结果不够好，我会在更改滑块值时更改对比度。

用户仍然需要手动指示太阳黑子的位置，目前是通过在图像上绘制矩形（Image控件上的mousedown、mousemove、mouseup事件），并使用GetPixel/SetPixel方法（遍历绘制矩形的坐标而不是整个图像）。通过这些操作，我们可以突出显示属于太阳黑子的“脏”像素。

绘制该矩形证明是一件令人头疼的事。令人头疼的原因是您需要将鼠标坐标转换为图像上的坐标，否则您的矩形将不会出现在您期望的位置。我不确定我是否以最好的方式解决了这个问题，但它确实有效。

MouseDown事件处理程序

private void image_selected_MouseLeftButtonDown(object sender, System.Windows.Input.MouseButtonEventArgs e) {
    if(drawstate == DRAWSTATE.RECTANGLE){
        //get the value of the 0,0 point of the image on the screen.
        gaugepoint = image_selected.PointToScreen(new System.Windows.Point(0,0));
        //offset the clicked point with the gaugepoint and correct the X,Y factor.
        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);
    }                                                //end if
}

如您所见，鼠标坐标已通过校正因子进行转换，并将一个Rectangle形状添加到包含Image对象的canvas对象中。 sv变量是滚动条位置的校正。

稍后将在本文中讨论DRAWSTATE条件语句。

MouseMove处理程序

private void image_selected_MouseMove(object sender, System.Windows.Input.MouseEventArgs e) {
    if(drawstate == DRAWSTATE.RECTANGLE){
        //Here we just correct the Rectangle we created in image_selected_MouseLeftButtonUp()
        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;
        }                                            //end if
    }                                                //end if
}

这只更新矩形的大小。

最后是MouseUp事件处理程序

private void image_selected_MouseLeftButtonUp(object sender, System.Windows.Input.MouseButtonEventArgs e) {
    if(drawstate == DRAWSTATE.RECTANGLE){
        //offset the clicked point with the gaugepoint and correct the Y factor.
        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));
        //horizontal/vertical offset is used for correction of the scroll bars.
        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));
        //add the rectangle to the controls. (you can do that in WPF ;-))
        ((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);
        //add the rectangle to the Square list.
        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);

        //highlight the pixels considered dirty.
        HighlightDirtyPixels(bmprect, square);
    }                                                //end if
    else{
        if(drawstate == DRAWSTATE.POLYLINE){
            //offset the clicked point with the gaugepoint and correct the Y factor.
            double x = e.GetPosition((System.Windows.Controls.Image)sender).X + gaugepoint.X;
            double y = e.GetPosition((System.Windows.Controls.Image)sender).Y + gaugepoint.Y;// + correction.Y;
            if(tmppolygon == null){
                tmppolygon = new System.Windows.Shapes.Polygon();
                c.Children.Add(tmppolygon);
        }                                        //end if
            tmppolygon.Points.Add(new System.Windows.Point(x, y));
            tmppolygon.Stroke = new SolidColorBrush(System.Windows.Media.Colors.Red);
        }                                            //end if
    }                                                //end else
}

请注意，我在这里将矩形显示为绿色（最初是红色），这只是为了方便用户，让他知道他的矩形已“准备就绪”。（也出于调试目的，请参见下一段）

这里有一个问题是，MouseUp事件大多数时候没有触发，结果我得到了多个（红色）矩形，这些矩形没有按预期工作。原因是我们在mousemove期间更新绘制的矩形，当您释放鼠标按钮时，您位于矩形上方，而不是image_selected控件上方。因此，我添加了

private void newrect_MouseLeftButtonUp(object sender, System.Windows.Input.MouseButtonEventArgs e) {
    image_selected_MouseLeftButtonUp(image_selected, e);
}

您还会注意到有一个drawstate polyline，这是为了解决稍后将提到的一个问题。

这里是使用GetPixel/SetPixel方法的函数

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;
        //Also check that the pixels are INSIDE the circle!
        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;
                }                                    //end if
                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;
                    }                                //end if
                    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;
                    }                                //end if
                }                                    //end else
                pixelindex++;
                square.Pixels.Add(p);
            }                                        //end for
        }                                            //end for
        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);
    }                                                //end try
    catch(Exception ex){
        WriteInfo(msg + "\tError: " + ex.Message);
    }                                                //end catch
    bmp_contrastfix = new System.Drawing.Bitmap(bmp);
    SetModifiedImage(bmp, false);
}

这需要一些解释

bmp_contrastfix变量是对比度修复后保存图像的对象。这里还使用了不同的对象：

sun (BusinessObjects.Objects.Sun)
ssa (BusinessObjects.Objects.SunspotArea)
p (BusinessObjects.Shapes.Pixel)

它们用于在整个过程中计算不同的属性，类定义在下面的“参考文献”部分。

那么，让我们花点时间反思一下。我们将图像加载到内存中，修复了其像素格式以使用lockbits/unlockbits方法进行对比度修复，然后我们可以在图像上绘制一个矩形，在较小的区域内使用GetPixel/SetPixel方法来突出显示太阳黑子（组）。在此过程中，我们修复了几个问题：

确定了“脏”像素的定义
解决了PixelFormat问题。
修复了一个非常恼人和微妙的bug，即在MouseUp事件期间鼠标位于矩形而不是图像上。

如何计算它们的表面积？

所有先前的工作都是为了准确地完成这一项。计算太阳黑子（组）的面积。所以，让我们开始计算绘制矩形内的“脏”像素数量，看看圆（=太阳表面）中有多少像素，我们就算出来了！太棒了！

等等，...圆？呃，哪里有那个圆？是的，您可以在图像上清楚地看到它，但您的程序并不知道圆的边缘是否是太阳黑子的一部分还是其他什么东西...

欢迎来到

问题1 - 检测定义太阳边缘的圆

这出奇地简单，. . . 如果你知道怎么做的话。我尝试了两个库：AForge.Net和Emgu（OpenCV的封装器）。 Emgu的方法明显更好，尽管我必须提供一些手动方法来移动和调整检测到的圆的大小，使其“精确”地拟合绘图。

这是检测方法

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);
    }                                                //end if
    else{
        img = new Image<Bgr,byte>(bmp_ori);
    }                                                //end else

    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));

    //http://www.emgu.com/wiki/files/2.4.2/document/Index.html
    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);
        }                                            //end if
    }                                                //end for
    SetModifiedImage(img.ToBitmap(), false);
}

大部分是Emgu特有的对象，更多信息请参考其页面：http://www.emgu.com/wiki/index.php/Main_Page 并且可以在互联网上找到许多示例。

使用HoughCircles方法可以在图像中找到圆。是所有圆，但您可以稍微调整属性来清除大部分。（例如，圆的最小尺寸，...）它还能找到部分在图像外的圆，所以我使用一个IsCircleInDrawing方法，它看起来像这样：

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;
    }                                                //end if
    return indrawing;
}

这只包括完全在图像边界内的圆。

实际上，它几乎总是能找到正确的圆。在极少数情况下找不到时，您可以更改对比度修复并再次尝试检测。

对比度修复和形状检测看起来是这样的

Contrast Fix And Shape Detection

有了这些，您就可以获得太阳的中心点及其半径。您需要这些才能计算表面积。有趣的部分从这里开始。图纸是一个圆，因此逻辑上的想法是计算太阳黑子面积覆盖圆盘的百分比。既然您知道太阳的直径（http://en.wikipedia.org/wiki/Sun）和像素半径，您就可以通过累加所有脏像素来计算平方米，对吧？

错误！

问题2 - 太阳是一个球体

您在那里计算的表面积是一个投影表面积。太阳是一个球体而不是一个圆盘，因此表面积随着与圆盘中心距离的增加而“增加”。这部分在互联网上不容易找到。

这是代码部分

public double CalculatePixelSurface(BusinessObjects.Shapes.Square square, BusinessObjects.Objects.Sun sun, UNIT unit){
    double surface = 0;
    double X = 0;
    double Y = 0;
    //Formula: (1/N^2) * (1/(2*PI*cos(alpha)) = (1/N^2) * (1/(2*PI*SQRT(1-x^2-y^2)) (N = # pixels of the radius)
    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));            // = SQRT(1-x²-y²)
            cosalpha = (cosalpha > 16 ? 16 : cosalpha);                                //if it is too close to the limb, take a treshold.
            surface += ((1/Math.Pow(sun.RadiusInPixels, 2)) * (1/(2*Math.PI*cosalpha)));
        }                                            //end if
    }                                                //end for
    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;
    }                                                //end switch
    
    return surface;
}

这里我们遍历一个包含太阳黑子的Square对象，并为每个“脏”像素计算其去投影表面积。一旦您有了公式，该方法就不是很困难。然而，得出那个公式并不容易。此外，我取正方形的中心，并计算它相对于中心的角度，以检查它是否与我选择的大致正确。例如。如果我取一个中心附近的方形，最终得到78度，那是不对的。通常，面积表示为太阳半球百万分之一（Msh）的圆盘的比例。

数学。

原始公式看起来是这样的：

Formula

对于非数学家（比如我）来说，这是无用的，所以我采访了一位同事，他向我解释了。

用更编程的方式来看，这看起来像

(1/N^2) * (1/(2*PI*cos(alpha)) = (1/N^2) * (1/(2*PI*SQRT(1-x^2-y^2))

如何得到它？

您需要做的第一件事是将球体缩放到1,1。（这意味着，例如，450像素的半径=1）。

N = 半径的像素数。

我们需要这个缩放是因为它简化了查找公式的过程。例如。如果一个圆的表面积= Pi.r²，您可以说对于半径为1的圆，表面积只是Pi（r² = 1² = 1）。

Scaled Sphere

太阳黑子面积（在圆上）是所有像素的总和。要在缩放后的圆上获得相同的面积，您得到：

像素总和 / N²。

作为圆盘的比例，它变成：

像素总和 / (Pi*N²)

前面的公式很简单，但仍然是二维的。我们不是在二维工作，而是在三维工作，球体的表面积公式是：

4*Pi*r²

太好了。我们有r = 1，所以它简化为4*Pi，但也要注意我们不是在处理一个完整的球体，只是半个。所以我们看到的图画中的圆有一个（缩放的）表面积：

2*Pi

如何获得我们球体的一部分表面积？

为此，您需要知道太阳黑子相对于中心的距离。用通俗的话来说，圆盘（圆）中心的像素的表面积是1像素（也就是说，该像素在二维和三维中具有相等的“表面值”）。您的像素离中心越远，它就越被“投影”，因此需要进行更多的校正。为了获得这个“去投影”值，您需要知道像素/太阳中心线与Z轴的角度（alpha）。是的，Z轴垂直于（或进入）页面。我们目前还没有关于这个第三轴（以及图纸的3D特性）的任何信息。

它看起来会像这样：

Scaled Sphere (3D)

（请注意，轴已与之前的（二维）图纸发生了变化！）

红十字是我们想要计算的点。线C（实线）是从中心到位于球体上的像素。 A是投影像素到XY平面的线。 B是从该投影像素到中心的距离。 alpha是Z轴与直线C之间的角度。

如果您还记得学校里学的

cos(alpha) = B / C

其中B是直角三角形的邻边，C是斜边。（http://en.wikipedia.org/wiki/Trigonometric_functions）

我们知道C，这是球体的半径。缩放后为1。因此，我们只需要找到B的值。

在我们的例子中，公式中的B是线A。这只是Z轴上的值。由于我们只有x,y信息，我们需要找到z值。

在球体中，有：

X² + Y² + Z² = r²

由于半径为1，我们有：

X² + Y² + Z² = 1

Z² = 1 - X² - Y²

Z = SQRT(1 - X² - Y²)

因此，alpha的余弦变为：

cos(alpha) = SQRT(1 - X² - Y²) / 1

因此，太阳黑子（组）在球体上的表面积比例为：

sumofpixels/N² * 1/2*Pi*cos(alpha)

sumofpixels/N² * 1/2*Pi*SQRT(1-X²-Y²)

所以我们在代码中为每个像素所做的是：

1/N² * 1/2*Pi*SQRT(1-X²-Y²)

并将它们（每个像素的表面值）累加起来得到总数。

在代码中，这变成了（CalculatePixelSurface函数）：

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));            // = SQRT(1-x^2-y^2)
        cosalpha = (cosalpha > 16 ? 16 : cosalpha);                                //if it is too close to the limb, take a treshold.
        surface += ((1/Math.Pow(sun.RadiusInPixels, 2)) * (1/(2*Math.PI*cosalpha)));
    }                                        //end if
}                                                //end for

请注意，对余弦值进行了检查（cosalpha = (cosalpha > 16 ? 16 : cosalpha)，这是因为当像素接近东西边缘时，Area AtDiskCenter将趋于无穷大。为了避免这种情况，我们使用Paul Higgins的方法将校正因子（即求和项）限制为16，这对应于86度左右的LOS角度。

更新 2016/03/14：经过一些比较，我们注意到结果中存在一个错误，这导致了

cosalpha = (cosalpha > 16 ? 16 : cosalpha);

我们将其更改为：

double tresholdradians = (2*Math.PI/360.0)*86.0;
cosalpha = (cosalpha < Math.Cos(tresholdradians) ? Math.Cos(tresholdradians) : cosalpha);

这解决了问题。

如何去除不需要的特征？

1.点在多边形内算法

使用矩形会限制您的选择方法。如果太阳黑子组位于圆的边缘（您也选择了边缘，并且那些是“脏”像素）。您还有不需要的特征，如线条、日珥、十字星，...有时会混入选择中。

我们可以移除圆外的任何东西（在HightLightDirtyPixels方法中）

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;

基本上，我们计算所选矩形中的每个像素是否在太阳的半径范围内。

这仍然给我们留下了其他不需要的特征，因此我们研究了使用多边形代替矩形。用户可以手动定义一个多边形并计算该多边形内的脏像素。检测边界矩形以缩短循环周期很容易。在我的原型中，我使用了一个按钮在多边形和矩形之间切换。这是代码：

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;
    }                                                //end if
    else{
        if(    tmpbtn.Content.ToString() == "Polyline On"){
            tmpbtn.Content = "Polyline Off";

            //gets bounding rectangle
            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);
            //Make it slightly larger in width, because of offset effects ...
            r.X -= 20;
            r.Width = += 20;  //make it larger, not move it.
            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;
        }                                            //end if
    }                                                //end else
}

请注意，我们将边界矩形左侧做得稍宽一些。这是因为算法总是在多边形左侧跳过几个像素（见下图）。将矩形做得稍宽一些可以解决这个问题。

因此，下一个难题是确定一个像素是否在一个多边形内。这里有一篇关于此的文章：http://alienryderflex.com/polygon/

一旦找到代码，它就很简单：

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);
        }                                            //end if
        j = i;
    }                                                //end for
    return pinp;
}

int sides =polygon.Points.Count 如果边交叉，则不完全正确，但在我们的例子中，用户应该定义一个围绕某物的多边形，因此边不应该交叉。

这是一张最终的图像，其中太阳黑子是用多边形选择的：

Sunspot Highlight

（注意：多边形问题右侧的小黑尖通过将边界矩形做得稍宽一些来解决）

2.图像遮罩方法

寻找多边形内像素的另一种方法是复制图像并在此图像上绘制一个填充的多边形。（Graphics.FillPolygon方法）。通过用奇怪的颜色值（比如青色）填充多边形，您可以遍历相同的边界矩形，并检查颜色是否为青色（多边形内部）或其他颜色（多边形外部）。它是这样的...：

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;
        //Also check that the pixels are INSIDE the circle!
        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;
                }                                    //end if
                else{
                    if(insidepolygon){
                        bmp.SetPixel(x, y, System.Drawing.Color.White);
                        p.PixelColor = System.Drawing.Color.White;
                        msg = "white\t" + "\tx:" + x + "\ty:" + y;
                    }                                //end if
                }                                    //end else
                pixelindex++;
                square.Pixels.Add(p);
            }                                        //end for
        }                                            //end for
        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);
    }                                                //end try
    catch(Exception ex){
        WriteInfo(msg + "\tError: " + ex.Message);
    }                                                //end catch
    bmp_contrastfix = new System.Drawing.Bitmap(bmp);
    SetModifiedImage(bmp, false);
}

请注意，这里我们将bmp对象复制到bmp_mask中，并使用Grapics对象将填充的多边形绘制到图像上（而不是绘制在图像之上）。我们可以遍历边界矩形的坐标，并检查遮罩的颜色，以了解原始图像上的像素是否位于遮罩上绘制的多边形内部。

PointInPolygon2方法与PointInPolygon方法相比得到了简化。

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){    //aqua
        pinp = true;
    }
}

“点在多边形内算法”和“图像遮罩方法”似乎都同样有效，尽管第二种方法占用的内存更多。

关注点

当我开始这个项目时，我没有图像处理经验，我是一边学一边做的。
在做像这样的复杂事情时，始终先从原型开始，直到您满意为止。然后，也只有在那时，才将其放入实际应用程序中。
您需要明白，对于软件来说，当我们进行对比度修复时，所有“白色部分”都被视为“脏”的。这包括圆盘外部的任何区域、圆盘边缘本身、所有太阳黑子以及所有不需要的特征。
这仍然是一个原型。我需要许多人的帮助（包括在这里Codeproject上）才能使其正常工作。但是，代码中仍然可能存在错误。（例如，在选择太阳黑子组的polyline方法中仍然有一个错误）
一次处理一件事。正如您所见，一些更改对其之前编写的代码产生了影响。例如。 CreateNonIndexedImage函数是在我遇到“对比度修复”中的错误之后才出现的，并在图像加载时使用（几乎是我写的第一件事）。因此，请确保一次只处理一件事。即使这样，事情也可能很快变得非常复杂。试图一次做太多事情是灾难的根源。
在此过程中我（我们）学到了什么？
- 图像的PixelFormat可能会扰乱lockbits/unlockbits方法。即使它没有扰乱（引发异常），您仍然需要知道此信息才能进行for循环。
- 鼠标事件按预期工作，但不一定按您预期的那样工作。（这里的例子是mouseup事件在square对象而不是image对象上触发。
- 矩形“选择”是最简单的方法，但不够准确。
- 使用Emgu进行特征检测。我们通过更改搜索属性并检查完全在绘图内的圆来限制结果，从而成功地获得了不需要的（检测到的）特征。
- 大量的数学，因为我们知道太阳是一个球体（而不是一个圆盘）。我非常享受向最终公式推导过程中的推理的优雅。

有趣的事实

在我们的数据集中大约有23000张图纸。
最古老的图纸可追溯到1940年3月5日。
至今仍有图纸被绘制（不只是我们研究所的）。

历史

2014年9月：第一个版本。

参考文献

使用的对象

以下是我在代码中使用的对象。我只是添加了它们作为参考。

类 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;
            }                                                //end set
            get{
                return radiusinmeter;
            }                                                //end get
        }
        
        public int RadiusInPixels{
            set{
                radiusinpixels = value;
            }                                                //end set
            get{
                return radiusinpixels;
            }                                                //end get
        }

        public double MeterPerPixel{
            get{
                return (double)radiusinmeter / (double)radiusinpixels;
            }                                                //end get
        }
    
        public double SquareMeterPerPixel{
            get{
                double surfacep = TotalSolarSurface(UNIT.Pixels);
                double surfacem = TotalSolarSurface(UNIT.Meter);
                return surfacem / surfacep;
            }                                                //end get
        }

        public double TotalSolarSurface(UNIT unit){
            //Surface of a sphere = 4 * PI * r²
            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;
            }                                                //end switch
            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;
            }                                                //end switch
            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();
        }

    }
}

类 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,            //fraction of the total surface
            Meter,    
            Msd,                //millionth of the solar disc
            Msh,                //millionth of the solar hemisphere
            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;
            //Formula: (1/N²) * (1/(2*PI*cos(alpha)) = (1/N²) * (1/(2*PI*SQRT(1-x²-y²)) (N = # pixels of the radius)
            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));            // = SQRT(1-x²-y²)
                    //cosalpha = (cosalpha > 16 ? 16 : cosalpha);                                //if it is too close to the limb, take a treshold.
                    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)));
                }                                            //end if
            }                                                //end for
            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;
            }                                                //end switch
            
            return surface;
        }

    }
}

类 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; }
    }
}

类 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; }
    }
}

类 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; }
    }
}

结论

这是我第一次创建如此复杂的原型，正是因为经验（以及一些帮助），我才成功了。最好的决定是不要开始在现有应用程序中编码。对于更有经验的程序员来说，这可能很明显，但对于初学者来说，这是一个很好的建议。