LIBMAS/DedicatedAnalysis/GrainSizeStandard/IntegrationClass/OLDcellsegm.cs

using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using System.Runtime.InteropServices;
using System.Text;
using System.Threading.Tasks;
using OpenCvSharp;
using Point = System.Drawing.Point;
using Size = OpenCvSharp.Size;
using ZhaoKThinAndConnectionDLL;

namespace PaintDotNet.DedicatedAnalysis.GrainSizeStandard.IntegrationClass
{
    internal class InParameter
    {
        public Mat MGryImg;

        public Mat MHesOutImage;

        public Mat MValley;

        public int iDimA;

        public int iDimB;
    }
    public class NewSeg
    {
        private Mat m_MOrgImg;

        private Mat m_MOutImg;

        private double[] m_adOutParam;

        private InParameter m_InParam;

        public int DoSeg(Mat _MOrgImg, Mat _MOutImg, double[] _adOutParam, Vec4b value)
        {
            PrepareParameter();
            m_MOrgImg = _MOrgImg;
            m_MOutImg = _MOutImg;
            m_adOutParam = _adOutParam;
            if (3 == m_MOrgImg.Channels())
            {
                Cv2.CvtColor(m_MOrgImg, m_InParam.MGryImg, ColorConversionCodes.RGB2GRAY);
            }
            else
            {
                if (1 != m_MOrgImg.Channels())
                {
                    return -1;
                }

                m_InParam.MGryImg = m_MOrgImg.Clone();
            }

            ImgHessian();
            FindValley();
            FindShed(value);
            return 1;
        }

        private int FindShed(Vec4b value)
        {
            Mat mat = new Mat(m_InParam.MValley.Size(), 4);
            Cv2.FindContours(m_InParam.MValley, out OpenCvSharp.Point[][] contours, out HierarchyIndex[] hierarchy, RetrievalModes.CComp, ContourApproximationModes.ApproxSimple);
            for (int i = 0; i < hierarchy.Length; i++)
            {
                Cv2.DrawContours(mat, contours, i, Scalar.All(i + 1), 1, LineTypes.Link8, hierarchy);
            }

            Cv2.Watershed(m_MOrgImg, mat);
            int rows = mat.Rows;
            int cols = mat.Cols;
            int[] array = new int[2];
            //Vec3b value = new Vec3b(0, 0, byte.MaxValue);
            array[0] = 0;
            while (array[0] < rows)
            {
                array[1] = 0;
                while (array[1] < cols)
                {
                    if (0 >= mat.At<int>(array[0], array[1]))
                    {
                        m_MOutImg.Set(array, value);
                    }

                    array[1]++;
                }

                array[0]++;
            }

            return 1;
        }

        private int FindValley()
        {
            Mat structuringElement = Cv2.GetStructuringElement(MorphShapes.Rect, new Size(m_InParam.iDimA, m_InParam.iDimA), new OpenCvSharp.Point(-1, -1));
            Cv2.MorphologyEx(m_InParam.MHesOutImage, m_InParam.MHesOutImage, MorphTypes.Close, (InputArray)structuringElement);
            Cv2.BitwiseNot(m_InParam.MHesOutImage, m_InParam.MValley);
            Size size = m_InParam.MValley.Size();
            Mat mat = (Mat)Mat.Zeros(size.Height + 2, size.Width + 2, m_InParam.MValley.Type());
            m_InParam.MValley.CopyTo(mat.SubMat(new Range(1, size.Height + 1), new Range(1, size.Width + 1)));
            Cv2.FloodFill(mat, new OpenCvSharp.Point(0, 0), new Scalar(255.0));
            Mat mat2 = new Mat();
            mat.SubMat(new Range(1, size.Height + 1), new Range(1, size.Width + 1)).CopyTo(mat2);
            Cv2.BitwiseNot(mat2, mat2);
            m_InParam.MValley = (Mat)(m_InParam.MValley | mat2);
            return 1;
        }

        private int PrepareParameter()
        {
            m_InParam = new InParameter();
            m_InParam.MGryImg = new Mat();
            m_InParam.MHesOutImage = new Mat();
            m_InParam.MValley = new Mat();
            m_InParam.iDimA = 3;
            m_InParam.iDimB = 1;
            return 1;
        }

        private int ImgHessian()
        {
            Mat mat = new Mat();
            m_InParam.MGryImg.ConvertTo(mat, 5);
            int rows = mat.Rows;
            int cols = mat.Cols;
            int num = 5;
            double num2 = 1.2;
            Mat mat2 = new Mat(2 * num + 1, 2 * num + 1, MatType.CV_32FC1, Scalar.All(0.0));
            Mat mat3 = new Mat(2 * num + 1, 2 * num + 1, MatType.CV_32FC1, Scalar.All(0.0));
            Mat mat4 = new Mat(2 * num + 1, 2 * num + 1, MatType.CV_32FC1, Scalar.All(0.0));
            for (int i = -num; i <= num; i++)
            {
                for (int j = -num; j <= num; j++)
                {
                    double num3 = 1.0 - (double)(i * i) / (num2 * num2);
                    double num4 = Math.Exp((double)(-1 * (i * i + j * j)) / (2.0 * num2 * num2));
                    double num5 = -1.0 / (Math.PI * 2.0 * Math.Pow(num2, 4.0));
                    double num6 = num3 * num4 * num5;
                    float value = (float)num6;
                    mat2.Set(i + num, j + num, value);
                    mat4.Set(i + num, j + num, (float)((1.0 - (double)(j * j) / (num2 * num2)) * Math.Exp((double)(-1 * (i * i + j * j)) / (2.0 * num2 * num2)) * (-1.0 / (Math.PI * 2.0 * Math.Pow(num2, 4.0)))));
                    mat3.Set(i + num, j + num, (float)((double)(i * j) * Math.Exp((double)(-1 * (i * i + j * j)) / (2.0 * num2 * num2)) * (1.0 / (Math.PI * 2.0 * Math.Pow(num2, 6.0)))));
                }
            }

            Mat mat5 = new Mat(rows, cols, MatType.CV_32FC1, Scalar.All(0.0));
            Mat mat6 = new Mat(rows, cols, MatType.CV_32FC1, Scalar.All(0.0));
            Mat mat7 = new Mat(rows, cols, MatType.CV_32FC1, Scalar.All(0.0));
            Cv2.Filter2D(mat, mat5, mat5.Depth(), mat2);
            Cv2.Filter2D(mat, mat6, mat6.Depth(), mat4);
            Cv2.Filter2D(mat, mat7, mat7.Depth(), mat3);
            byte[,] array = new byte[rows, cols];
            float[] array2 = new float[rows * cols];
            float[] array3 = new float[rows * cols];
            float[] array4 = new float[rows * cols];
            Marshal.Copy(mat5.Data, array2, 0, rows * cols);
            Marshal.Copy(mat6.Data, array3, 0, rows * cols);
            Marshal.Copy(mat7.Data, array4, 0, rows * cols);
            for (int k = 0; k < rows; k++)
            {
                for (int l = 0; l < cols; l++)
                {
                    float num7 = array2[k * cols + l];
                    float num8 = array3[k * cols + l];
                    float num9 = array4[k * cols + l];
                    double num10 = (num7 + num8) * (num7 + num8) - 4f * (num7 * num8 - num9 * num9);
                    if (!(num10 < 0.0))
                    {
                        num10 = Math.Sqrt(num10);
                        double num11 = ((double)(num7 + num8) + num10) / 2.0;
                        double value2 = ((double)(num7 + num8) - num10) / 2.0;
                        if (num11 > 0.0 && Math.Abs(num11) > 1.0 + Math.Abs(value2))
                        {
                            array[k, l] = byte.MaxValue;
                        }
                    }
                }
            }

            m_InParam.MHesOutImage = new Mat(rows, cols, MatType.CV_8UC1, array, 0L);
            return 0;
        }
    }

