LIBMAS/DedicatedAnalysis/GrainSize/GBT6394_2017/OLDcellsegm_0916版本.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 OUTmaxfilt3
    {
        public List<float>[,] maxim;
        public float[,] sumim;
        public int numpoints;
    }
    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(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();
            imsegm = smoothim(imsegm,prm.smoothim.method,"prm",prm);
            Cv2.ImShow("smoothim", imsegm);//kkkkk
            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);
            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();
            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);
            //Cv2.Threshold(imsegm, minimacell, 0, 255, ThresholdTypes.Otsu);
            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);
            //Cv2.Threshold(matTempB, minimacell, 0, 255, ThresholdTypes.Otsu);
            int diam;//mmmmm
            Mat kernel;
            diam = 9;
            kernel = Cv2.GetStructuringElement(MorphShapes.Cross, new Size(diam, diam), new Point(-1, -1));
            Cv2.Dilate(minimacell, minimacell, kernel);
            diam = 3;
            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 = 9;
            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;
        }
        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),rst);
            return rst;
        }
        Mat dircohenh(Mat im,int iPrmD,double[] adPrmH)
        {
            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);
            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
                OUTmaxfilt3 fltrst = maxfilt3(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]);
                        fTemp = Math.Max(fTemp,
                            (fltrst.sumim[point[0], point[1]] - (fltrst.maxim[point[0], point[1]]).Sum()) / (fltrst.numpoints - numremele)
                                        );
                        filtim.Set<float>(point, fTemp);
                    }
                }
            }      
            return filtim;
        }
        OUTmaxfilt3 maxfilt3(Mat im, MyFilter chere, int numremele/*, int[] dim*/)
        {
            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);
            int BigRows = imRs.Rows;
            int BigCols = imRs.Cols;
            /*
            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];
            rst.sumim = new float[imFt.Rows, imFt.Cols];
            rst.maxim = new List<float>[imFt.Rows, imFt.Cols];
            for (int i = 0; i < imFt.Rows; i++)
            {
                for (int j = 0; j < imFt.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;
            //Mat imhere = new Mat(imFt.Rows, imFt.Cols, MatType.CV_32F, 0.0);
            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] < imFt.Rows; point[0]++)
                {
                    for (point[1] = 0; point[1] < imFt.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();
                    }
                }
            }
            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;
        }
    }
}