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; }
///
/// 小托辊不锈钢支架 宽度单位 mm, 默认15mm
///
public int HolderWidth { get; set; } = 15;
///
/// 编码器1 的比例 mm/脉冲
///
public double Mmpp { get; set; } = 0.16;
///
/// 小托辊不锈钢支架 中间grid位置
///
public int HolderMidGrid { get; set; } = 40;
///
/// 小托辊不锈钢支架 开始grid位置
///
public int HolderBeginGrid { get; set; } = 20;
///
/// 小托辊不锈钢支架 结束grid位置
///
public int HolderEndGrid { get; set; } = 60;
public int HolderLenGrid => HolderEndGrid - HolderBeginGrid + 1;
///
/// 空气AD
///
public int AirAd { get; set; } = 25900;
///
/// 小托辊不锈钢支架 grid图
///
public int[] HolderDatas { get; set; }
///
/// 小托辊不锈钢支架 滤波后的 grid图, 就是AD盒采集的
///
public int[] FilterDatas { get; set; }
///
/// 把滤波后的数据,反卷积后的 grid图
///
public int[] DeconvDatas { get; set; }
///
/// 卷积核 大小,单位grid
///
public ObservableCollection ConvCore { get; } = new ObservableCollection();
///
/// 卷积核 大小,单位grid
///
public int CoreLen { get; set; } = 17;
///
/// DeconvDatas 与 HolderDatas 的相关性
///
public double R { get; set; }
///
/// FilterDatas 与 HolderDatas 的相关性
///
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 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 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;
}
}