using System;
using System.Collections.Generic;
using System.Text;
using System.IO;
using System.Drawing;
using System.Windows.Forms;
using System.Drawing.Drawing2D;
using System.Drawing.Imaging;
using MapProj;
namespace GeoFly
{
	/// <summary>
	/// 主要用于放置一些公用函数
	/// </summary>
	public static class MatrixFuncs
	{
		public static double[,] Float2Double(float[,] data)
		{
			int rows=data.GetLength(0);
			int cols=data.GetLength(1);
			double[,] result=new double[rows,cols];
			for(int row=0;row<rows;row++)
			{
				for(int col=0;col<cols;col++)
				{
					result[row,col]=(double)data[row,col];
				}
			}
			return result;
		
		}
		public static double[,] Int2Double(int[,] data)
		{
			int rows=data.GetLength(0);
			int cols=data.GetLength(1);
			double[,] result=new double[rows,cols];
			for(int row=0;row<rows;row++)
			{
				for(int col=0;col<cols;col++)
				{
					result[row,col]=(double)data[row,col];
				}
			}
			return result;
		
		}
		public static double[,] ReGridBilinear(double[,] data, double[] latsSource, double[] lonsSource, double[] newlats, double[] newlons)
		{
			int nlat = newlats.Length;
			int nlon = newlons.Length;
           
			double[,] result = new double[nlat, nlon];

			int[] lonIndice = new int[nlon];
			double dlat = latsSource[1] - latsSource[0];
			double dlon = lonsSource[1] - lonsSource[0];
			
			for (int i = 0; i < nlon; i++) {
				double lon = newlons[i];                
				lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);
			}
			for (int ilat = 0; ilat < nlat; ilat++) {
				double lat = newlats[ilat];
				int latIndex = (int)((lat - latsSource[0]) / dlat);
				for (int ilon = 0; ilon < nlon; ilon++) {
					//Calculate bilinear interpolation using nearest four grid cells
					double lon = newlons[ilon];
					int lonIndex = lonIndice[ilon];

					double A = data[latIndex, lonIndex];
					double B = A;
					if (lonIndex + 1 < lonsSource.Length)
						B = data[latIndex, lonIndex + 1];
					double C = A;
					double D = A;
					if (latIndex + 1 < latsSource.Length) {
						C = data[latIndex + 1, lonIndex];
						if (lonIndex + 1 < lonsSource.Length)
							D = data[latIndex + 1, lonIndex + 1];
					}
					
					double[] abcd = new double[]{ A, B, C, D };
					bool flag = false;
					for (int ind = 0; ind < 4; ind++) {
						if (abcd[ind] < -9990) {
							result[ilat, ilon] = Math.Max(A, Math.Max(B, Math.Max(C, D)));
							flag = true;
						}
					}
					if (flag)
						continue;
					double lonA = lonsSource[lonIndex];
					double lonB = lonA;
					if (lonIndex + 1 < lonsSource.Length)
						lonB = lonsSource[lonIndex + 1];
					double E = LinearInter(A, B, lonA, lonB, lon);
					double F = LinearInter(C, D, lonA, lonB, lon);

					double latA = latsSource[latIndex];
					double latB = latA;
					if (latIndex + 1 < latsSource.Length)
						latB = latsSource[latIndex + 1];
					double vvv = LinearInter(E, F, latA, latB, lat);
					result[ilat, ilon] = vvv;					
				}
			}
			return result;
		}

