#include "Feature_extraction.h"

//濃度特徴
cv::Mat Feature_extraction::Density_histgram(Raw_image<short>& ct, array3<unsigned int>& labeling, size_t labelcount)
{

	//ヒストグラム作成　-100〜 1200;
	auto r_Hist = std::vector<DensityData>(labelcount + 1);
	for (int L = 1; L < labelcount + 1; L++) {

		r_Hist[L].Hist = new int[1300];
		r_Hist[L].NormalizeHist = new double[1300];
		r_Hist[L].CPD = new double[1300];


		for (int i = 0; i < 1300; i++) {
			r_Hist[L].Hist[i] = 0;
			r_Hist[L].NormalizeHist[i] = 0;
		}
	}

	for (int i = 0; i < ct.data.VOXELS; i++) {
		if (labeling[i] != 0) {

			if (ct.image[i] < -100) {
				r_Hist[labeling[i]].Hist[0] ++;
			}
			else if (ct.image[i] > 1199) {
				r_Hist[labeling[i]].Hist[1299] ++;
			}
			else
				r_Hist[labeling[i]].Hist[ct.image[i] + 100] ++;



			if (ct.image[i] > r_Hist[labeling[i]].Max || r_Hist[labeling[i]].Max == -999) {
				if (ct.image[i] < 1299)
					r_Hist[labeling[i]].Max = ct.image[i];
				else
					r_Hist[labeling[i]].Max = 1299;

			}

			if (ct.image[i] < r_Hist[labeling[i]].Min || r_Hist[labeling[i]].Min == -999)
				if (ct.image[i] > -100)
					r_Hist[labeling[i]].Min = ct.image[i];
				else
					r_Hist[labeling[i]].Min = -100;

			r_Hist[labeling[i]].HistVOXELS++;
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < 1300; i++) {
			if (i == 0 || i == 1299) {
				r_Hist[k].NormalizeHist[i] = 0;
			}
			else {
				r_Hist[k].NormalizeHist[i] = r_Hist[k].Hist[i] / (double)r_Hist[k].HistVOXELS;
			}

			if (i != 0)
				r_Hist[k].CPD[i] += r_Hist[k].CPD[i - 1] + r_Hist[k].NormalizeHist[i];
			else
				r_Hist[k].CPD[i] = 0;
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		r_Hist[k].Average = 0;
		for (int i = 1; i < 1300; i++) {
			r_Hist[k].Average += (i - 100) * r_Hist[k].Hist[i];
		}
		r_Hist[k].Average = r_Hist[k].Average / r_Hist[k].HistVOXELS;

	}


	for (int k = 1; k < labelcount + 1; k++) {
		r_Hist[k].Variance = 0;
		for (int i = 1; i < 1300; i++) {
			r_Hist[k].Variance += (double)std::pow(((i - 100) - r_Hist[k].Average), 2) * r_Hist[k].NormalizeHist[i];
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		r_Hist[k].Skewness = 0;
		r_Hist[k].Kurtosis = 0;
		for (int i = 1; i < 1300; i++) {
			r_Hist[k].Skewness += (double)std::pow(((i - 100) - r_Hist[k].Average), 3) * r_Hist[k].NormalizeHist[i];
			r_Hist[k].Kurtosis += (double)std::pow(((i - 100) - r_Hist[k].Average), 4) * r_Hist[k].NormalizeHist[i];
		}
		if (r_Hist[k].Variance == 0) {
			r_Hist[k].Skewness = 0;
			r_Hist[k].Kurtosis = 2;
		}
		else {
			r_Hist[k].Skewness = r_Hist[k].Skewness / (double)std::pow(r_Hist[k].Variance, 1.5);
			r_Hist[k].Kurtosis = r_Hist[k].Kurtosis / (double)std::pow(r_Hist[k].Variance, 2);
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		delete[] r_Hist[k].CPD;
		delete[] r_Hist[k].Hist;
		delete[] r_Hist[k].NormalizeHist;
	}

	cv::Mat samples;
	samples.create((int)labelcount, 7, CV_32F);

	for (int i = 0; i < labelcount; i++) {
		samples.at<float>(i, 0) = (float)r_Hist[i + 1].HistVOXELS;
		samples.at<float>(i, 1) = (float)r_Hist[i + 1].Max;
		samples.at<float>(i, 2) = (float)r_Hist[i + 1].Min;
		samples.at<float>(i, 3) = (float)r_Hist[i + 1].Average;
		samples.at<float>(i, 4) = (float)r_Hist[i + 1].Variance;
		samples.at<float>(i, 5) = (float)r_Hist[i + 1].Skewness;
		samples.at<float>(i, 6) = (float)r_Hist[i + 1].Kurtosis;
	}

	return (samples);
}

//テクスチャ特徴
array3<uchar> Feature_extraction::Histogram_equalization(Raw_image<short>& ct, int len)
{
	//ヒストグラム作成　-100〜 1200;
	array3<uchar> gets(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0);
	unsigned int* s_Hist = new unsigned int[1300];
	unsigned int* s_CHist = new unsigned int[1300];
	double s_VOXELS = 0;
	uchar* FrattarHist = new uchar[1300];

	//ヒストグラム平坦化
#pragma omp parallel for
	for (int i = 0; i < 1300; i++) {
		s_Hist[i] = 0;
		s_CHist[i] = 0;
	}

	for (int i = 0; i < ct.data.VOXELS; i++) {
		if (ct.image[i] < -100) {
			s_Hist[0] ++;
		}
		else if (ct.image[i] > 1199) {
			s_Hist[1299] ++;
		}
		else {
			s_Hist[ct.image[i] + 100] ++;
			s_VOXELS++;
		}
	}

	for (int i = 0; i < 1300; i++) {
		if (i == 0)
			s_CHist[i] = 0;
		else if (i == 1299)
			s_CHist[i] = s_CHist[i - 1];
		else
			s_CHist[i] = s_CHist[i - 1] + s_Hist[i];

	}
	int Counta = -999;


	for (int i = 0; i < 1300; i++) {

		if (i == 0)
			FrattarHist[i] = 0;
		else if (i == 1299)
			FrattarHist[i] = len - 1;
		else if ((((double)s_CHist[i] / s_VOXELS * (double)len) - 1) < 0)
			FrattarHist[i] = 0;
		else
			FrattarHist[i] = (uchar)(((double)s_CHist[i] / s_VOXELS * (double)len) - 1);

		if (FrattarHist[i] < 0) {
			FrattarHist[i] = 0;
		}

		if (Counta != FrattarHist[i]) {
			Counta = FrattarHist[i];
		}

	}



#pragma omp parallel for
	for (int i = 0; i < ct.data.VOXELS; i++) {
		if (ct.image[i] < -100)
			gets[i] = 0;
		else if (ct.image[i] > 1199)
			gets[i] = len - 1;
		else
			gets[i] = FrattarHist[ct.image[i] + 100];
	}

	delete[] s_Hist;
	delete[] s_CHist;
	delete[] FrattarHist;
	return(gets);
}
cv::Mat Feature_extraction::GLCM(array3<uchar>& image, array3<unsigned int>& label, size_t& labelcount,int len)
{
	std::vector<GLCMData> s_Sim(labelcount + 1);

	array3<uchar> region(image.width(), image.height(), image.depth(), image.reso1(), image.reso1(), image.reso1(), 0);
	array3<unsigned int> Labeling(image.width(), image.height(), image.depth(), image.reso1(), image.reso1(), image.reso1(), 0);
	for (int i = 0; i < image.size(); i++) {
		Labeling[i] = label[i];
		if (label[i] != 0)
			region[i] = 255;
	}
	//mist::erosion(region,2);
	//for (int i = 0; i < Data.VOXELS; i++) {
	//	if (region[i] != 0)
	//		Labeling[i] = 0;
	//}
	int ImageSize = image.width() * image.height();
	//構造要素の生成
	int FilterStructure[27];
	int StructureCount = 0;
	for (int dz = -1; dz <= 1; dz++) {
		for (int dy = -1; dy <= 1; dy++)
			for (int dx = -1; dx <= 1; dx++) {
				FilterStructure[StructureCount] = ImageSize * dz + image.width() * dy + dx;
				StructureCount++;
			}
	}
	int pass = ((3 * 3 * 3) / 2);

	for (int k = 0; k < labelcount + 1; k++) {
		s_Sim[k].simMatrinx = new float* [len];
		s_Sim[k].Q = new float* [len];


		for (int i = 0; i < len; i++) {
			s_Sim[k].simMatrinx[i] = new float[len];
			s_Sim[k].Q[i] = new float[len];
		}
	}

	for (int k = 0; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].simMatrinx[i][j] = 0;
				s_Sim[k].Q[i][j] = 0;
			}
	}

	//同時生起行列計算
	for (int z = 1; z < image.depth() - 1; z++)
		for (int y = 1; y < image.height() - 1; y++)
			for (int x = 1; x < image.width() - 1; x++) {
				int point = ImageSize * z + y * image.width() + x;

				if (Labeling[point] != 0) {
					for (int i = 0; i < 27; i++) {
						if (i != pass) {
							s_Sim[Labeling[point]].MatrixCount++;
							s_Sim[Labeling[point]].simMatrinx[image[point]][image[point + FilterStructure[i]]] ++;
						}
					}
				}
			}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].simMatrinx[i][j] = s_Sim[k].simMatrinx[i][j] / s_Sim[k].MatrixCount;
			}
	}

	//特徴量抽出処
	for (int k = 1; k < labelcount + 1; k++) {
		s_Sim[k].Py = new float[len];
		s_Sim[k].Pxsy = new float[len];
		s_Sim[k].Px = new float[len];
		s_Sim[k].Pxpy = new float[2 * len - 1];
	}

	//初期化
	for (int k = 1; k < labelcount + 1; k++) {
		for (int j = 0; j < len; j++) {
			s_Sim[k].Py[j] = 0;
			s_Sim[k].Px[j] = 0;
			s_Sim[k].Pxsy[j] = 0;
		}
	}


	for (int k = 1; k < labelcount + 1; k++) {
		for (int j = 0; j < 2 * len - 1; j++) {
			s_Sim[k].Pxpy[j] = 0;
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].Py[i] += s_Sim[k].simMatrinx[j][i];
				s_Sim[k].Px[i] += s_Sim[k].simMatrinx[i][j];
				s_Sim[k].Pxpy[i + j] += s_Sim[k].simMatrinx[i][j];
				s_Sim[k].Pxsy[abs(i - j)] += s_Sim[k].simMatrinx[i][j];
			}
	}

	//全14種
	//1.角度の二次モーメント
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].ASM += pow(s_Sim[k].simMatrinx[i][j], 2);
			}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			s_Sim[k].DASM += pow(s_Sim[k].Pxsy[i], 2);

	}

	//2.コントラスト
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++) {
			s_Sim[k].Contrast += i * i * s_Sim[k].Pxsy[i];
		}
	}

	//3.相関
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++) {
			s_Sim[k].Ax += i * s_Sim[k].Px[i];
			s_Sim[k].Ay += i * s_Sim[k].Py[i];
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++) {
			s_Sim[k].Sx += pow(i - s_Sim[k].Ax, 2) * s_Sim[k].Px[i];
			s_Sim[k].Sy += pow(i - s_Sim[k].Ay, 2) * s_Sim[k].Py[i];
		}
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].Correlation += i * j * s_Sim[k].simMatrinx[i][j] - s_Sim[k].Ax * s_Sim[k].Ay;
			}

		if (s_Sim[k].Sx * s_Sim[k].Sy != 0)
			s_Sim[k].Correlation = s_Sim[k].Correlation / (s_Sim[k].Sx * s_Sim[k].Sy);
	}

	//4.分散
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].Variance += pow(i - s_Sim[k].Ax, 2) * s_Sim[k].simMatrinx[j][i];
			}
	}

	//5.差分モーメント
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				s_Sim[k].IDM += (float)(1 / (1 + pow(i - j, 2))) * s_Sim[k].simMatrinx[j][i];
			}
	}

	//6.平均
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < 2 * len - 1; i++)
			s_Sim[k].SAverage += i * s_Sim[k].Pxpy[i];
	}

	//7.分散
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < 2 * len - 1; i++)
			s_Sim[k].SVariance += pow(i - s_Sim[k].SAverage, 2) * s_Sim[k].Pxpy[i];
	}
	//8.合計エントロピー
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < 2 * len - 1; i++) {
			if (s_Sim[k].Pxpy[i] > 0) {
				s_Sim[k].SEntropy += s_Sim[k].Pxpy[i] * log2(s_Sim[k].Pxpy[i]);
			}
		}
		s_Sim[k].SEntropy = s_Sim[k].SEntropy * -1;
	}

	//9.エントロピー
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				if (s_Sim[k].simMatrinx[j][i] > 0) {
					s_Sim[k].Entropy += s_Sim[k].simMatrinx[j][i] * log2(s_Sim[k].simMatrinx[j][i]);
				}
			}

		s_Sim[k].Entropy = s_Sim[k].Entropy * -1;
	}

	//差平均
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			s_Sim[k].SAvey += i * s_Sim[k].Pxsy[i];
	}


	//10.差分散
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			s_Sim[k].DVariance += pow(i - s_Sim[k].SAvey, 2) * s_Sim[k].Pxsy[i];
	}

	//11.差エントロピー
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			if (s_Sim[k].Pxsy[i] > 0) {
				s_Sim[k].DEntropy += s_Sim[k].Pxsy[i] * log2(s_Sim[k].Pxsy[i]);
			}
		s_Sim[k].DEntropy = s_Sim[k].DEntropy * -1;
	}


	//12.IMOC1
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			if (s_Sim[k].Px[i] > 0) {
				s_Sim[k].HX += s_Sim[k].Px[i] * log2(s_Sim[k].Px[i]);
			}
		s_Sim[k].HX = s_Sim[k].HX * -1;
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			if (s_Sim[k].Py[i] > 0)
				s_Sim[k].HY += s_Sim[k].Py[i] * log2(s_Sim[k].Py[i]);

		s_Sim[k].HY = s_Sim[k].HY * -1;
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				if (s_Sim[k].Px[j] * s_Sim[k].Py[i] > 0)
					s_Sim[k].HXY1 += s_Sim[k].simMatrinx[j][i] * log2(s_Sim[k].Px[j] * s_Sim[k].Py[i]);
			}
		s_Sim[k].HXY1 = s_Sim[k].HXY1 * -1;
	}

	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++)
			for (int j = 0; j < len; j++) {
				if (s_Sim[k].Px[j] * s_Sim[k].Py[i] > 0)
					s_Sim[k].HXY2 += s_Sim[k].Px[j] * s_Sim[k].Py[i] * log2(s_Sim[k].Px[j] * s_Sim[k].Py[i]);
			}
		s_Sim[k].HXY2 = s_Sim[k].HXY2 * -1;
	}

	for (int k = 1; k < labelcount + 1; k++) {
		double MAX = s_Sim[k].HX > s_Sim[k].HY ? s_Sim[k].HX : s_Sim[k].HY;
		if (MAX != 0)
			s_Sim[k].IMOC1 = (s_Sim[k].Entropy - s_Sim[k].HXY1) / (float)MAX;
		else
			s_Sim[k].IMOC1 = 0;
	}

	//13.IMOC2
	for (int k = 1; k < labelcount + 1; k++) {
		s_Sim[k].IMOC2 = (float)pow(1 - exp(-2.0 * (s_Sim[k].HXY2 - s_Sim[k].Entropy)), 0.5);
	}

	//14.MCC
	for (int k = 1; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++) {
			for (int j = 0; j < len; j++) {
				for (int s = 0; s < len; s++) {
					if (s_Sim[k].Px[i] * s_Sim[k].Py[j] > 0)
						s_Sim[k].Q[i][j] += s_Sim[k].simMatrinx[i][s] * s_Sim[k].simMatrinx[s][j] / (s_Sim[k].Px[i] * s_Sim[k].Py[j]);
				}

				if (s_Sim[k].Q[i][j] > s_Sim[k].Mf || s_Sim[k].Mf == -999) {
					if (s_Sim[k].Mn == -999) {
						s_Sim[k].Mf = s_Sim[k].Q[i][j];
						s_Sim[k].Mn = s_Sim[k].Q[i][j];
					}
					else {
						s_Sim[k].Mn = s_Sim[k].Mf;
						s_Sim[k].Mf = s_Sim[k].Q[i][j];
					}
				}
			}

		}
	}


	for (int k = 1; k < labelcount + 1; k++) {
		s_Sim[k].MCC = (float)pow(s_Sim[k].Mn, 0.5);
	}


	cv::Mat samples;
	samples.create((int)labelcount, 16, CV_32F);

	for (int i = 0; i < labelcount; i++) {
		samples.at<float>(i, 0) = (float)s_Sim[i + 1].ASM;
		samples.at<float>(i, 1) = (float)s_Sim[i + 1].Contrast;
		samples.at<float>(i, 2) = (float)s_Sim[i + 1].Correlation;
		samples.at<float>(i, 3) = (float)s_Sim[i + 1].Variance;
		samples.at<float>(i, 4) = (float)s_Sim[i + 1].IDM;
		samples.at<float>(i, 5) = (float)s_Sim[i + 1].SAverage;
		samples.at<float>(i, 6) = (float)s_Sim[i + 1].SVariance;
		samples.at<float>(i, 7) = (float)s_Sim[i + 1].SEntropy;
		samples.at<float>(i, 8) = (float)s_Sim[i + 1].Entropy;
		samples.at<float>(i, 9) = (float)s_Sim[i + 1].DVariance;
		samples.at<float>(i, 10) = (float)s_Sim[i + 1].DEntropy;
		samples.at<float>(i, 11) = (float)s_Sim[i + 1].IMOC1;
		samples.at<float>(i, 12) = (float)s_Sim[i + 1].IMOC2;
		samples.at<float>(i, 13) = (float)s_Sim[i + 1].MCC;
		samples.at<float>(i, 14) = (float)s_Sim[i + 1].DASM;
		samples.at<float>(i, 15) = (float)s_Sim[i + 1].SAvey;
	}

	for (int k = 0; k < labelcount + 1; k++) {
		for (int i = 0; i < len; i++) {
			delete[] s_Sim[k].simMatrinx[i];
			delete[] s_Sim[k].Q[i];
		}
		delete[] s_Sim[k].simMatrinx;
		delete[] s_Sim[k].Q;
		delete[] s_Sim[k].Px;
		delete[] s_Sim[k].Py;
		delete[] s_Sim[k].Pxpy;
		delete[] s_Sim[k].Pxsy;
	}

	return samples;
}

