WdConvolutionVm.cs

using FlyADBase;
using GalaSoft.MvvmLight.Command;
using MathNet.Numerics.LinearAlgebra.Double;
using System; 
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.ComponentModel;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Flyad7_WPF
{
    public class WdConvolutionVm : INotifyPropertyChanged
    {
        public event PropertyChangedEventHandler PropertyChanged;


        public RelayCommand CreateHolderCmd { get; private set; }
        public RelayCommand CalConvCmd { get; private set; }
        public RelayCommand CalConv2Cmd { get; private set; }
        public RelayCommand UpdateConvCmd { get; private set; }
        /// <summary>
        /// 小托辊不锈钢支架 宽度单位 mm,  默认15mm 
        /// </summary>
        public int HolderWidth { get; set; } = 15;

        /// <summary>
        /// 编码器1 的比例 mm/脉冲
        /// </summary>
        public double Mmpp { get; set; } = 0.16;

        
        /// <summary>
        /// 小托辊不锈钢支架 中间grid位置
        /// </summary>
        public int HolderMidGrid { get; set; } = 40;
        /// <summary>
        /// 小托辊不锈钢支架 开始grid位置
        /// </summary>
        public int HolderBeginGrid { get; set; } = 20;
        /// <summary>
        /// 小托辊不锈钢支架 结束grid位置
        /// </summary>
        public int HolderEndGrid { get; set; } = 60;
        public int HolderLenGrid => HolderEndGrid - HolderBeginGrid + 1;
        /// <summary>
        /// 空气AD
        /// </summary>
        public int AirAd { get; set; } = 25900;
        /// <summary>
        /// 小托辊不锈钢支架  grid图
        /// </summary>
        public int[] HolderDatas { get; set; }

        /// <summary>
        /// 小托辊不锈钢支架  滤波后的 grid图， 就是AD盒采集的
        /// </summary>
        public int[] FilterDatas { get; set; }

        /// <summary>
        /// 把滤波后的数据，反卷积后的 grid图
        /// </summary>
        public int[] DeconvDatas { get; set; }


        /// <summary>
        /// 卷积核 大小，单位grid
        /// </summary>
        public ObservableCollection<CoreValue> ConvCore { get; } = new ObservableCollection<CoreValue>();

        /// <summary>
        /// 卷积核 大小，单位grid
        /// </summary>
        public int CoreLen { get; set; } = 17;


        /// <summary>
        /// DeconvDatas 与 HolderDatas 的相关性
        /// </summary>
        public double R { get; set; }
        /// <summary>
        /// FilterDatas 与 HolderDatas 的相关性
        /// </summary>
        public double HolderR { get; set; }

        System.Windows.Forms.DataVisualization.Charting.Chart chart0;
        System.Windows.Forms.DataVisualization.Charting.Chart chart1;
        System.Windows.Forms.DataVisualization.Charting.Chart chart2;
        IFlyADClientAdv flyad;
        
        public WdConvolutionVm()
        {
            CreateHolderCmd = new RelayCommand(CreateHolder);
            CalConvCmd = new RelayCommand(CalConv);
            CalConv2Cmd = new RelayCommand(CalConv2);
            UpdateConvCmd = new RelayCommand(UpdateConv);
        }



        public void Init(
            IFlyADClientAdv flyad,
            System.Windows.Forms.DataVisualization.Charting.Chart chart0,
            System.Windows.Forms.DataVisualization.Charting.Chart chart1,
            System.Windows.Forms.DataVisualization.Charting.Chart chart2,
            int[] filterDatas)
        {
            this.flyad = flyad;
            this.chart0 = chart0;
            this.chart1 = chart1;
            this.chart2 = chart2;
            FilterDatas = filterDatas;
            DrawFilter();
        }
        private async void DrawFilter()
        {
            await Draw0(0, FilterDatas);
        }

        async Task Draw0(int index, int[] data)
        {
            var chart = chart0;
            await Task.Yield();
            //画图
            await Task.Factory.StartNew(() =>
            {
                App.Current.Dispatcher.Invoke(() =>
                {
                    System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[index];
                    series.Points.Clear();
                    for (int i = 0; i < data.Count(); i++)
                    {
                        if (Misc.MyBase.ISVALIDATA(data[i]))
                            series.Points.AddXY(i, data[i]);
                    }
                    //if (series.Points.Count() > 0)
                    //{
                    //    double max = series.Points.Max(p => p.YValues[0]);
                    //    double min = series.Points.Min(p => p.YValues[0]);
                    //    if (max != min)
                    //    {
                    //        double range = (max - min) / 10;
                    //        chart.ChartAreas[0].AxisY.Minimum = min - range;
                    //        chart.ChartAreas[0].AxisY.Maximum = max + range;
                    //    }
                    //}
                });
            });
        }

        async Task Draw2(int index, double r)
        {
            var chart = chart2;
            //画图
            await Task.Yield();
            await Task.Factory.StartNew(() =>
            {
                App.Current.Dispatcher.Invoke(() =>
                {
                    System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[index];
                    series.Points.AddY(r * 10000);
                    double max = series.Points.Max(p => p.YValues[0]);
                    double min = series.Points.Min(p => p.YValues[0]);
                    if (max != min)
                    {
                        double range = (max - min) / 10;
                        chart.ChartAreas[0].AxisY.Minimum = min - range;
                        chart.ChartAreas[0].AxisY.Maximum = max + range;
                    }
                });
            });
        }

        async Task Draw1(int index,IEnumerable<double> data)
        {
            var chart = chart1;
            //画图
            await Task.Yield();
            await Task.Factory.StartNew(() =>
            {
                App.Current.Dispatcher.Invoke(() =>
                {
                    System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[index];

                    series.Points.Clear();
                    for (int i = 0; i < data.Count(); i++)
                    {
                        series.Points.AddXY(i, data.ElementAt(i));
                    }

                    double max = series.Points.Max(p => p.YValues[0]);
                    double min = series.Points.Min(p => p.YValues[0]);
                    if (max != min)
                    {
                        double range = (max - min) / 10;
                        chart.ChartAreas[0].AxisY.Minimum = min - range;
                        chart.ChartAreas[0].AxisY.Maximum = max + range;
                    }
                });
            });
        }
        private async void CalConv()
        {
            //解卷积核
            // Holder * M = Filter;
            // A * X = B;
            DenseMatrix matrixA = new DenseMatrix(HolderLenGrid, CoreLen);
            DenseVector vectorB = new DenseVector(HolderLenGrid);

            for (int i = 0; i < HolderLenGrid; i++)
            {
                DenseVector vector = new DenseVector(CoreLen);
                for (int j = 0; j < CoreLen; j++) 
                {
                    vector[j] = HolderDatas[HolderBeginGrid + i + j - CoreLen / 2];
                }
                matrixA.SetRow(i, vector);
                vectorB[i] = FilterDatas[HolderBeginGrid + i];
            }


            //A~AX = A~B
            DenseMatrix ATA = (DenseMatrix)matrixA.TransposeThisAndMultiply(matrixA);
            DenseVector ATB = (DenseVector)matrixA.TransposeThisAndMultiply(vectorB);

            //DenseVector X = (DenseVector)matrixA.LU().Solve(vectorB);
            DenseVector X = (DenseVector)ATA.LU().Solve(ATB);
            ConvCore.Clear();

            //导出到 ConvCore
            for (int i = 0; i < CoreLen; i++)
            {
                ConvCore.Add(new CoreValue() { Value = X[i] });
            }
            await Draw1(0, ConvCore.Select(c=>c.Value));

            var filterdatas2 = Conv(HolderDatas, ConvCore.Select(c => c.Value));
            await Draw0(2, filterdatas2);
            R = Misc.MyMath.Correl(FilterDatas, filterdatas2);
            Draw2Clear();
            await Draw2(0, R);
        }

        async void Draw2Clear()
        {
            var chart = chart2;
            await Task.Yield();
            //画图
            App.Current.Dispatcher.Invoke(() =>
            {
                System.Windows.Forms.DataVisualization.Charting.Series series = chart.Series[0];
                series.Points.Clear();
                //chart.ChartAreas[0].AxisX.Minimum = double.NaN;
                //chart.ChartAreas[0].AxisX.Maximum = double.NaN;
            });
        }
        private async void CalConv2()
        {
            double step = 0.00001;
            int decimals = 5;
            //在 ConvCore 基础上微调
            
            //必须为单数
            if (ConvCore.Count % 2 != 1)
                ConvCore.RemoveAt(ConvCore.Count - 1);

            //必须都大于等于0


            //ConvCore 必须是对称的
            for (int i = 0; i < ConvCore.Count() / 2; i++)
            {
                if (ConvCore[i].Value < 0)
                    ConvCore[i].Value = 0;
                ConvCore[i].Value = Math.Round(ConvCore[i].Value, decimals);

                ConvCore[ConvCore.Count - 1 - i].Value = ConvCore[i].Value;
            }



            while (true)
            {
                var filterdatas2 = Conv(HolderDatas, ConvCore.Select(c => c.Value));
                double r0 = Misc.MyMath.Correl(FilterDatas, filterdatas2);

                bool hasChanged = false;

                double[] ks = ConvCore.Select(c => c.Value).ToArray();
                double r = r0;
                await Task.Factory.StartNew(() =>
                {
                    //每个参数改变 step, 确定方向
                    for (int i = 0; i <(ks.Count()/2+1); i++)
                    {
                        double curr_step = step;
                        double r1;
                        double r2;
                        while (true)
                        {
                            ks[i] = ConvCore[i].Value + curr_step;
                            ks[ks.Count() - 1 - i] = ks[i];
                            filterdatas2 = Conv(HolderDatas, ks);
                            r1 = Misc.MyMath.Correl(FilterDatas, filterdatas2);


                            ks[i] = ConvCore[i].Value - curr_step;
                            ks[ks.Count() - 1 - i] = ks[i];
                            filterdatas2 = Conv(HolderDatas, ks);
                            r2 = Misc.MyMath.Correl(FilterDatas, filterdatas2);

                            if ((r1 == r0) || (r2 == r0))
                            {
                                //居然没有变化，把step翻倍
                                curr_step *= 1.5;
                            }
                            else
                            {
                                //有变化了
                                break;
                            }
                        }

                        if (r1 > r0)
                        {
                            //方向正确
                            ks[i] = ConvCore[i].Value + curr_step;
                            ks[ks.Count() - 1 - i] = ks[i];
                            hasChanged = true;
                            r = r1;
                            continue;
                        }
                        else if (r2 > r0)
                        {
                            //方向正确
                            ks[i] = ConvCore[i].Value - curr_step;
                            ks[ks.Count() - 1 - i] = ks[i];
                            hasChanged = true;
                            r = r2;
                            continue;
                        }
                        else
                        {
                            //方向错误，回到以前
                            ks[i] = ConvCore[i].Value;
                        }
                    }
                });

                R = r;
                await Draw2(0, r);

                if (!hasChanged)
                {
                    //已经计算完毕
                    break;
                }
                else 
                {
                    for (int i = 0; i < ks.Count(); i++)
                        ConvCore[i].Value = ks[i];
                    await Draw1(0, ConvCore.Select(c => c.Value));
                }
            }

            UpdateConv();
        }
        private async void UpdateConv()
        {

            await Draw1(0, ConvCore.Select(c => c.Value));
            var filterdatas2 = Conv(HolderDatas, ConvCore.Select(c => c.Value));
            await Draw0(2, filterdatas2);
            R = Misc.MyMath.Correl(FilterDatas, filterdatas2);
            await Draw2(0, R);
        }
        
        int[] Conv(int[] a, IEnumerable<double> core) 
        {
            int[] b = new int[a.Length];
            for (int i = 0; i < a.Length; i++) 
            {
                double sum=0;
                double cnt=0;
                for (int j = 0; j < core.Count(); j++) 
                {
                    int index = i -  core.Count() / 2 + j;
                    if (index < 0)
                        continue;
                    if (index > a.Length - 1)
                        break;
                    if (Misc.MyBase.ISVALIDATA(a[index])){
                        cnt += core.ElementAt(j);
                        sum += core.ElementAt(j) * a[index]; 
                    }
                }
                if (cnt > 0)
                {
                    b[i] = (int)(sum / cnt);
                }
                else 
                {
                    b[i] = Misc.MyBase.NULL_VALUE;
                }
            }
            return b;
        }
        private async void CreateHolder()
        {
            HolderDatas = new int[FilterDatas.Count()];

            //把 Target 的 测到不锈钢部分删除， 剩下求平均值

            int sensorGridWidth = (int)(25 / Mmpp / flyad.PosOfGrid);
            int holderGridWidth = (int)(HolderWidth / Mmpp / flyad.PosOfGrid);

            for (int i = 0; i < HolderDatas.Count(); i++)
            {
                if (Misc.MyBase.ISVALIDATA(FilterDatas[i]))// && i>= HolderBeginGrid && i<= HolderEndGrid)
                    HolderDatas[i] = AirAd;
                else
                    HolderDatas[i] = Misc.MyBase.NULL_VALUE;
            }

            for (int i = 0; i < holderGridWidth; i++)
            {
                int index = HolderMidGrid + i - holderGridWidth / 2;
                HolderDatas[index] = 0;
            }

            //找到最大相关性对应的位置
            HolderR = Misc.MyMath.Correl(FilterDatas, HolderDatas);

            await Draw0(1, HolderDatas);
        }
    }

    public class CoreValue :INotifyPropertyChanged
    {
        public double Value { get; set; }

        public event PropertyChangedEventHandler PropertyChanged;
    }
}