		public static float[,] ReGridBilinear(double[,] data, double[] latsSource, double[] lonsSource, float[] newlats, float[] newlons)
		{
			int nlat = newlats.Length;
			int nlon = newlons.Length;
           
			float[,] result = new float[nlat, nlon];

			int[] lonIndice = new int[nlon];
			double dlat = latsSource[1] - latsSource[0];
			double dlon = lonsSource[1] - lonsSource[0];
			
			for (int i = 0; i < nlon; i++) {
				double lon = newlons[i];                
				lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);
			}
			for (int ilat = 0; ilat < nlat; ilat++) {
				double lat = newlats[ilat];
				int latIndex = (int)((lat - latsSource[0]) / dlat);
				for (int ilon = 0; ilon < nlon; ilon++) {
					//Calculate bilinear interpolation using nearest four grid cells
					double lon = newlons[ilon];
					int lonIndex = lonIndice[ilon];

					double A = data[latIndex, lonIndex];
					double B = A;
					if (lonIndex + 1 < lonsSource.Length)
						B = data[latIndex, lonIndex + 1];
					double C = A;
					double D = A;
					if (latIndex + 1 < latsSource.Length) {
						C = data[latIndex + 1, lonIndex];
						if (lonIndex + 1 < lonsSource.Length)
							D = data[latIndex + 1, lonIndex + 1];
					}
					
					double[] abcd = new double[]{ A, B, C, D };
					bool flag = false;
					for (int ind = 0; ind < 4; ind++) {
						if (abcd[ind] < -9990) {
							result[ilat, ilon] = (float)Math.Max(A, Math.Max(B, Math.Max(C, D)));
							flag = true;
						}
					}
					if (flag)
						continue;
					double lonA = lonsSource[lonIndex];
					double lonB = lonA;
					if (lonIndex + 1 < lonsSource.Length)
						lonB = lonsSource[lonIndex + 1];
					double E = LinearInter(A, B, lonA, lonB, lon);
					double F = LinearInter(C, D, lonA, lonB, lon);

					double latA = latsSource[latIndex];
					double latB = latA;
					if (latIndex + 1 < latsSource.Length)
						latB = latsSource[latIndex + 1];
					double vvv = LinearInter(E, F, latA, latB, lat);
					result[ilat, ilon] = (float)vvv;
				}
			}
			return result;
		}
		private static float AbsDist(float x0, float y0, float x1, float y1)
		{
			return Math.Abs(x0 - x1) + Math.Abs(y0 - y1);
		}
		public static float[,] ReGridNearest(double[,] data, double[] latsSource, double[] lonsSource, float[] newlats, float[] newlons)
		{
			int nlat = newlats.Length;
			int nlon = newlons.Length;
           
			float[,] result = new float[nlat, nlon];

			int[] lonIndice = new int[nlon];
			float dlat = (float)(latsSource[1] - latsSource[0]);
			float dlon = (float)(lonsSource[1] - lonsSource[0]);
			
			for (int i = 0; i < nlon; i++) {
				float lon = newlons[i];                
				lonIndice[i] = (int)((lon - lonsSource[0]) / dlon);// Search(lon, lonsSource);
			}
			for (int ilat = 0; ilat < nlat; ilat++) {
				float lat = newlats[ilat];
				int latIndex = -(int)((lat - latsSource[0]) / dlat);// Search(lat, latsSource);
				//latIndex=latsSource.Length-1-latIndex;
				for (int ilon = 0; ilon < nlon; ilon++) {
					//Calculate bilinear interpolation using nearest four grid cells
					float lon = newlons[ilon];
					//lonIndex=Search(lon,lons1);
					int lonIndex = lonIndice[ilon];

					float A = (float)data[latIndex, lonIndex];
					float ADist = AbsDist(lat, lon, (float)latsSource[latIndex], (float)lonsSource[lonIndex]);
					float B = A;
					float BDist = ADist;
					if (lonIndex + 1 < lonsSource.Length) {
						B = (float)data[latIndex, lonIndex + 1];
						BDist = AbsDist(lat, lon, (float)latsSource[latIndex], (float)lonsSource[lonIndex + 1]);
					}
					float C = A;
					float	CDist = ADist;
					float D = A;
					float DDist = ADist;
					if (latIndex + 1 < latsSource.Length) {
						C = (float)data[latIndex + 1, lonIndex];
						CDist = AbsDist(lat, lon, (float)latsSource[latIndex + 1], (float)lonsSource[lonIndex]);
						if (lonIndex + 1 < lonsSource.Length) {
							D = (float)data[latIndex + 1, lonIndex + 1];
							DDist = AbsDist(lat, lon, (float)latsSource[latIndex + 1], (float)lonsSource[lonIndex + 1]);
						}
					}

					float[] aaa = new float[]{ A, B, C, D };
					float[] bbb = new float[]{ ADist, BDist, CDist, DDist };
					float min = 1e20f;
					float minValue = 0;
					for (int i = 0; i < 4; i++) {
						if (min > bbb[i]) {
							min = bbb[i];
							minValue = aaa[i];
						}
					}
					result[nlat - 1 - ilat, ilon] = minValue;				
					
				}
			}
			return result;
		}
		public static float[,] ReGridLambert(double[,] data, double[] latsSource, double[] lonsSource, float startLat, float startLon, float cellSize, int rowCount, int colCount
		                                     , float refLat, float refLon, float stdLat1, float stdLat2)
		{
			Array.Reverse(latsSource);
			float[,] result = new float[rowCount, colCount];
			MapProj.GaussKruger proj = new GaussKruger();
			double[] values0 = proj.ToXY_Lambert(startLon, startLat, refLon, refLat, stdLat1, stdLat2);
			float[] values=Array.ConvertAll<double,float>(values0,d=>(float)d);
			float lux = values[0];
			float luy = values[1];
			float X0 = 0;
			float Y0 = 0;
			for (int row = 0; row < rowCount; row++) {
				for (int col = 0; col < colCount; col++) {
					X0 = lux + cellSize * col;
					Y0 = luy - cellSize * row;
					values0= proj.ToLonLa_Lambert(X0, Y0, refLon, refLat, stdLat1, stdLat2);
					float[] vv=Array.ConvertAll<double,float>(values0,d=>(float)d);
					
					float lat = vv[1];
					float lon = vv[0];
					result[row, col] = (float)Interp(lat, lon, data, latsSource, lonsSource);
				}
			}
			return result;
		}
		private static float Interp(float curLat, float curLon, double[,] data, double[] latsSource, double[] lonsSource)
		{
			float dlat = (float)(latsSource[1] - latsSource[0]);
			float dlon = (float)(lonsSource[1] - lonsSource[0]);
			int latIndex = -(int)((curLat - latsSource[0]) / dlat);
			int lonIndex = (int)((curLon - lonsSource[0]) / dlon);
			if (dlat < 0) {
				
				dlat *= -1;
			}
			if (lonIndex + 1 >= lonsSource.Length) {
				return -9999;
			}
			float A = (float)data[latIndex, lonIndex];
			float B = (float)data[latIndex, lonIndex + 1];
			float C = (float)data[latIndex + 1, lonIndex];
			float D = (float)data[latIndex + 1, lonIndex + 1];
			float lonA = (float)lonsSource[lonIndex];
			float lonB = (float)lonsSource[lonIndex + 1];
			float	E = (float)LinearInter(A, B, lonA, lonB, curLon);
			float F = (float)LinearInter(C, D, lonA, lonB, curLon);
			float latA = (float)latsSource[latIndex];
			float latB = (float)latsSource[latIndex + 1];
			if (A < -9990 || B < -9990 || C < -9990 || D < -9990)
				return Math.Max(A, Math.Max(B, Math.Max(C, D)));
			return (float)LinearInter(E, F, latA, latB, curLat);
		}
		
		private static double LinearInter(double AValue, double BValue, double ACoor, double BCoor, double ECoor)
		{
			if (BCoor - ACoor == 0)
				return (AValue + BValue) / 2;
			return (AValue * (BCoor - ECoor) + BValue * (ECoor - ACoor)) / (BCoor - ACoor);
		}
		public static double[,] Smooth3x3(double[,] matrix)
		{
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			double[,] m = new double[rows, cols];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (i == 0 || j == 0 || i == rows - 1 || j == cols - 1)
						m[i, j] = matrix[i, j];
					else {
						m[i, j] = matrix[i, j] + matrix[i, j - 1] + matrix[i, j + 1] +
						matrix[i + 1, j] + matrix[i + 1, j - 1] + matrix[i + 1, j + 1] +
						matrix[i - 1, j] + matrix[i - 1, j - 1] + matrix[i - 1, j + 1];
						m[i, j] /= 9;
					}
				}
			}
			return m;
		}
       
		/// <summary>
		/// 从文本文件中读入所有数据行，数据行需要进一步分解为可用的数据
		/// </summary>
		/// <param name="FileName">文本文件名</param>
		/// <returns>数据行列表</returns>
		public static List<string> FileRead(string FileName)
		{
			//检查文件
			if (!File.Exists(FileName)) {
				throw new Exception("文件" + FileName + "不存在！");

			}
			ProgressBar bar = new ProgressBar();
			bar.Show();

			bar.Text = "正在读取数据";
			List<string> slist = new List<string>();
			StreamReader sr = new StreamReader(FileName);
			while (!sr.EndOfStream) {
				string s = sr.ReadLine();
				if (s != "")
					slist.Add(s);
			}
			sr.Close();
			bar.Close();
			return slist;
		}
        
		/// <summary>
		/// 将矩阵转化为黑白灰度位图
		/// </summary>
		/// <param name="matrix">要转化的矩阵</param>
		/// <param name="rows">总行数</param>
		/// <param name="cols">总列数</param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix(double[,] matrix)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = matrix[i, j];
					if (matrix[i, j] < min)
						min = matrix[i, j];
				}
			}
            
			double scope = max - min;
			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
					} else {
						byte value = (byte)((matrix[i, j] - min) * 255.0 / scope);
						rgb[(i * cols + j) * 4] = value;
						rgb[(i * cols + j) * 4 + 1] = value;
						rgb[(i * cols + j) * 4 + 2] = value;
						rgb[(i * cols + j) * 4 + 3] = 255;
					}
                
				}
			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);
			return bitmap;
		}
		public static Bitmap BitmapFromMatrix(int[,] matrix)
		{

			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = (int)matrix[i, j];
					if (matrix[i, j] < min)
						min = (int)matrix[i, j];
				}
			}
           
			double scope = max - min;
			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
					} else {
						byte value = (byte)((matrix[i, j] - min) * 255.0 / scope);

						rgb[(i * cols + j) * 4] = value;
						rgb[(i * cols + j) * 4 + 1] = value;
						rgb[(i * cols + j) * 4 + 2] = value;
						rgb[(i * cols + j) * 4 + 3] = 255;
					}
				}
			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);
			return bitmap;
		}
		/// <summary>
		/// 将矩阵转化为黑-蓝-绿-红-白渐变的彩色图像
		/// </summary>
		/// <param name="matrix"></param>
		/// <param name="rows"></param>
		/// <param name="cols"></param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix5(double[,] matrix)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = matrix[i, j];
					if (matrix[i, j] < min)
						min = matrix[i, j];
				}
			}
            
			double scope = max - min;

			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {                
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					double gray = (matrix[i, j] - min) / scope * 4;
                    
					byte r, g, b;
					if (gray < 1) {
						r = 0;
						g = 0;
						b = (byte)(gray * 255);
					} else if (gray < 2) {
						gray -= 1.0;
						r = 0;
						g = (byte)(gray * 255);
						b = (byte)((1 - gray) * 255);
					} else if (gray < 3) {
						gray -= 2.0;
						r = (byte)(gray * 255);
						g = (byte)((1 - gray) * 255);
						b = 0;
					} else {
						gray -= 3.0;
						r = 255;
						g = (byte)(gray * 255);
						b = (byte)(gray * 255);
					}
					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}
                
			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);
            
			return bitmap;
		}
        
		/// <summary>
		/// 将矩阵转化为黑-蓝-绿-红-白渐变的彩色图像
		/// </summary>
		/// <param name="matrix"></param>
		/// <param name="rows"></param>
		/// <param name="cols"></param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix5(int[,] matrix)
		{

			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = matrix[i, j];
					if (matrix[i, j] < min)
						min = matrix[i, j];
				}
			}

			double scope = max - min;

			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					double gray = (matrix[i, j] - min) / scope * 4;

					byte r, g, b;
					if (gray < 1) {
						r = 0;
						g = 0;
						b = (byte)(gray * 255);
					} else if (gray < 2) {
						gray -= 1.0;
						r = 0;
						g = (byte)(gray * 255);
						b = (byte)((1 - gray) * 255);
					} else if (gray < 3) {
						gray -= 2.0;
						r = (byte)(gray * 255);
						g = (byte)((1 - gray) * 255);
						b = 0;
					} else {
						gray -= 3.0;
						r = 255;
						g = (byte)(gray * 255);
						b = (byte)(gray * 255);
					}
					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}

			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);

			return bitmap;

		}
		/// <summary>
		/// 生成黑-蓝-绿-黄-红-白渐变的彩色
		/// </summary>
		/// <param name="min"></param>
		/// <param name="max"></param>
		/// <param name="value"></param>
		/// <returns></returns>
		public static Color ToColor6(double min, double max, double value)
		{
			if (max <= min)
				return Color.Black;
			if (max < min)
				throw new Exception("Max value must larger than min value:min=" + min + " max=" + max);
			value = (value > max) ? max : value;
			value = (value < min) ? min : value;
			double scope = max - min;
			double gray = (value - min) / scope * 5;
                    
			byte r, g, b;
			if (gray < 1) { //从黑至蓝
				r = 0;
				g = 0;
				b = (byte)(gray * 255);    
			} else if (gray < 2) {   //从蓝到绿
				gray -= 1.0;
				r = 0;
				g = (byte)(gray * 255);
				b = (byte)((1 - gray) * 255);
			} else if (gray < 3) {   //从绿至黄
				gray -= 2.0;
				r = (byte)(gray * 255);
				g = 255;
				b = 0;
			} else if (gray < 4) {  //从黄至红
				gray -= 3.0;
				r = 255;
				g = (byte)((1 - gray) * 255);
				b = 0; 
			} else {
				gray -= 4.0;
				r = 255;
				g = (byte)(gray * 255);
				b = (byte)(gray * 255);
			}
			return Color.FromArgb(r, g, b);
		}
		/// <summary>
		/// 将矩阵转化为黑-蓝-绿-黄-红-白渐变的彩色图像
		/// </summary>
		/// <param name="matrix"></param>
		/// <param name="rows"></param>
		/// <param name="cols"></param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix6(double[,] matrix)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = matrix[i, j];
					if (matrix[i, j] < min)
						min = matrix[i, j];
				}
			}            
			return BitmapFromMatrix6(matrix, max, min);        
		}
		public static Bitmap BitmapFromMatrix_Pallete(double[,] matrix, double max, double min,int count)
		{
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);            
			double scope = max - min;
			double step=scope/count;
			Color[] colors=new Color[count];
			for(int i=0;i<count;i++)
			{
				colors[i]= ToColor6(min,max,step*i+min);
			}
						
			Bitmap bitmap = new Bitmap(cols, rows);
			
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {                
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9998) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					
					int index=(int)((matrix[i,j]-min)/step);
					if(index<0)
					{
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					if(index>=count)
						index=count-1;
					Color color =colors[index];
					
					rgb[(i * cols + j) * 4] =color.B;
					rgb[(i * cols + j) * 4 + 1] = color.G;
					rgb[(i * cols + j) * 4 + 2] = color.R;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}                
			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);
            
			return bitmap;
		}
		public static Bitmap BitmapFromMatrix6(double[,] matrix, double max, double min)
		{
			//找出最大值和最小值
           
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);            
			double scope = max - min;

			Bitmap bitmap = new Bitmap(cols, rows);
			
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {                
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9998) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					double gray = (matrix[i, j] - min) / scope * 5;
                    
					byte r, g, b;
					if (scope == 0) {
						r = 0;
						g = 0;
						b = 0;
					} else {
						if (gray < 1) { //从黑至蓝
							r = 0;
							g = 0;
							b = (byte)(gray * 255);    
						} else if (gray < 2) {   //从蓝到绿
							gray -= 1.0;
							r = 0;
							g = (byte)(gray * 255);
							b = (byte)((1 - gray) * 255);
						} else if (gray < 3) {   //从绿至黄
							gray -= 2.0;
							r = (byte)(gray * 255);
							g = 255;
							b = 0;
						} else if (gray < 4) {  //从黄至红
							gray -= 3.0;
							r = 255;
							g = (byte)((1 - gray) * 255);
							b = 0; 
						} else {
							gray -= 4.0;
							r = 255;
							g = (byte)(gray * 255);
							b = (byte)(gray * 255);
						}
					}
					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}                
			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);
            
			return bitmap;
        
		}
        
		/// <summary>
		/// 将矩阵转化为黑-蓝-绿-黄-红-白渐变的彩色图像
		/// </summary>
		/// <param name="matrix"></param>
		/// <param name="rows"></param>
		/// <param name="cols"></param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix6(int[,] matrix)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = matrix[i, j];
					if (matrix[i, j] < min)
						min = matrix[i, j];
				}
			}

			double scope = max - min;

			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					double gray = (matrix[i, j] - min) / scope * 5;

					byte r, g, b;
					if (gray < 1) { //从黑至蓝
						r = 0;
						g = 0;
						b = (byte)(gray * 255);
					} else if (gray < 2) {   //从蓝到绿
						gray -= 1.0;
						r = 0;
						g = (byte)(gray * 255);
						b = (byte)((1 - gray) * 255);
					} else if (gray < 3) {   //从绿至黄
						gray -= 2.0;
						r = (byte)(gray * 255);
						g = 255;
						b = 0;
					} else if (gray < 4) {  //从黄至红
						gray -= 3.0;
						r = 255;
						g = (byte)((1 - gray) * 255);
						b = 0;
					} else {
						gray -= 4.0;
						r = 255;
						g = (byte)(gray * 255);
						b = (byte)(gray * 255);
					}
					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}

			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);

			return bitmap;

		}
		/// <summary>
		/// 将矩阵转变为三色混合图像
		/// </summary>
		/// <param name="matrix"></param>
		/// <param name="color1"></param>
		/// <param name="color2"></param>
		/// <param name="color3"></param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix3(double[,] matrix, Color color1, Color color2, Color color3)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = (int)matrix[i, j];
					if (matrix[i, j] < min)
						min = (int)matrix[i, j];
				}
			}
			byte r1 = color1.R;
			byte g1 = color1.G;
			byte b1 = color1.B;
			byte r2 = color2.R;
			byte g2 = color2.G;
			byte b2 = color2.B;
			byte r3 = color3.R;
			byte g3 = color3.G;
			byte b3 = color3.B;
			double scope = max - min;
			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					byte r, g, b;
					double gray = (matrix[i, j] - min) / scope * 2;
					if (gray < 1) {
						r = (byte)(gray * (r2 - r1) + r1);
						g = (byte)(gray * (g2 - g1) + g1);
						b = (byte)(gray * (b2 - b1) + b1);
					} else {
						gray -= 1;
						r = (byte)(gray * (r3 - r2) + r2);
						g = (byte)(gray * (g3 - g2) + g2);
						b = (byte)(gray * (b3 - b2) + b2);
					}

					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}

			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);

			return bitmap;
		}
		public static Bitmap BitmapFromMatrix3(int[,] matrix, Color color1, Color color2, Color color3)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = (int)matrix[i, j];
					if (matrix[i, j] < min)
						min = (int)matrix[i, j];
				}
			}
			byte r1 = color1.R;
			byte g1 = color1.G;
			byte b1 = color1.B;
			byte r2 = color2.R;
			byte g2 = color2.G;
			byte b2 = color2.B;
			byte r3 = color3.R;
			byte g3 = color3.G;
			byte b3 = color3.B;
			double scope = max - min;
			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					byte r, g, b;
					double gray = (matrix[i, j] - min) / scope * 2;
					if (gray < 1) {
						r = (byte)(gray * 2 * (r2 - r1) + r1);
						g = (byte)(gray * 2 * (g2 - g1) + g1);
						b = (byte)(gray * 2 * (b2 - b1) + b1);
					} else {
						gray -= 1;
						r = (byte)(gray * 2 * (r3 - r2) + r2);
						g = (byte)(gray * 2 * (g3 - g2) + g2);
						b = (byte)(gray * 2 * (b3 - b2) + b2);
					}

					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}

			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);

			return bitmap;
		}
		/// <summary>
		/// 将矩阵转变为二色混合图像
		/// </summary>
		/// <param name="matrix"></param>
		/// <param name="color1"></param>
		/// <param name="color2"></param>
		/// <returns></returns>
		public static Bitmap BitmapFromMatrix2(double[,] matrix, Color color1, Color color2)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = (int)matrix[i, j];
					if (matrix[i, j] < min)
						min = (int)matrix[i, j];
				}
			}
			byte r1 = color1.R;
			byte g1 = color1.G;
			byte b1 = color1.B;
			byte r2 = color2.R;
			byte g2 = color2.G;
			byte b2 = color2.B;
           
			double scope = max - min;
			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					byte r, g, b;
					double gray = (matrix[i, j] - min) / scope;
					r = (byte)(gray * (r2 - r1) + r1);
					g = (byte)(gray * (g2 - g1) + g1);
					b = (byte)(gray * (b2 - b1) + b1);

					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}

			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);

			return bitmap;
		}
		public static Bitmap BitmapFromMatrix2(int[,] matrix, Color color1, Color color2)
		{
			//找出最大值和最小值
			double max = -99999999;
			double min = 99999999;
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990)
						continue;
					if (matrix[i, j] > max)
						max = (int)matrix[i, j];
					if (matrix[i, j] < min)
						min = (int)matrix[i, j];
				}
			}
			byte r1 = color1.R;
			byte g1 = color1.G;
			byte b1 = color1.B;
			byte r2 = color2.R;
			byte g2 = color2.G;
			byte b2 = color2.B;
            
			double scope = max - min;
			Bitmap bitmap = new Bitmap(cols, rows);
			BitmapData bmData = bitmap.LockBits(new Rectangle(0, 0, cols, rows), ImageLockMode.ReadWrite, PixelFormat.Format32bppArgb);
			IntPtr ptr = bmData.Scan0;
			int n = bmData.Stride * bitmap.Height;
			byte[] rgb = new byte[n];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] < -9990) {
						rgb[(i * cols + j) * 4 + 3] = 0;
						continue;
					}
					byte r, g, b;
					double gray = (matrix[i, j] - min) / scope;
					r = (byte)(gray * (r2 - r1) + r1);
					g = (byte)(gray * (g2 - g1) + g1);
					b = (byte)(gray * (b2 - b1) + b1);

					rgb[(i * cols + j) * 4] = b;
					rgb[(i * cols + j) * 4 + 1] = g;
					rgb[(i * cols + j) * 4 + 2] = r;
					rgb[(i * cols + j) * 4 + 3] = 255;
				}

			}
			System.Runtime.InteropServices.Marshal.Copy(rgb, 0, ptr, n);
			bitmap.UnlockBits(bmData);

			return bitmap;
		}
		//深度复制一个矩阵
		public static double[,] CopyMatrix(double[,] matrix)
		{
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			double[,] m = new double[rows, cols];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					m[i, j] = matrix[i, j];
				}
			}
			return m;
		}
		//深度复制一个矩阵
		public static int[,] CopyMatrix(int[,] matrix)
		{
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			int[,] m = new int[rows, cols];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					m[i, j] = matrix[i, j];
				}
			}
			return m;
		}
		public static void SaveMatrix(int[,] matrix, string FileName)
		{
			StreamWriter sw = new StreamWriter(FileName);
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					sw.Write(matrix[i, j].ToString("0.00") + "\t");
				}
				sw.WriteLine();
			}
			sw.Close();
		}
		public static void SaveMatrix(double[,] matrix, string FileName)
		{
			StreamWriter sw = new StreamWriter(FileName);
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					sw.Write(matrix[i, j].ToString("0.00") + "\t");
				}
				sw.WriteLine();
			}
			sw.Close();
		}
		public static double[,] LogTransform(double[,] matrix)
		{
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			double[,] ma = new double[rows, cols];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] <= 0)
						ma[i, j] = -9999;
					else
						ma[i, j] = Math.Log(matrix[i, j]);
				}
			}
			return ma;
		}
		public static double[,] LogTransform(int[,] matrix)
		{
			int rows = matrix.GetLength(0);
			int cols = matrix.GetLength(1);
			double[,] ma = new double[rows, cols];
			for (int i = 0; i < rows; i++) {
				for (int j = 0; j < cols; j++) {
					if (matrix[i, j] <= 0)
						ma[i, j] = -9999;
					else
						ma[i, j] = Math.Log(matrix[i, j]);
				}
			}
			return ma;
		}
		public static double[,] A = new double[3, 3] {
			{ 0.567, 0.086, -0.002 }, {
				0.253,
				0.504,
				-0.050
			}, {
				-0.006,
				-0.039,
				0.244
			}
		};
		public static double[,] B = new double[3, 3] {
			{ 0.781, 0, 0 }, {
				0.328,
				0.637,
				0
			}, {
				0.238,
				-0.341,
				0.873
			}
		};
	}
}