//勾配特徴
cv::Mat Feature_extraction::Laplace_feature(Raw_image<short>& ct, array3<unsigned int>& label, size_t& labelcount)
{
	mist::array3<short> CT(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, ct.image, ct.data.VOXELS);
	mist::array3<float> CTb, laplace;
	mist::gaussian_filter(CT, CTb);
	mist::laplacian_filter(CTb, laplace);
	auto r_Hist = std::vector<Extravasation_detect::LabelData>(labelcount + 1);
	for (int i = 0; i < ct.data.VOXELS; i++)
		if (label[i] != 0) {
			r_Hist[label[i]].Area++;
			r_Hist[label[i]].CTValue += laplace[i];
		}

	for (int i = 1; i < labelcount+ 1; i++) {
		r_Hist[i].CTValue = r_Hist[i].CTValue / r_Hist[i].Area;
	}

	cv::Mat samples;
	samples.create((int)labelcount, 1, CV_32F);

	for (int i = 0; i < labelcount; i++) {
		samples.at<float>(i, 0) = (float)r_Hist[i + 1].CTValue;
	}

	return (samples);
}

//Matの操作
cv::Mat Feature_extraction::MatLink(cv::Mat& firstMat, cv::Mat& secondMat)
{
	if (firstMat.rows != secondMat.rows) {
		printf("サンプル数が等しくありません\n");
		return cv::Mat();
	}

	int Firstcol = firstMat.cols;
	int Secondcol = secondMat.cols;

	cv::Mat samples;
	samples.create(firstMat.rows, firstMat.cols + secondMat.cols, CV_32F);

	for (int i = 0; i < firstMat.rows; i++) {
		for (int k = 0; k < Firstcol + Secondcol; k++) {
			if (k < Firstcol)
				samples.at<float>(i, k) = (float)firstMat.at<float>(i, k);
			else
				samples.at<float>(i, k) = (float)secondMat.at<float>(i, k - Firstcol);
		}
	}

	return samples;
}
cv::Mat Feature_extraction::Addsample(cv::Mat& firstMat, cv::Mat& secondMat)
{
	if (firstMat.cols != secondMat.cols) {
		printf("特徴量数が等しくありません.F:%d S:%d\n", firstMat.cols, secondMat.cols);
		return cv::Mat();
	}
	cv::Mat samples;
	samples.create(firstMat.rows + secondMat.rows, firstMat.cols, CV_32F);

	for (int i = 0; i < firstMat.rows + secondMat.rows; i++) {
		for (int k = 0; k < firstMat.cols; k++) {
			if (i < firstMat.rows)
				samples.at<float>(i, k) = (float)firstMat.at<float>(i, k);
			else
				samples.at<float>(i, k) = (float)secondMat.at<float>(i - firstMat.rows, k);
		}
	}
	return samples;
}
cv::Mat Feature_extraction::LabelMat(cv::Mat& input, int Label)
{
	cv::Mat samples;
	samples.create(input.rows, input.cols + 1, CV_32F);
	for (int i = 0; i < input.rows; i++) {
		for (int k = 0; k < input.cols + 1; k++) {
			if (k == 0)
				samples.at<float>(i, k) = (float)Label;
			else
				samples.at<float>(i, k) = (float)input.at<float>(i, k - 1);
		}
	}

	return samples;
}