    /////
    /////
    class OUTsegmsurf
    {
        public Mat cellbw;
        public Mat wat;
        public Mat imsegm;
        public Mat minima;
        public Mat minimacell;
        //kkkkk public NEW_TYPE_SegInf info;
    }
    class In5segmsurf
    {
        public string smoothim_method;
        public int filterridges;
        public double gaussian_stdev;
        public string LEVEL;
        public double classifycells_convexarea;
        public double classifycells_convexperim;
        public Mat OrgImg;
        public In5segmsurf Clone()
        {
            return this.MemberwiseClone() as In5segmsurf;
        }
    }
    class In6segmsurf
    {
        public string prmfile;//Path to parameter file
        public Mat minima;//All markers, for instance from manual drawings
        public Mat minimacell;//; for cells, for instance from manual drawings
        public Mat imnucl;//Nucleus image for defining markers
    }
    class PMergefragments//区域合并
    {
        public int optlog;
        public int optint;
        public double intint;
        public double conv;
        public PMergefragments Clone()
        {
            return this.MemberwiseClone() as PMergefragments;
        }
    }
    class PSmoothim//平滑
    {
        public string method;
        public double[] h;
        public int gpu;
        public int dim;
        public int ced_maxniter;
        public PSmoothim Clone()
        {
            return this.MemberwiseClone() as PSmoothim;
        }
    }
    class PClassifycells//区域分类
    {
        public double convexarea;
        public double convexperim;
        public double[] h;
        public double minvolfull;
        public double maxvolfull;
        public string method;

        public PClassifycells Clone()
        {
            return this.MemberwiseClone() as PClassifycells;
        }

    }
    class In4smoothim
    {
        public int filterridges;
        public string LEVEL;

        public double gaussian_stdev;

        public double[] getminima_h;
        public double getminima_minvolfull;
        public double getminima_maxvolfull;
        public double getminima_level;
        //public string getminima_method;

        public PMergefragments mergefragments;
        public PSmoothim smoothim;
        public PClassifycells classifycells;

        public int[] dim;
        public double minvol;
        public double minvolvox;
        public double maxvol;
        public double maxvolvox;

        public double watminvolfull;
        public double watmaxvolfull;
        public double minvolfull;
        public double maxvolfull;

        public double[] h;
        public double just;
        public int gpu;
        public int illum;
        public double illumdiameter;
        public int merge;
        public int iInterpolationMethod;//zzzzz new, 1 for Cubic spline interpolation, 2 for Bilinear interpolation, 3 for Nearest neighbor interpolation
    }

    class OUTgetminima
    {
        public Mat minimacell;
        public Mat minima;
    }
    class OUTclassifycells
    {
        public Mat cellbw;
        public PClassifycells classifycells;
    }
    class OUTmaxfilt3
    {
        public List<float>[,] maxim;
        public float[,] sumim;
        public int numpoints;
        public float[][][] minim;
    }
    class MyFilter
    {
        public double[] x;
        public double[] y;
        public double[] z;
    }
    class OLDcellsegm
    {
        public const double m_PI = 3.1415926535897931;
        void mf_mergeinputpar(In4smoothim prm, In5segmsurf prmin)
        {
            prm.smoothim.method = prmin.smoothim_method;
            prm.filterridges = prmin.filterridges;
            prm.gaussian_stdev = prmin.gaussian_stdev;
            prm.LEVEL = (string)prmin.LEVEL.Clone();
            prm.classifycells.convexarea = prmin.classifycells_convexarea;
            prm.classifycells.convexperim = prmin.classifycells_convexperim;
        }
        void mf_cellsize(double minvol, double maxvol, In4smoothim prm, int ThreeDimNO)
        {
            if (0 == prm.h[0] * prm.h[1] * prm.h[2])
            {
                Console.WriteLine("There is zero voxel size");
            }
            double minr = Math.Pow((3 * minvol) / (4 * m_PI), 1 / 3);
            //maxr = ((3*maxvol)/(4*pi))^(1/3);
            //maxz = maxr/h(3);
            double minz = minr / prm.h[2];
            if (ThreeDimNO < minz)
            {
                //shrinink a sphere to a circle
                double r = Math.Pow((3 * maxvol) / (4 * m_PI), 1 / 3);
                maxvol = Math.Pow(ThreeDimNO, prm.just) * m_PI * r * r;
                r = Math.Pow((3 * minvol) / (4 * m_PI), 1 / 3);
                minvol = Math.Pow(ThreeDimNO, prm.just) * m_PI * r * r;
            }
            double voxelvol = prm.h[0] * prm.h[1] * prm.h[2];
            double minvolvox = Math.Round(minvol / voxelvol, MidpointRounding.AwayFromZero);
            double maxvolvox = Math.Round(maxvol / voxelvol, MidpointRounding.AwayFromZero);
            prm.minvol = minvol;
            prm.minvolvox = minvolvox;
            prm.maxvol = maxvol;
            prm.maxvolvox = maxvolvox;
        }

