LIBMAS/DedicatedAnalysis/GrainSize/GBT6394_2017/OLDcellsegm_0909版本.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using OpenCvSharp;

namespace ConsoleApplication1
{
    /////
    /////
    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;

    }

    class OUTgetminima
    {
        public Mat minimacell;
        public Mat minima;
    }
    class OUTclassifycells
    {
        public Mat cellbw;
        public PClassifycells classifycells;
    }
    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(Mat GryImg, double minvol, double maxvol, 
            string[] INnamehere, In5segmsurf IN5varhere,In6segmsurf IN6varhere)
        {
            OUTsegmsurf rtn = new OUTsegmsurf();
            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;

            /*Cell segmentation start*/
            //Correcting for uneven illumination
            //Smoothing
            PSmoothim prminsmt = prm.smoothim.Clone();
            //kkkkk imsegm = cellsegm.smoothim(imsegm,prm.smoothim.method,'prm',prm);
            //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);
            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_8UC3);
            //ttttt Cv2.DrawContours(rtn.wat, contour, -1, Scalar.Red);
            int Rows = watermark.Rows;
            int Cols = watermark.Cols;
            int[] point = new int[2];
            Vec3b redcol = new Vec3b(0,0,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<Vec3b>(point, redcol);
                    }
                }
            }
            //Cv2.ImShow("contour", rtn.wat);//kkkkk
            //Cv2.WaitKey(0);//kkkkk

            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();
            Cv2.MeanStdDev(imsegm, mat_mean, mat_stddev);
            double m, s;
            m = mat_mean.At<double>(0, 0);
            s = mat_stddev.At<double>(0, 0);
            double d = 1.0;//mmmmm
            double th = m + d * s;
            Mat minima = new Mat(), minimacell = new Mat();
            Cv2.Threshold(imsegm, minimacell, th, 255, ThresholdTypes.Binary);
            int diam;//mmmmm
            Mat kernel;
            diam = 11;
            kernel = Cv2.GetStructuringElement(MorphShapes.Cross, new Size(diam, diam), new Point(-1, -1));
            Cv2.Dilate(minimacell, minimacell, kernel);
            diam = 5;
            kernel = Cv2.GetStructuringElement(MorphShapes.Rect, new Size(diam, diam), new Point(-1, -1));
            Cv2.Erode(minimacell, minimacell, kernel);
            //Cv2.MorphologyEx(minimacell, minimacell, MorphTypes.Open, kernel);
            //Cv2.ImShow("bianjie", minimacell);
            //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 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 = 10;
            kernel = Cv2.GetStructuringElement(MorphShapes.Rect, new Size(diam, diam), new Point(-1, -1));
            Cv2.MorphologyEx(minima, minima, MorphTypes.Open, kernel);
            //Cv2.ImShow("neibu", minima);//kkkkk
            //Cv2.WaitKey(0);//kkkkk
            return minima;
        }
    }
}