//スケーリング
void Feature_extraction::Scaling(cv::Mat& feature)
{
	double mean[24] = { 1272.082474,228.1726804,46.04896907,136.4137956
	,1591.962119,0.013715459,2.886998773,0.023232985,40.59967549,-222509.8116,55.82264455,
	0.35392532,99.37873564,105.9691725,4.644675446,6.800864593,23.03388621,3.277484601,-0.114886773,0.755248402,0.494389606,
	0.158205928,3.749504322,-84.98969072 };
	double std[24] = { 4743.678168,105.5522472,39.9497129,47.78214929,2334.156954,0.454731201,0.919405181,0.020641736,35.97252044,1303192.787,
		45.77337404,0.116288456,9.207859944,99.59759599,0.674429203,1.190320582,18.6302418,0.624693816,0.038442077,0.087941346,
		0.095944749,0.066310616,1.87270045,59.92641447 };

	for (int i = 0; i < feature.rows; i++) {
		auto p = feature.row(i).ptr<float>(0);
		for (int j = 0; j < feature.cols; j++) {
			p[j] = (p[j] - mean[j]) / std[j];
		}
	}
}

//スクリーニング
cv::Mat Feature_extraction::Screening(cv::Mat& samples)
{
	int rows = samples.rows;
	int count = 0;
	std::vector<cv::Mat> s;
	for (int i = 0; i < rows; i++) {
		if (samples.at<float>(i, 0) >= 15) {
			count++;
			s.push_back(samples.row(i).clone());

		}
	}
	cv::Mat Return;
	Return.create(count, samples.cols, CV_32F);

	for (int i = 0; i < count; i++) {
		for (int j = 0; j < samples.cols; j++) {
			Return.at<float>(i, j) = s[i].at<float>(0, j);
		}
	}

	return Return;
}