        public OUTsegmsurf mf_segmsurf_progress_auto(Mat GryImg, double minvol, double maxvol,
            string[] INnamehere, In5segmsurf IN5varhere, In6segmsurf IN6varhere, Color phaseColor)
        {
            long start = Cv2.GetTickCount();
            OUTsegmsurf rtn = new OUTsegmsurf();
            //ttttt Cv2.FindContours(watermark, out contour, out hier, RetrievalModes.CComp, ContourApproximationModes.ApproxSimple, null);
            rtn.wat = new Mat(GryImg.Size(), MatType.CV_8UC4, new Scalar(255, 0, 0, 0));
            //Mat OutImg = new Mat(GryImg.Rows, GryImg.Cols, MatType.CV_8UC3, Scalar.All(0));
            NewSeg NSeg = new NewSeg();
            //转为为灰度?
            if (/*3 == GryImg.Channels() || */4 == GryImg.Channels())
            {
                Cv2.CvtColor(GryImg, GryImg, ColorConversionCodes.BGRA2BGR);
            }
            else if (1 == GryImg.Channels())
            {
                Cv2.CvtColor(GryImg, GryImg, ColorConversionCodes.GRAY2BGR);
            }
            //Cv2.ImShow("oriImg", GryImg);
            //Cv2.WaitKey();
            double time_v2 = (Cv2.GetTickCount() - start) / Cv2.GetTickFrequency();
            System.Console.WriteLine("晶界重现_V2_执行时间1：" + time_v2);
            NSeg.DoSeg(GryImg.Clone()/*OrgImg*/, rtn.wat, null, new Vec4b(phaseColor.B, phaseColor.G, phaseColor.R, 255));
            ////ttttt Cv2.DrawContours(rtn.wat, contour, -1, Scalar.Red);
            //int Rows = OutImg.Rows;
            //int Cols = OutImg/*watermark*/.Cols;
            //int[] point = new int[2];
            //Vec4b redcol;
            //if (phaseColor != null)
            //{
            //    redcol = new Vec4b(phaseColor.B, phaseColor.G, phaseColor.R, 255);
            //}
            //else
            //    redcol = new Vec4b(0, 0, 255/*255*/, 255);
            //for (point[0] = 0; point[0] < Rows; point[0]++)
            //{
            //    for (point[1] = 0; point[1] < Cols; point[1]++)
            //    {
            //        if (0 >= watermark.At<int>(point[0], point[1]))
            //        {
            //            rtn.wat.Set<Vec4b>(point[0], point[1], redcol);//###
            //        }
            //    }
            //}
            time_v2 = (Cv2.GetTickCount() - start) / Cv2.GetTickFrequency();
            System.Console.WriteLine("晶界重现_V2_执行时间2：" + time_v2);
            ////Cv2.ImShow("contour", rtn.wat);//kkkkk
            ////Cv2.WaitKey(0);//kkkkk
            //control.Cursor = cursor;// System.Windows.Forms.Cursors.Default;

            return rtn;

            //System.Windows.Forms.Cursor cursor = control.Cursor;//.GetHashCode();//.CopyHandle();
            //control.Cursor = System.Windows.Forms.Cursors.WaitCursor;
            string mfilename = System.Diagnostics.Process.GetCurrentProcess().MainModule.FileName;
            Console.WriteLine("This is" + mfilename + "for segmentation of cells");
            //image for segmentation
            Mat imsegm = GryImg.Clone();
            //cell volumes in cubic microns
            double minvolfull = 1000 * minvol;
            double maxvolfull = 1000 * maxvol;
            Mat imsegmini = imsegm.Clone();
            Mat minimacell = null;
            Mat minima = null;
            Mat imnucl = null;
            string prmfile = null;
            /*
            //for visualilzation 
            int vis = 0;
            int visfinal = 0;
            */

            In4smoothim prm = new In4smoothim();
            prm.mergefragments = new PMergefragments();
            prm.smoothim = new PSmoothim();
            prm.classifycells = new PClassifycells();
            /*Standard parameters*/
            //voxel size
            prm.h = new double[3] { 0.5, 0.5, 1.5 };
            //adjust for not fully circular shape of 3D cells
            prm.just = 0.9;
            //filter ridges
            prm.filterridges = 0;
            //GPU, requires
            prm.gpu = 0;
            //correct illumination
            prm.illum = 0;
            prm.illumdiameter = 25;
            /*Minima parameters*/
            //level
            prm.getminima_level = 0.30;
            /*
            //method
            prm.getminima_method = "automated";
            */
            /*Smoothing parameters*/
            //directional coherence enhancing diffusion
            prm.smoothim.method = "dirced";
            //2D or 3D smoothing
            prm.smoothim.dim = 2;
            //number of iterations of smoothing in CED
            prm.smoothim.ced_maxniter = 100;
            //dimension image
            int[] dim = new int[3] { imsegm.Rows, imsegm.Cols, 1 };
            /*Classification parameters*/
            prm.classifycells.method = "threshold";
            /*Watershed parameters*/
            //To test intensity of wat lines; merging
            prm.merge = 0;
            //methods for merging
            prm.mergefragments.optlog = 2;
            prm.mergefragments.optint = 2;
            /* The signicance of watershed lines, stronger if larger, done on ridgeim!!
            if the raw image, then use much smaller, around 1
            set to 1.05-1.10
            prm.watthint = 1.06;*/
            prm.mergefragments.intint = 1.3;
            /* The conexity measure for the merging, we should not allow a very much
            smaller convexity after merging
            thconv = 1 if no change in convexity
            thconv > 1 if increase in convexity
            thconv < 1 if decrease in convexity */
            prm.mergefragments.conv = 1.0;

            if (maxvolfull < minvolfull)
            {
                Console.WriteLine(mfilename + ": The upper cell volume is lower than the lower, check the input");
            }
            /* 无用matlab代码       
            % read from parameter file?
            for i = 4 : 2 : nargin
                namehere = varargin{4};
                varhere = varargin{5};
                switch(namehere)
                    case('prmfile')
                        prmfile = varhere;            
                    otherwise
                        % nothing
                end;
            end;

            %
            % Load the parameter file if one is given
            %
            if ~isempty(prmfile)
                isfolder = ~isempty(strfind(prmfile,'/'));
                [a b c] = fileparts(prmfile);
                if isfolder        
                    A = pwd;
                    cd(a);        
                end;    
                msg = sprintf('Reading from parameter file %s',prmfile);
                disp(msg);
                cmd = [b c];
                prmin = eval(cmd);
    
                if isfolder
                    cd(A);
                end;
                % merge the structs by overriding the existing values
                prm = mergestruct(prm,prmin);
            end;
            */

            /*The command line has the highest priority, merge the variables here*/
            In5segmsurf prmin = null;
            for (int i = 0; i < INnamehere.GetLength(0); i++)
            {
                //this variable
                string namehere = INnamehere[i];
                switch (namehere)
                {
                    case "minima":
                        minima = IN6varhere.minima;
                        break;
                    case "prm":
                        prmin = IN5varhere.Clone();
                        mf_mergeinputpar(prm, prmin);
                        break;
                    case "minimacell":
                        minimacell = IN6varhere.minimacell;
                        break;
                    case "imnucleus":
                        imnucl = IN6varhere.imnucl;
                        break;
                    default:
                        Console.WriteLine("Wrong option " + namehere);
                        break;
                }
            }
            prm.dim = new int[2] { imsegm.Rows, imsegm.Cols };

            /*  无用matlab代码
            % if you give in a nucleus image you want to use it!
            if ~isempty(imnucl)
                prm.getminima.method = 'nucleus';
                prm.classifycells.method = 'minimacell';
            end;

            % test for credibility!
            if isequal(prm.getminima.method,'manual')
                if isempty(minima) && isempty(minimacell)
                    error('You must specify either minima eller minimacell')
                end;
            end;
            if isequal(prm.getminima.method,'nucleus')
                if isempty(imnucl)
                    error('You must specify a nucleus image');
                end;
            end;
            */
            /*Computed settings*/
            //get the ADJUSTED 3D cell size, in case its not a full 3D volume. 
            //NB: This is not always linearly scalable, try to avoid using reduced 3D.
            mf_cellsize(minvolfull, maxvolfull, prm, dim[2]);

            //give the voxel sizes in thousand, since all programs are tuned for this
            //prm.classifycells.minvolfull = prm.classifycells.minvolfull;
            //prm.classifycells.maxvolfull = prm.classifycells.maxvolfull;
            prm.getminima_minvolfull = minvolfull / 1000;
            prm.getminima_maxvolfull = maxvolfull / 1000;
            prm.watminvolfull = minvolfull / 1000;
            prm.watmaxvolfull = maxvolfull / 1000;
            prm.minvolfull = minvolfull;
            prm.maxvolfull = maxvolfull;
            prm.classifycells.minvolfull = minvolfull / 1000;
            prm.classifycells.maxvolfull = maxvolfull / 1000;

            Console.WriteLine("Using settings:");
            /* 无用matlab代码
            printstructscreen(prm);
            */

            //pixel size
            prm.getminima_h = (double[])prm.h.Clone();
            prm.classifycells.h = (double[])prm.h.Clone();
            prm.smoothim.h = (double[])prm.h.Clone();

            //gpu for smoothing
            prm.smoothim.gpu = prm.gpu;

            //zzzzz
            prm.iInterpolationMethod = 1;

            /*Cell segmentation start*/
            //Correcting for uneven illumination
            //Smoothing
            PSmoothim prminsmt = prm.smoothim.Clone();
            imsegm = smoothim(imsegm, prm.smoothim.method, "prm", prm);
            //zzzzz Cv2.ImShow("smoothim", imsegm);//kkkkk
            //zzzzz Cv2.WaitKey(0);//kkkkk
            //Ridge enhancement
            if (1 == prm.filterridges)
            {
                Console.WriteLine("Ridge enhancement");
                //kkkkk imsegm = cellsegm.ridgeenhhessian(imsegm,prm.h);
            }
            else
            {
                Console.WriteLine("No ridge enhancement");
            }

            //Finding minima
            OUTgetminima gmout = getminima(imsegm, minimacell, prmin);
            //kkkkk minimacell = gmout.minimacell;
            minima = gmout.minima;

            //3D watershed segmentation
            Console.WriteLine("3D watershed segmentation");

            //Iimpose = imimposemin(imsegm,minima);
            //wat = watershed(Iimpose);

            //vector<vector<Point>> contour;
            Mat watermark = new Mat(imsegm.Size(), MatType.CV_32S);
            OpenCvSharp.Point[][] contour;
            HierarchyIndex[] hier;
            Cv2.FindContours(minima, out contour, out hier, RetrievalModes.CComp, ContourApproximationModes.ApproxSimple, null);
            for (int i = 0; i < hier.Length; i++)
            {
                Cv2.DrawContours(watermark, contour, i, Scalar.All(i + 1), 1, LineTypes.Link8, hier);
            }
            //Mat Temp = new Mat(imsegm.Size(), MatType.CV_8UC3);
            //imsegm.CopyTo(Temp);
            Cv2.Watershed(IN5varhere.OrgImg, watermark);//mmmmm

            //Mat afterWatershed = new Mat();//kkkkk
            //Cv2.ConvertScaleAbs(watermark, afterWatershed);//kkkkk
            //Cv2.ImWrite("D:\\afterWatershed.jpg", afterWatershed);//kkkkk
            //Cv2.ImShow("watermark", afterWatershed);//kkkkk
            //Cv2.WaitKey(0);//kkkkk

            PClassifycells prmclf = prm.classifycells.Clone();

            OUTclassifycells cfout;
            if (null != minimacell)//mmmmm
            {
                cfout = classifycells(watermark, imsegmini, prmclf, minimacell);
            }
            else
            {
                cfout = classifycells(watermark, imsegmini, prmclf);
            }
            Mat cellbw = cfout.cellbw;
            //kkkkk info.classifycells = cfout.classifycells;
            //kkkkk info.prm = prm;

            //9. 绘制轮廓 drawContours
            //Mat maskers = Mat::zeros(imgDist8U.size(), CV_32SC1);
            //for (size_t i=0; i<contour.size(); i++){
            //drawContours(maskers, contour, static_cast<int>(i), Scalar::all(static_cast<int>(i) + 1)); }
            //imshow("maskers", maskers);

            //ttttt Cv2.FindContours(watermark, out contour, out hier, RetrievalModes.CComp, ContourApproximationModes.ApproxSimple, null);
            rtn.wat = new Mat(watermark.Size(), MatType.CV_8UC4, new Scalar(255, 0, 0, 0));
            //ttttt Cv2.DrawContours(rtn.wat, contour, -1, Scalar.Red);
            int Rows = watermark.Rows;
            int Cols = watermark.Cols;
            int[] point = new int[2];
            Vec4b redcol;
            if (phaseColor != null)
            {
                redcol = new Vec4b(phaseColor.B, phaseColor.G, phaseColor.R, 255);
            }
            else
                redcol = new Vec4b(0, 0, 255/*255*/, 255);
            for (point[0] = 0; point[0] < Rows; point[0]++)
            {
                for (point[1] = 0; point[1] < Cols; point[1]++)
                {
                    if (0 >= watermark.At<int>(point[0], point[1]))
                    {
                        rtn.wat.Set<Vec4b>(point[0], point[1], redcol);//###
                    }
                }
            }
            double time = (Cv2.GetTickCount() - start) / Cv2.GetTickFrequency();
            System.Console.WriteLine("晶界重现_执行时间：" + time);
            ////Cv2.ImShow("contour", rtn.wat);//kkkkk
            ////Cv2.WaitKey(0);//kkkkk
            //control.Cursor = cursor;// System.Windows.Forms.Cursors.Default;

            return rtn;
        }
        