//すべての特徴を一つのマットに保存する関数
cv::Mat Feature_extraction::All_feature_exstract(Raw_image<short> ct, array3<unsigned int>& label, size_t& labelcount, int len)
{
	auto HE = Histogram_equalization(ct, len);
	cv::Mat SimMat = GLCM(HE, label, labelcount,len);
	cv::Mat DensityMat = Density_histgram(ct, label, labelcount);
	cv::Mat FestureMat = MatLink(DensityMat, SimMat);
	cv::Mat laplaceMat = Laplace_feature(ct, label, labelcount);
	FestureMat = MatLink(FestureMat, laplaceMat);
	FestureMat = Screening(FestureMat);
	Scaling(FestureMat);
	return FestureMat;
}

//学習データの正解特徴と誤検出特徴の数の調整
cv::Mat Feature_extraction::Feature_Cut(cv::Mat& samples, bool shafful, float Falserate)
{
	int Label0 = 0;
	int Label1 = 0;
	for (int i = 0; i < samples.rows; i++) {
		if (samples.at<float>(i, 0) == 1) {
			Label1++;

		}
		else {
			Label0++;

		}
	}
	int Falsenum;

	if (Falserate == -1) {
		Falsenum = Label0;
	}
	else
		Falsenum = (int)(Label1 * Falserate);

	printf("Label1:%d,Label0:%d\n", Label1, Label0);
	cv::Mat TrueMat, FalseMat;
	TrueMat.create(Label1, samples.cols, CV_32F);
	FalseMat.create((int)(Falsenum), samples.cols, CV_32F);

	int TCount = 0;
	int FCount = 0;
	int Rate = (int)(Label0 / Label1 / Falserate);

	printf("発生Rate:%d\n", Rate);
	std::random_device rnd;     // 非決定的な乱数生成器を生成
	std::mt19937 mt(rnd());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
	std::uniform_int_distribution<> rand100(0, Rate - 1);

	for (int i = 0; i < samples.rows; i++) {
		if (samples.at<float>(i, 0) == 1) {
			if (TCount < Label1) {
				for (int k = 0; k < samples.cols; k++) {
					TrueMat.at<float>(TCount, k) = samples.at<float>(i, k);
				}
				TCount++;
			}

		}
		else {
			if (Falserate == -1) {
				if (FCount < FalseMat.rows) {
					for (int k = 0; k < samples.cols; k++) {
						FalseMat.at<float>(FCount, k) = samples.at<float>(i, k);
					}
					FCount++;
				}

			}
			else {
				if (rand100(mt) < 5) {
					if (FCount < FalseMat.rows) {
						for (int k = 0; k < samples.cols; k++) {
							FalseMat.at<float>(FCount, k) = samples.at<float>(i, k);
						}
						FCount++;
					}
				}

			}

		}
	}
	printf("選定結果:1:%d,0:%d\n", TCount, FCount);

	cv::Mat Return = Addsample(TrueMat, FalseMat);

	std::random_device rnda;     // 非決定的な乱数生成器を生成
	std::mt19937 mta(rnda());     //  メルセンヌ・ツイスタの32ビット版、引数は初期シード値
	std::uniform_int_distribution<> randes(0, Return.rows - 1); // [0, 99] 範囲の一様乱数

	if (shafful) {
		printf("シャッフル開始\n");
		for (int i = 0; i < Return.rows; i++) {
			int randres = randes(mta);
			cv::Mat tmp = Return.row(randres).clone();
			cv::Mat tmp1 = Return.row(i).clone();

			for (int k = 0; k < Return.cols; k++) {
				Return.at<float>(i, k) = tmp.at<float>(0, k);
				Return.at<float>(randres, k) = tmp1.at<float>(0, k);
			}
		}
	}

	return Return;
}