        OUTclassifycells classifycells(Mat watermark, Mat imsegmini, PClassifycells prmclf, Mat minimacell = null)
        {
            OUTclassifycells rtn = new OUTclassifycells();
            return rtn;
        }
        OUTgetminima getminima(Mat imsegm, Mat minimacell, In5segmsurf prmin)
        {
            OUTgetminima rtn = new OUTgetminima();
            rtn.minima = null;
            rtn.minimacell = null;
            if (null == minimacell)
            {
                rtn.minima = automated(imsegm, prmin);
            }
            return rtn;
        }
        Mat automated(Mat imsegm, In5segmsurf prmin)
        {
            Mat mat_mean = new Mat(), mat_stddev = new Mat();
            double m, s, d;
            //Cv2.MeanStdDev(imsegm, mat_mean, mat_stddev);
            //m = mat_mean.At<double>(0, 0);
            //s = mat_stddev.At<double>(0, 0);
            //d = 0;//mmmmm
            double th;
            //th = Math.Min(Math.Max(m+d*s, m/2),1.5*m);
            Mat minima = new Mat(), minimacell = new Mat(imsegm.Size(), MatType.CV_8UC3);
            //Cv2.Threshold(imsegm, minimacell, th, 255, ThresholdTypes.Binary);
            int iRad = 39;
            Mat matTempA = new Mat();
            Mat matTempB = new Mat();
            Cv2.MedianBlur(imsegm, matTempA, iRad);
            Cv2.AddWeighted(imsegm, 1, matTempA, -1, 0, matTempB);
            Cv2.MeanStdDev(matTempB, mat_mean, mat_stddev);
            m = mat_mean.At<double>(0, 0);
            s = mat_stddev.At<double>(0, 0);
            d = 0;//mmmmm
            th = m + d * s;
            Cv2.Threshold(matTempB, minimacell, th, 255, ThresholdTypes.Binary);
            int diam;//mmmmm
            Mat kernel;
            diam = 9;
            kernel = Cv2.GetStructuringElement(MorphShapes.Cross, new Size(diam, diam), new OpenCvSharp.Point(-1, -1));
            Cv2.Dilate(minimacell, minimacell, kernel);
            diam = 3;
            kernel = Cv2.GetStructuringElement(MorphShapes.Rect, new Size(diam, diam), new OpenCvSharp.Point(-1, -1));
            Cv2.Erode(minimacell, minimacell, kernel);
            Cv2.MorphologyEx(minimacell, minimacell, MorphTypes.Open, kernel);
            //zzzzz Cv2.ImShow("bianjie", minimacell);
            //zzzzz Cv2.WaitKey(0);
            Cv2.BitwiseNot(minimacell, minima);
            Size bianjieSize = minima.Size();
            Mat Temp = Mat.Zeros(bianjieSize.Height + 2, bianjieSize.Width + 2, minimacell.Type());//延展图像
            minima.CopyTo(Temp.SubMat(new Range(1, bianjieSize.Height + 1), new Range(1, bianjieSize.Width + 1)));
            Cv2.FloodFill(Temp, new OpenCvSharp.Point(0, 0), new Scalar(255));
            Mat cutImg = new Mat();
            Temp.SubMat(new Range(1, bianjieSize.Height + 1), new Range(1, bianjieSize.Width + 1)).CopyTo(cutImg);
            Cv2.BitwiseNot(cutImg, cutImg);
            minima = minima | cutImg;
            diam = 9;
            kernel = Cv2.GetStructuringElement(MorphShapes.Rect, new Size(diam, diam), new OpenCvSharp.Point(-1, -1));
            Cv2.MorphologyEx(minima, minima, MorphTypes.Open, kernel);
            //zzzzz Cv2.ImShow("neibu", minima);//kkkkk
            //zzzzz Cv2.WaitKey(0);//kkkkk
            return minima;
        }
        Mat smoothim(Mat im, string method, string prmflg, In4smoothim prm)
        {
            //diameter of filter
            int iPrmD = 13;
            Mat rst = new Mat();
            Cv2.ConvertScaleAbs(dircohenh(im, iPrmD, prm.h, prm.iInterpolationMethod), rst);
            return rst;
        }
        Mat dircohenh(Mat im, int iPrmD, double[] adPrmH, int iInterpolationMethod)
        {
            int[] dim = new int[3] { im.Rows, im.Cols, 1 };
            if (3 == im.Dims())
            {
                dim[2] = im.Height;
            }
            //number of elements to remove before one takes the Olympic average. Set to
            //1-->3
            int numremele = 2;
            int stepxy = 30;
            int numdeghere = (180 - stepxy) / stepxy + 1;
            int[] deghere = new int[numdeghere];
            deghere[0] = 0;
            for (int i = 1; i < numdeghere; i++)
            {
                deghere[i] = deghere[i - 1] + stepxy;
            }
            double[,] deg2D = new double[numdeghere, 2];
            for (int i = 1; i < numdeghere; i++)
            {
                deg2D[i, 0] = (double)(deghere[i]);
                deg2D[i, 1] = 90;
            }
            double[,] deg;
            deg = deg2D;
            if (1 != dim[2])
            {
                double[,] deg3D = new double[numdeghere + 1, 2];
                for (int i = 1; i < numdeghere; i++)
                {
                    deg3D[i, 0] = deg2D[i, 0];
                    deg3D[i, 1] = deg2D[i, 1];
                }
                deg3D[numdeghere, 0] = 0;
                deg3D[numdeghere, 1] = 0;
                deg = deg3D;
            }
            //the half filter size
            int p = iPrmD / 2;
            //make filter
            MyFilter[] c = makefilter(deg, p);
            //number of directions
            int numdir = c.Length;
            //inisization
            Mat filtim = new Mat(im.Rows, im.Cols, MatType.CV_32F, 0.0);
            //zzzzz
            //201118 for Acceleration   Mat imRs = new Mat();
            Mat imRs = null;
            float[] afImRs = null;
            int iProportion = 5;
            if (1 == iInterpolationMethod)
            {
                imRs = new Mat();//201118 for Acceleration
                Mat imFt = new Mat();
                im.ConvertTo(imFt, MatType.CV_32F);
                Cv2.Resize(imFt, imRs, new Size(imFt.Cols * iProportion, imFt.Rows * iProportion), 0, 0, InterpolationFlags.Cubic);
                //afImRs = new float[imRs.Rows, imRs.Cols];
                afImRs = new float[imRs.Rows * imRs.Cols];
                Marshal.Copy(imRs.Data, afImRs, 0, imRs.Rows * imRs.Cols);
                /*
            for (int i = 0; i < imRs.Rows; i++)
            {
                for (int j = 0; j < imRs.Cols; j++)
                {
                    afImRs[i, j] = imRs.At<float>(i, j);
                }

                 */
                imRs = null;
            }
            for (int i = 0; i < numdir; i++)
            {
                //coordinates of points in this direction after rotation of filter
                MyFilter chere = c[i];
                //filter image for max value
                //zzzzz
                OUTmaxfilt3 fltrst;
                if (1 == iInterpolationMethod)
                {
                    //201118 for Acceleration   fltrst = maxfilt3ByCubicSpline(im, chere, numremele/*, dim*/, imRs, iProportion);
                    fltrst = maxfilt3ByCubicSpline(im, chere, numremele/*, dim*/, afImRs, iProportion);
                }
                else if (2 == iInterpolationMethod)
                {
                    fltrst = maxfilt3ByBilinear(im, chere, numremele/*, dim*/);
                }
                else
                {
                    fltrst = maxfilt3ByNearest(im, chere, numremele/*, dim*/);
                }

                //structural filtering   
                //filtim = max(filtim,(sumim - sum(maxim,4))/ (numpoints-numremele));
                int[] point = new int[2];
                for (point[0] = 0; point[0] < filtim.Rows; point[0]++)
                {
                    for (point[1] = 0; point[1] < filtim.Cols; point[1]++)
                    {
                        float fTemp = filtim.At<float>(point[0], point[1]);
                        /*201118 for Acceleration 
                        fTemp = Math.Max(fTemp,
                            (fltrst.sumim[point[0], point[1]] - (fltrst.maxim[point[0], point[1]]).Sum()) / (fltrst.numpoints - numremele)
                                        );
                        */
                        //(fltrst.minim[point[0], point[1]]).Sum()
                        fTemp = fTemp = Math.Max(fTemp, (fltrst.minim[point[0]][point[1]].Take<float>(fltrst.numpoints - numremele)).Sum() / (fltrst.numpoints - numremele));
                        filtim.Set<float>(point, fTemp);
                    }
                }
            }
            return filtim;
        }
        OUTmaxfilt3 maxfilt3ByNearest(Mat im, MyFilter chere, int numremele/*, int[] dim*/)
        {
            /////Mat imRs = new Mat();
            /////im.ConvertTo(imRs, MatType.CV_32F);
            int BigRows = im.Rows;
            int BigCols = im.Cols;
            OUTmaxfilt3 rst = new OUTmaxfilt3();
            rst.sumim = new float[im.Rows, im.Cols];
            rst.maxim = new List<float>[im.Rows, im.Cols];
            for (int i = 0; i < im.Rows; i++)
            {
                for (int j = 0; j < im.Cols; j++)
                {
                    (rst.sumim)[i, j] = 0;
                    List<float> lisTmp = new List<float>();
                    for (int k = 0; k < numremele; k++)
                    {
                        lisTmp.Add(0);
                    }
                    (rst.maxim)[i, j] = lisTmp;
                }
            }
            rst.numpoints = chere.x.Length;
            for (int i = 0; i < rst.numpoints; i++)
            {
                double xhere;
                double yhere;
                int[] point = new int[2];
                for (point[0] = 0; point[0] < im.Rows; point[0]++)
                {
                    for (point[1] = 0; point[1] < im.Cols; point[1]++)
                    {
                        xhere = (point[0] + (chere.x)[i]);
                        yhere = (point[1] + (chere.y)[i]);
                        int iX = Math.Min(Math.Max((int)(xhere + 0.5), 0), BigRows - 1);
                        int iY = Math.Min(Math.Max((int)(yhere + 0.5), 0), BigCols - 1);
                        float fTemp = im.At<byte>(iX, iY);
                        (rst.sumim)[point[0], point[1]] += fTemp;
                        (rst.maxim)[point[0], point[1]].Add(fTemp);
                        (rst.maxim)[point[0], point[1]].Sort();
                        (rst.maxim)[point[0], point[1]].RemoveAt(0);
                        (rst.maxim)[point[0], point[1]].Reverse();
                    }
                }
            }
            return rst;
        }
        OUTmaxfilt3 maxfilt3ByBilinear(Mat im, MyFilter chere, int numremele/*, int[] dim*/)
        {
            /////Mat imRs = new Mat();
            /////im.ConvertTo(imRs, MatType.CV_32F);
            int BigRows = im.Rows;
            int BigCols = im.Cols;

            OUTmaxfilt3 rst = new OUTmaxfilt3();
            rst.sumim = new float[im.Rows, im.Cols];
            rst.maxim = new List<float>[im.Rows, im.Cols];
            for (int i = 0; i < im.Rows; i++)
            {
                for (int j = 0; j < im.Cols; j++)
                {
                    (rst.sumim)[i, j] = 0;
                    List<float> lisTmp = new List<float>();
                    for (int k = 0; k < numremele; k++)
                    {
                        lisTmp.Add(0);
                    }
                    (rst.maxim)[i, j] = lisTmp;
                }
            }
            rst.numpoints = chere.x.Length;
            for (int i = 0; i < rst.numpoints; i++)
            {
                double xhere;
                double yhere;
                int[] point = new int[2];
                for (point[0] = 0; point[0] < im.Rows; point[0]++)
                {
                    for (point[1] = 0; point[1] < im.Cols; point[1]++)
                    {
                        xhere = point[0] + (chere.x)[i];
                        yhere = point[1] + (chere.y)[i];
                        int iXS = Math.Min(Math.Max((int)(xhere), 0), BigRows - 1);
                        int iYS = Math.Min(Math.Max((int)(yhere), 0), BigCols - 1);
                        int iXB = Math.Min(Math.Max((int)(xhere + 1.0), 0), BigRows - 1);
                        int iYB = Math.Min(Math.Max((int)(yhere + 1.0), 0), BigCols - 1);
                        double dWXB = Math.Min(Math.Max((xhere - (double)iXS), 0.0), 1.0);
                        double dWXS = 1.0 - dWXB;
                        double dWYB = Math.Min(Math.Max((yhere - (double)iYS), 0.0), 1.0);
                        double dWYS = 1.0 - dWYB;
                        double dTempA = im.At<byte>(iXS, iYS);
                        double dTempB = im.At<byte>(iXS, iYB);
                        double dTempC = im.At<byte>(iXB, iYS);
                        double dTempD = im.At<byte>(iXB, iYB);
                        double dTempAB = dTempA * dWYS + dTempB * dWYB;
                        double dTempCD = dTempC * dWYS + dTempD * dWYB;
                        float fTemp = (float)(dTempAB * dWXS + dTempCD * dWXB);
                        (rst.sumim)[point[0], point[1]] += fTemp;
                        (rst.maxim)[point[0], point[1]].Add(fTemp);
                        (rst.maxim)[point[0], point[1]].Sort();
                        (rst.maxim)[point[0], point[1]].RemoveAt(0);
                        (rst.maxim)[point[0], point[1]].Reverse();
                    }
                }
            }
            return rst;
        }
        //201118 for Acceleration   OUTmaxfilt3 maxfilt3ByCubicSpline(Mat im, MyFilter chere, int numremele/*, int[] dim*/, Mat imRs, int iProportion)
        OUTmaxfilt3 maxfilt3ByCubicSpline(Mat im, MyFilter chere, int numremele, float[] afImRs, int iProportion)
        {
            /////Mat imFt = new Mat();
            /////Mat imRs = new Mat();
            /////im.ConvertTo(imFt, MatType.CV_32F);
            /////int iProportion = 5;
            /////Cv2.Resize(imFt, imRs, new Size(imFt.Cols * iProportion, imFt.Rows * iProportion), 0, 0, InterpolationFlags.Cubic);
            //201118 for Acceleration   int BigRows = imRs.Rows;
            //201118 for Acceleration   int BigCols = imRs.Cols;
            int BigRows = im.Rows * iProportion;
            int BigCols = im.Cols * iProportion;
            /*
            double[,] x = new double[dim[0], dim[1]];
            double[,] y = new double[dim[0], dim[1]];
            for (int i = 0; i < dim[0]; i++)
            {
                for (int j = 0; j < dim[1]; j++)
                {
                    x[i, j] = (double)i;
                    y[i, j] = (double)j;
                }
            }
            */
            OUTmaxfilt3 rst = new OUTmaxfilt3();
            //rst.maxim = new Mat[numremele];
            //201118 for Acceleration   rst.sumim = new float[im.Rows, im.Cols];
            //201118 for Acceleration   rst.maxim = new List<float>[im.Rows, im.Cols];
            rst.numpoints = chere.x.Length;
            rst.minim = new float[im.Rows][][];
            for (int i = 0; i < im.Rows; i++)
            {
                rst.minim[i] = new float[im.Cols][];
                for (int j = 0; j < im.Cols; j++)
                {
                    rst.minim[i][j] = new float[rst.numpoints];
                    //201118 for Acceleration   (rst.sumim)[i, j] = 0;
                    //201118 for Acceleration   List<float> lisTmp = new List<float>();

                    //201118 for Acceleration   for (int k = 0; k < numremele; k++)
                    //201118 for Acceleration   {
                    //201118 for Acceleration       lisTmp.Add(0);
                    //201118 for Acceleration   }

                    //201118 for Acceleration   (rst.maxim)[i, j] = lisTmp;
                    //201118 for Acceleration   List<float> lisTmpMin = new List<float>();
                    //201118 for Acceleration   (rst.minim)[i, j] = lisTmpMin;
                }
            }

            //Mat imhere = new Mat(imFt.Rows, imFt.Cols, MatType.CV_32F, 0.0);
            /* 201118 for Acceleration
            for (int i = 0; i < rst.numpoints; i++)
            {
                //the coordinate where we want to find the values
                //NB!!!!!!!!!
                //Must not do this for z direction, then it goes wrong.
                //Do NOOOOOT replacestr any Inf or NaN by image values, then the filter
                //does not work anymore!!!!!!
                //xhere(xhere < h(1)) = h(1);xhere(xhere > dim(1)*h(1)) = dim(1)*h(1);
                //yhere(yhere < h(1)) = h(1);yhere(yhere > dim(2)*h(1)) = dim(2)*h(2);

                //double[,] xhere = new double[dim[0], dim[1]];
                //double[,] yhere = new double[dim[0], dim[1]];
                double xhere;
                double yhere;
                int[] point = new int[2];
                for (point[0] = 0; point[0] < im.Rows; point[0]++)
                {
                    for (point[1] = 0; point[1] < im.Cols; point[1]++)
                    {
                        //xhere[i, j] = Math.Min(Math.Max(x[i, j]+(chere.x)[t],1),dim[0]);
                        //yhere[i, j] = Math.Min(Math.Max(y[i, j]+(chere.y)[t],1),dim[1]);
                        xhere = (point[0] + (chere.x)[i]) * iProportion;
                        yhere = (point[1] + (chere.y)[i]) * iProportion;
                        int iX = Math.Min(Math.Max((int)(xhere + 0.5), 0), BigRows - 1);
                        int iY = Math.Min(Math.Max((int)(yhere + 0.5), 0), BigCols - 1);
                        float fTemp = imRs.At<float>(iX, iY);
                        //imhere.Set<float>(point, fTemp);
                        (rst.sumim)[point[0], point[1]] += fTemp;
                        (rst.maxim)[point[0], point[1]].Add(fTemp);
                        (rst.maxim)[point[0], point[1]].Sort();
                        (rst.maxim)[point[0], point[1]].RemoveAt(0);
                        (rst.maxim)[point[0], point[1]].Reverse();
                    }
                }
            }
            */
            for (int ir = 0; ir < im.Rows; ir++)
            {
                for (int ic = 0; ic < im.Cols; ic++)
                {
                    double xhere;
                    double yhere;
                    for (int ip = 0; ip < rst.numpoints; ip++)
                    {
                        xhere = (ir + (chere.x)[ip]) * iProportion;
                        yhere = (ic + (chere.y)[ip]) * iProportion;
                        int iX = Math.Min(Math.Max((int)(xhere + 0.5), 0), BigRows - 1);
                        int iY = Math.Min(Math.Max((int)(yhere + 0.5), 0), BigCols - 1);
                        //201118 for Acceleration   float fTemp = imRs.At<float>(iX, iY);
                        float fTemp = afImRs[iX * BigCols + iY];
                        //201118 for Acceleration   (rst.sumim)[ir, ic] += fTemp;
                        //201118 for Acceleration   (rst.maxim)[ir, ic].Add(fTemp);
                        (rst.minim)[ir][ic][ip] = fTemp;
                    }
                    //201118 for Acceleration   (rst.maxim)[ir, ic].Sort((x, y) => -x.CompareTo(y));
                    //201118 for Acceleration   (rst.maxim)[ir, ic].RemoveRange(numremele, rst.numpoints - numremele);
                    Array.Sort((rst.minim)[ir][ic]);
                    //201118 for Acceleration   (rst.minim)[ir, ic].Sort();
                    //201118 for Acceleration   (rst.minim)[ir, ic].RemoveRange(rst.numpoints - numremele,numremele);

                }
            }
            return rst;
        }
        MyFilter[] makefilter(double[,] deg, int p)
        {
            int iNum = deg.GetLength(0);
            MyFilter[] c = new MyFilter[iNum];
            //r = linspace(-p,p,7);
            int iArrLength = 7;
            double[] r = new double[iArrLength];
            double dStep = p * 2.0 / (iArrLength - 1.0);
            for (int i = 0; i < 7; i++)
            {
                r[i] = -1.0 * p + i * dStep;
            }
            /*
             for i = 1 : size(deg,1)
                % the degrees
                theta = deg(i,1);
                phi = deg(i,2);        
                % the positions
                c(i).x = r*sind(phi)*cosd(theta);
                c(i).y = r*sind(phi)*sind(theta);
                c(i).z = r*cosd(phi);
            end; 
            */
            for (int i = 0; i < iNum; i++)
            {
                double theta = deg[i, 0];
                double phi = deg[i, 1];
                c[i] = new MyFilter();
                c[i].x = new double[iArrLength];
                c[i].y = new double[iArrLength];
                c[i].z = new double[iArrLength];
                for (int j = 0; j < iArrLength; j++)
                {
                    (c[i].x)[j] = r[j] * Math.Sin(phi) * Math.Cos(theta); ;
                    (c[i].y)[j] = r[j] * Math.Sin(phi) * Math.Sin(theta); ;
                    (c[i].z)[j] = r[j] * Math.Cos(phi);
                }
            }
            return c;
        }
        public Mat ProFroStandardImage(Mat GryImg, int iConnect, Color phaseColor)
        {
            //Mat GryImg = new Mat();
            ////转为为灰度?
            //if (3 == OrgImg.Channels())
            //{
            //    Cv2.CvtColor(OrgImg, GryImg, ColorConversionCodes.RGB2GRAY);
            //}
            //else
            //{
            //    if (1 == OrgImg.Channels())
            //    {
            //        GryImg = OrgImg.Clone();
            //    }
            //    else
            //    {
            //        return null;
            //    }
            //}
            Mat CrystalImg = new Mat();
            Mat BoundaryImg = new Mat();
            Mat mat_mean = new Mat(), mat_stddev = new Mat();
            double m, s, d = 0.5;
            Cv2.MeanStdDev(GryImg, mat_mean, mat_stddev);
            m = mat_mean.At<double>(0, 0);
            s = mat_stddev.At<double>(0, 0);
            //Cv2.Threshold(GryImg, CrystalImg, (int)Math.Min(m+d*s,230), 255, ThresholdTypes.Binary);
            //Cv2.Threshold(GryImg, CrystalImg, 0, 1, ThresholdTypes.Otsu);
            Cv2.Threshold(GryImg, CrystalImg, m, 255, ThresholdTypes.Binary);
            //Cv2.ImShow("src", CrystalImg);//kkkkk
            //Cv2.WaitKey(0);//kkkkk
            int diam = 3;//mmmmm
            Mat kernel = Cv2.GetStructuringElement(MorphShapes.Ellipse, new Size(diam, diam), new OpenCvSharp.Point(-1, -1));
            //Cv2.Erode(CrystalImg, CrystalImg, kernel);
            Cv2.MorphologyEx(CrystalImg, CrystalImg, MorphTypes.Open, kernel);
            Cv2.BitwiseNot(CrystalImg, BoundaryImg);
            //kernel = Cv2.GetStructuringElement(MorphShapes.Ellipse, new Size(diam, diam), new Point(-1, -1));
            //Cv2.MorphologyEx(BoundaryImg, BoundaryImg, MorphTypes.Open, kernel);
            //Cv2.MorphologyEx(BoundaryImg, BoundaryImg, MorphTypes.Close, kernel);
            //Cv2.Dilate(BoundaryImg, BoundaryImg, kernel);
            Class1 TCAlgo = new Class1();
            TCAlgo.ImgThin(BoundaryImg);
            //string filename = "D:\\temp.jpg";
            //Cv2.ImWrite(filename, BoundaryImg);
            if (1 == iConnect)
            {
                TCAlgo.ConnectionBoundary(BoundaryImg, 108);
            }
            //Cv2.ImShow("wai", BoundaryImg);//kkkkk
            //Cv2.WaitKey(0);//kkkkk
            //Cv2.BitwiseNot(BoundaryImg, CrystalImg);
            //diam = 9;
            //kernel = Cv2.GetStructuringElement(MorphShapes.Cross, new Size(diam, diam), new Point(-1, -1));
            //Cv2.MorphologyEx(BoundaryImg, BoundaryImg, MorphTypes.Close, kernel);
            //Cv2.Erode(CrystalImg, CrystalImg, kernel);
            //Cv2.MorphologyEx(CrystalImg, CrystalImg, MorphTypes.Open, kernel);
            //Cv2.MorphologyEx(CrystalImg, CrystalImg, MorphTypes.Close, kernel);
            //Cv2.ImShow("nei", CrystalImg);//kkkkk
            //Cv2.WaitKey(0);//kkkkk
            /*
            Mat watermark = new Mat(CrystalImg.Size(), MatType.CV_32S);
            Point[][] contour;
            HierarchyIndex[] hier;
            Cv2.FindContours(CrystalImg, out contour, out hier, RetrievalModes.CComp, ContourApproximationModes.ApproxSimple, null);
            for (int i = 0; i < hier.Length; i++)
            {
                Cv2.DrawContours(watermark, contour, i, Scalar.All(i + 1), 1, LineTypes.Link8, hier);
            }
            Cv2.Watershed(OrgImg, watermark);
            BoundaryImg = new Mat(watermark.Size(), MatType.CV_8UC3);
            //ttttt Cv2.DrawContours(rtn.wat, contour, -1, Scalar.Red);
            */
            /*
            Mat NewWat = new Mat();
            Cv2.BitwiseNot(CrystalImg, NewWat);
            Cv2.ImShow("bian", NewWat);//kkkkk
            Cv2.WaitKey(0);//kkkkk
            ImgThin(NewWat);
            int Rows = NewWat.Rows;
            int Cols = NewWat.Cols;
            */
            int[] point = new int[2];
            Vec4b redcol;
            if (phaseColor != null)
            {
                redcol = new Vec4b(phaseColor.B, phaseColor.G, phaseColor.R, 255);
            }
            else
                redcol = new Vec4b(0, 0, 255/*255*/, 255);
            Mat WAT = new Mat(BoundaryImg.Size(), MatType.CV_8UC4);
            for (point[0] = 0; point[0] < WAT.Rows; point[0]++)
            {
                for (point[1] = 0; point[1] < WAT.Cols; point[1]++)
                {
                    if (0 < BoundaryImg.At<byte>(point[0], point[1]))
                    {
                        WAT.Set<Vec4b>(point, redcol);
                    }
                }
            }
            return WAT;
        }
    }
}