//学習したモデルを用いての分類
array3<uchar> Feature_extraction::Predict(cv::Ptr<cv::ml::RTrees> r_tree, array3<unsigned int>& label, size_t& labelcount, cv::Mat samples)
{
	float* mtest = new float[labelcount + 1];
	float* test = new float[labelcount + 1];

	for (int k = 1; k < labelcount + 1; k++) {
		cv::Mat sam = samples.row(k - 1).clone();
		cv::Mat bsam = samples.row(k - 1).clone();
		mtest[k] = r_tree->predict(sam, cv::noArray(), cv::ml::StatModel::RAW_OUTPUT);
		mtest[k] = mtest[k] / 1000;

		if (mtest[k] > 0.4) {
			test[k] = 1;
		}
		else
			test[k] = 0;
	}

	array3<uchar> s(label.width(), label.height(), label.depth(), 0.5, 0.5, 0.5, 0);

#pragma omp parallel for
	for (int i = 0; i < label.size(); i++) {
		if (label[i] != 0) {
			if (test[label[i]] == 0) {
				s[i] = 0;
			}
			else {
				s[i] = 255;
			}
		}
	}

	delete[] mtest;
	delete[] test;
	return (s);
}

//CSVの入出力
void Feature_extraction::writeCSV(char* FileName, cv::Mat samples)
{
	FILE* file;
	fopen_s(&file, FileName, "w");
	int row = samples.rows;
	int col = samples.cols;

	for (int i = 0; i < row; i++) {
		for (int k = 0; k < col; k++) {
			fprintf(file, "%f", samples.at<float>(i, k));
			if (k != col)
				fprintf(file, ",");
		}
		fprintf(file, "\n");
	}
	fclose(file);
}
Feature_extraction::CSVData Feature_extraction::readCSV(char* dataset)
{
	char c[10000];
	char* c1;
	char* ends;
	char s[] = " ,\n";//カンマ、スペース、改行
					  /*ファイルオープン*/
	FILE* fp;
	if ((fp = fopen(dataset, "r")) == NULL) {
		printf("Can't open file.\n");
		abort();
	}
	std::vector<std::vector<float>> TS;
	//1行分を文字列としてcに読み込む
	while ((fgets(c, sizeof(c), fp)) != NULL) {
		//c1にカンマ、スペース、改行で区切った文字列を代入
		c1 = strtok(c, s);
		std::vector<float> GTS;
		while (c1 != NULL) {
			//TSにc1の文字列をdoubleに型変換して代入
			GTS.push_back((float)strtod(c1, &ends));
			//NULL文字で終わり
			c1 = strtok(NULL, s);
		}
		TS.push_back(GTS);
	}

	fclose(fp);

	int cols = (int)TS[0].size();
	int rows = (int)TS.size();;

	printf("rows:%d,cols%d\n", rows, cols);
	CSVData Res;
	Res.Labels.create(rows, 1, CV_32SC1);
	Res.samples.create(rows, cols - 1, CV_32F);

	for (int k = 0; k < rows; k++) {
		for (int i = 0; i < cols; i++) {
			if (i == 0) {
				Res.Labels.at<float>(k, 0) = TS[k][0];
			}
			else {
				Res.samples.at<float>(k, i - 1) = TS[k][i];
			}
		}
	}

	return Res;
}