diff --git a/Feature_extraction.cpp b/Feature_extraction.cpp new file mode 100644 index 0000000..90e33bc --- /dev/null +++ b/Feature_extraction.cpp @@ -0,0 +1,893 @@ +#include "Feature_extraction.h" + +//�Z�x���� +cv::Mat Feature_extraction::Density_histgram(Raw_image& ct, array3& labeling, size_t labelcount) +{ + + //�q�X�g�O�����쐬�@-100�` 1200; + auto r_Hist = std::vector(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(i, 0) = (float)r_Hist[i + 1].HistVOXELS; + samples.at(i, 1) = (float)r_Hist[i + 1].Max; + samples.at(i, 2) = (float)r_Hist[i + 1].Min; + samples.at(i, 3) = (float)r_Hist[i + 1].Average; + samples.at(i, 4) = (float)r_Hist[i + 1].Variance; + samples.at(i, 5) = (float)r_Hist[i + 1].Skewness; + samples.at(i, 6) = (float)r_Hist[i + 1].Kurtosis; + } + + return (samples); +} + +//�e�N�X�`������ +array3 Feature_extraction::Histogram_equalization(Raw_image& ct, int len) +{ + //�q�X�g�O�����쐬�@-100�` 1200; + array3 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]; + + //�q�X�g�O�������R�� +#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& image, array3& label, size_t& labelcount,int len) +{ + std::vector s_Sim(labelcount + 1); + + array3 region(image.width(), image.height(), image.depth(), image.reso1(), image.reso1(), image.reso1(), 0); + array3 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(); + //�\���v�f�̐��� + 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; + } + } + + //�������N�s��v�Z + 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; + } + } + + //�����ʒ��o�� + 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]; + } + } + + //�S14�� + //1.�p�x�̓񎟃��[�����g + 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.�R���g���X�g + 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.���U + 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.�������[�����g + 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.���U + 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.���v�G���g���s�[ + 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.�G���g���s�[ + 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.�����U + 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.���G���g���s�[ + 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(i, 0) = (float)s_Sim[i + 1].ASM; + samples.at(i, 1) = (float)s_Sim[i + 1].Contrast; + samples.at(i, 2) = (float)s_Sim[i + 1].Correlation; + samples.at(i, 3) = (float)s_Sim[i + 1].Variance; + samples.at(i, 4) = (float)s_Sim[i + 1].IDM; + samples.at(i, 5) = (float)s_Sim[i + 1].SAverage; + samples.at(i, 6) = (float)s_Sim[i + 1].SVariance; + samples.at(i, 7) = (float)s_Sim[i + 1].SEntropy; + samples.at(i, 8) = (float)s_Sim[i + 1].Entropy; + samples.at(i, 9) = (float)s_Sim[i + 1].DVariance; + samples.at(i, 10) = (float)s_Sim[i + 1].DEntropy; + samples.at(i, 11) = (float)s_Sim[i + 1].IMOC1; + samples.at(i, 12) = (float)s_Sim[i + 1].IMOC2; + samples.at(i, 13) = (float)s_Sim[i + 1].MCC; + samples.at(i, 14) = (float)s_Sim[i + 1].DASM; + samples.at(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; +} + +//���z���� +cv::Mat Feature_extraction::Laplace_feature(Raw_image& ct, array3& label, size_t& labelcount) +{ + mist::array3 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 CTb, laplace; + mist::gaussian_filter(CT, CTb); + mist::laplacian_filter(CTb, laplace); + auto r_Hist = std::vector(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(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("�T���v����������������܂���\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(i, k) = (float)firstMat.at(i, k); + else + samples.at(i, k) = (float)secondMat.at(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(i, k) = (float)firstMat.at(i, k); + else + samples.at(i, k) = (float)secondMat.at(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(i, k) = (float)Label; + else + samples.at(i, k) = (float)input.at(i, k - 1); + } + } + + return samples; +} + +//�X�P�[�����O +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(0); + for (int j = 0; j < feature.cols; j++) { + p[j] = (p[j] - mean[j]) / std[j]; + } + } +} + +//�X�N���[�j���O +cv::Mat Feature_extraction::Screening(cv::Mat& samples) +{ + int rows = samples.rows; + int count = 0; + std::vector s; + for (int i = 0; i < rows; i++) { + if (samples.at(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(i, j) = s[i].at(0, j); + } + } + + return Return; +} + +//���ׂĂ̓�������‚̃}�b�g�ɕۑ�����֐� +cv::Mat Feature_extraction::All_feature_exstract(Raw_image ct, array3& 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; +} + +//�w�K�f�[�^�̐�������ƌ댟�o�����̐��̒��� +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(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; // �񌈒�I�ȗ���������𐶐� + std::mt19937 mt(rnd()); // �����Z���k�E�c�C�X�^��32�r�b�g�ŁA�����͏����V�[�h�l + std::uniform_int_distribution<> rand100(0, Rate - 1); + + for (int i = 0; i < samples.rows; i++) { + if (samples.at(i, 0) == 1) { + if (TCount < Label1) { + for (int k = 0; k < samples.cols; k++) { + TrueMat.at(TCount, k) = samples.at(i, k); + } + TCount++; + } + + } + else { + if (Falserate == -1) { + if (FCount < FalseMat.rows) { + for (int k = 0; k < samples.cols; k++) { + FalseMat.at(FCount, k) = samples.at(i, k); + } + FCount++; + } + + } + else { + if (rand100(mt) < 5) { + if (FCount < FalseMat.rows) { + for (int k = 0; k < samples.cols; k++) { + FalseMat.at(FCount, k) = samples.at(i, k); + } + FCount++; + } + } + + } + + } + } + printf("�I�茋��:1:%d,0:%d\n", TCount, FCount); + + cv::Mat Return = Addsample(TrueMat, FalseMat); + + std::random_device rnda; // �񌈒�I�ȗ���������𐶐� + std::mt19937 mta(rnda()); // �����Z���k�E�c�C�X�^��32�r�b�g�ŁA�����͏����V�[�h�l + std::uniform_int_distribution<> randes(0, Return.rows - 1); // [0, 99] �͈͂̈�l���� + + if (shafful) { + printf("�V���b�t���J�n\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(i, k) = tmp.at(0, k); + Return.at(randres, k) = tmp1.at(0, k); + } + } + } + + return Return; +} + +//�w�K�������f����p���Ă̕��� +array3 Feature_extraction::Predict(cv::Ptr r_tree, array3& 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 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�̓��o�� +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(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";//�J���}�A�X�y�[�X�A���s + /*�t�@�C���I�[�v��*/ + FILE* fp; + if ((fp = fopen(dataset, "r")) == NULL) { + printf("Can't open file.\n"); + abort(); + } + std::vector> TS; + //1�s���𕶎���Ƃ���c�ɓǂݍ��� + while ((fgets(c, sizeof(c), fp)) != NULL) { + //c1�ɃJ���}�A�X�y�[�X�A���s�ŋ�؂������������ + c1 = strtok(c, s); + std::vector GTS; + while (c1 != NULL) { + //TS��c1�̕������double�Ɍ^�ϊ����đ�� + GTS.push_back((float)strtod(c1, &ends)); + //NULL�����ŏI��� + 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(k, 0) = TS[k][0]; + } + else { + Res.samples.at(k, i - 1) = TS[k][i]; + } + } + } + + return Res; +} diff --git a/Feature_extraction.h b/Feature_extraction.h new file mode 100644 index 0000000..53b82bb --- /dev/null +++ b/Feature_extraction.h @@ -0,0 +1,134 @@ +#pragma once +#include "raw_image_class.h" +#include "extravasation_detect.h" +#include +//using namespace cv; + +//���ފ��������o�Ɋւ���N���X +//���̑��ɂ�Mat�̘A���Ȃǂ̑�����܂܂�� +static class Feature_extraction +{ +public : + //�Z�x�����ʂ̍\���� + struct DensityData { + int* Hist; + double* NormalizeHist; + double* CPD; + const double t = 0.03; + + //�q�X�g�O�����͈͂̉�f�� + int HistVOXELS = 0; + //�ő�l + short Max = -999; + //�ŏ��l + short Min = -999; + //���� + double Average = 0; + //�c�x + double Skewness = 0; + //���U + double Variance = 0; + //��x + double Kurtosis = 0; + }; + //�Z�x�������o + static cv::Mat Density_histgram(Raw_image& ct, array3& labeling, size_t labelcount); + + //GLCM�̍\���� + struct GLCMData { + float** simMatrinx; + float** Q; + float Mf = -999; + float Mn = -999; + float* Py; + float* Px; + float* Pxpy; + float* Pxsy; + float SAvey = 0; + float Ax = 0; + float Ay = 0; + float Sx = 0; + float Sy = 0; + float HX = 0; + float HY = 0; + float HXY1 = 0; + float HXY2 = 0; + //�ŕp�l + float mode = 0; + + //�͈͂̉�f�� + float MatrixCount = 0; + + ////�擾��������� + //�p�x�̓񎟃��[�����g + float ASM = 0; + float DASM = 0; + + //�R���g���X�g + float Contrast = 0; + //���� + float Correlation = 0; + //���U + float Variance = 0; + //�t�����[�����g + float IDM = 0; + //���� + float SAverage = 0; + //���U + float SVariance = 0; + //�G���g���s�[ + float SEntropy = 0; + //�G���g���s�[ + float Entropy = 0; + //�����U + float DVariance = 0; + //���G���g���s�[ + float DEntropy = 0; + //���ւ̏��� + float IMOC1 = 0; + //���ւ̏��� + float IMOC2 = 0; + //�ő告�֌W�� + float MCC = 0; + }; + //�q�X�g�O�����ϓ��� + static array3 Histogram_equalization(Raw_image& ct, int len); + //GLCM + static cv::Mat GLCM(array3& image, array3& label, size_t& labelcount,int len); + + //���v���V�A������ + static cv::Mat Laplace_feature(Raw_image& ct, array3& label, size_t& labelcount); + + //Mat���� + static cv::Mat MatLink(cv::Mat& firstMat, cv::Mat& secondMat); + static cv::Mat Addsample(cv::Mat& firstMat, cv::Mat& secondMat); + static cv::Mat LabelMat(cv::Mat& input, int Label); + + //�X�P�[�����O + static void Scaling(cv::Mat& feature); + + //�X�N���[�j���O + static cv::Mat Screening(cv::Mat& samples); + + //�S���̓����̒��o + static cv::Mat All_feature_exstract(Raw_image ct, array3& label, size_t& labelcount, int len); + + //�f�[�^���̍팸 + static cv::Mat Feature_Cut(cv::Mat& samples, bool shafful, float Falserate); + //���� + static array3 Predict(cv::Ptr r_tree, array3& label, size_t& labelcount, cv::Mat samples); + + + + //csv�t�@�C���̎�舵�� + static void writeCSV(char* FileName, cv::Mat samples); + + //CSV�o�͂����f�[�^�̓ǂݍ��� + struct CSVData { + cv::Mat Labels; + cv::Mat samples; + }; + + static CSVData readCSV(char* dataset); +}; + diff --git a/Pass_management.h b/Pass_management.h new file mode 100644 index 0000000..5b3db07 --- /dev/null +++ b/Pass_management.h @@ -0,0 +1,59 @@ +#pragma once +////�p�X�Ɋւ���t�@�C��(�������‹��ɉ����ĕύX) +//���[�g�f�B���N�g�� +#define DIR "D:\\Ktest" +//#define DIR "D:\\YL\\Data\\BleedingDataNew" + +/*���̓t�@�C���̃p�X*/ +//�𑜓x +#define SIZEFILE "\\%d\\iso%d" +//CT +#define CTI "\\%d\\delay_isorotopic%d" +//�����摜 +#define MAK "\\%d\\maskiso%d" + +/*���ʂ̏o��*/ +#define OUTNAME "0228OLD" +//���̒��o��@�̕ύX +enum BONEMODE { OLD, WATER, ADAPTIVE }; +#define BE ADAPTIVE + +#define CASENUM "\\%d\\" +//���̏o�͐� +#define THREBON "BoneIsoNew%d" +#define NAME_DOUBLE "Double" +//���� +//#define RES "ResultIsoNew%d" +#define RES "OldThreshResult%d" +//�]���摜(���ޑO�ƌ�) +#define EVA "EvaIsoNew%d" +#define PREEVA "PreEvaIsoNew%d" + +/*�@�B�w�K���f��*/ +//�댟�o���ފ�(Random Forest:OpenCV) +#define MODEL "\\rtree_model_new" +//�X���C�X�ʒu���胂�f�� +#define CLASS "\\ClassModel.pb" +#define SHOULDER "\\ShoulderModel.pb" +#define WAIST "\\LowModel.pb" + +//�w�K�̃f�[�^�Z�b�g +#define DST "\\ClassiferDataIso" +#define RCSV "\\DataforROCEvalu" + +//�S�̌��ʂ�CSV�o�� +#define DOUBLEOUT "\\���ޑO����" +#define DOUBLEPRE "\\���ތ㌋��" + +#define FROC "\\FROC" +#define ROC "\\ROC" +#define REGIONDATA "\\REGIONDATA" + +//�g���q +#define RAW ".raw" +#define GZ ".gz" +#define CSV ".csv" +#define XML ".xml" +#define TXT ".txt" +#define FLO "%f" +//�^��?�ʔ�?�F31�C \ No newline at end of file diff --git a/extravasation_detect.cpp b/extravasation_detect.cpp new file mode 100644 index 0000000..3ceb8a3 --- /dev/null +++ b/extravasation_detect.cpp @@ -0,0 +1,1256 @@ +#include "extravasation_detect.h" +//using namespace cv; +#define MAXW 400 +#define MINW 2 +#define MINA 600 +#define THRESHOLDCENTER 425 +#define FAI -999 + +//�����o�֘A +array3 Extravasation_detect::Artificial_Segmentation(Raw_image& ct) +{ + auto Artficial = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + auto Artf = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + auto result = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + std::vector> Points; + + for (int z = 0; z < ct.data.Depth; z++) { + for (int y = 0; y < ct.data.Height; y++) + for (int x = 0; x < ct.data.Width; x++) { + int point = z * ct.data.IMAGESIZE + y * ct.data.Height + x; + + if (ct.image[point] > 600) { + Artficial[point] = 180; + } + + if (ct.image[point] > 3500) { + Artf[point] = 255; + Points.push_back(mist::vector3(x, y, z)); + + } + } + } + + mist::region_growing_utility::equal equal(180); + mist::region_growing_utility::sphere component(1); + array3 Artregion; + + auto count = mist::region_growing(Artficial, Artregion, Points, 255, component, equal); + + for (int i = 0; i < ct.data.VOXELS; i++) { + if (Artregion[i] != 0) + result[i] = 180; + } + return result; +} +int Extravasation_detect::Adaptive_Threshold(Raw_image& ct) +{ + + array3 art = Artificial_Segmentation(ct); + + short* tempCT = new short[ct.data.VOXELS]; + for (int i = 0; i < ct.data.VOXELS; i++) { + if (art[i] != 0) { + tempCT[i] = 0; + } + else { + tempCT[i] = ct.image[i]; + } + } + + array3 Upper = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, + ct.data.reso, ct.data.reso, + &tempCT[0], sizeof(short) * ct.data.IMAGESIZE * ct.data.Depth); + + int hist[500]; + for (int i = 0; i < 500; i++) + hist[i] = 0; + + unsigned long All = 0; + + for (int i = 0; i < Upper.size(); i++) { + if (Upper[i] > 50) { + All++; + if (Upper[i] >= 500) { + hist[499]++; + } + else { + hist[Upper[i]]++; + } + } + } + + int thre = 499; + + for (int i = thre; i > 50; i--) { + double bone = 0; + for (int j = 499; j >= i; j--) { + bone += hist[j]; + } + double p = bone * 100 / (double)All; + if (p > 8.5) { + thre = i; + break; + } + } + + delete[] tempCT; + return thre; +} +array3 Extravasation_detect::Watershed(short* ct, array3 prebone, int thre) +{ + array3 Return(prebone.width(), prebone.height(), prebone.depth(), prebone.reso1(), prebone.reso2(), prebone.reso3(), 0); + + + for (int slice = 0; slice < prebone.depth(); slice++) { + cv::Mat image(512, 512, CV_16S, (short*)(ct + slice * 512 * 512)); + cv::Mat image_w; + cv::Mat image_disp; + //printf("%d\n",slice); + cv::Mat front(512, 512, CV_8U, (uchar*)&prebone[slice * 512 * 512]); + cv::Mat backmarker(512, 512, CV_8U, cv::Scalar::all(255)); + + cv::Mat dfront; + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + morphologyEx(front, dfront, cv::MORPH_DILATE, elemental, cv::Point(-1, -1), 25); + //imshow("dfront", dfront); + backmarker = backmarker - dfront; + cv::Mat Bone = front.clone(); + array2 thinf(512, 512, prebone.reso1(), prebone.reso2(), (uchar*)front.ptr(0), sizeof(uchar) * 512 * 512); + mist::hilditch::thinning(thinf, thinf); + uchar* pf = front.ptr(0); + for (int i = 0; i < 512 * 512; i++, pf++) { + if (thinf[i] != 0) { + *pf = 255; + } + else + *pf = 0; + } + + cv::Mat fL = front.clone(); + + cv::Mat disp(512, 512, CV_8UC3, cv::Scalar::all(0)); + disp.setTo(cv::Scalar(255, 255, 255), backmarker); + + cv::Mat Marker; + int Label = cv::connectedComponents(fL, Marker, 8, CV_32S); + + float grad = 255.0 / 300.0; + float dispg = 255.0 / 400; + float rg = 255.0 / 600; + cv::Mat aho; + image.convertTo(image_w, CV_8U, grad, 127.5 - ((double)thre * grad)); + image.convertTo(image_disp, CV_8U, dispg, 127.5 - (155.0 * dispg)); + image.convertTo(aho, CV_8U, rg, 127.5 - (900.0 * rg)); + cv::GaussianBlur(image_w.clone(), image_w, cv::Size(3, 3), 1.2); + + short* imp = image.ptr(0); + uchar* imwp = image_w.ptr(0); + for (int i = 0; i < 512 * 512; i++) { + if (imp[i] <= -2040) { + imwp[i] = 255; + } + } + + cv::Mat Laplacian; + cv::Laplacian(image_w, Laplacian, CV_8U, 3); + cv::cvtColor(Laplacian.clone(), Laplacian, cv::COLOR_GRAY2BGR); + cv::cvtColor(image_w.clone(), image_w, cv::COLOR_GRAY2BGR); + + backmarker.convertTo(backmarker, CV_32S); + Marker = Marker + backmarker; + watershed(Laplacian, Marker); + + std::vector colors; + colors.push_back(cv::Vec3b(0, 0, 0)); + srand((unsigned int)time(NULL)); + for (int i = 1; i <= Label; i++) { + colors.push_back(cv::Vec3b(rand() % 255, rand() % 255, rand() % 255)); + } + + cv::Mat DstU(Marker.size(), CV_8U); + cv::Mat Dst(Marker.size(), CV_8UC3); + for (int i = 0; i < Dst.rows; ++i) { + int* lb = Marker.ptr(i); + cv::Vec3b* pix = Dst.ptr(i); + uchar* rix = DstU.ptr(i); + for (int j = 0; j < Dst.cols; ++j) { + if (lb[j] == 255) { + rix[j] = 0; + pix[j] = cv::Vec3b(255, 255, 255); + } + else { + rix[j] = 255; + pix[j] = colors[lb[j]]; + } + } + } + + uchar* Up = DstU.ptr(0); + + for (int i = 0; i < 512 * 512; i++, Up++) { + if (i % 512 == 0 || i < 512 || i > 512 * 511 || i % 512 == 511) + Return[slice * 512 * 512 + i] = 0; + else + Return[slice * 512 * 512 + i] = *Up; + } + + } + + return Return; +} +array3Extravasation_detect::Old_BoneExtractMethod(Raw_image& ct, BODYTYPE type) +{ + short* Bone = new short[ct.data.VOXELS]; + array3 gets(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + int exstaract, region; + if (type == LOWER) { + exstaract = 350; region = 340; + } + else if (type == HEAD) { + exstaract = 350; region = 340; + } + else { + exstaract = 350; region = 340; + } + + + for (int i = 0; i < ct.data.VOXELS; i++) { + if (ct.image[i] > exstaract) + Bone[i] = ct.image[i]; + else + Bone[i] = 0; + + if (ct.image[i] > region) + gets[i] = 180; + } + + array3 art = Artificial_Segmentation(ct); + + if (type != HEAD) { + for (int i = 0; i < ct.data.VOXELS; i++) { + if (art[i] != 0) { + Bone[i] = 0; + } + } + } + + auto BoneOnly = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + +#pragma omp parallel for + for (int i = 0; i < ct.data.VOXELS; i++) { + if (Bone[i] > 0) + BoneOnly[i] = 255; + else + BoneOnly[i] = 0; + } + + return BoneOnly; + + auto BoneLabel = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + size_t BoneCount = mist::labeling6(BoneOnly, BoneLabel); + + //���̑̐� + std::vector CTV = std::vector(BoneCount + 1); + + int* high = new int[BoneCount + 1]; + int* maxthre = new int[BoneCount + 1]; + + for (int i = 0; i < BoneCount + 1; i++) { + high[i] = 0; + maxthre[i] = 0; + } + + for (int slice = 0; slice < ct.data.Depth; slice++) + for (int j = 50; j < 450; j++) + for (int x = 0; x < ct.data.Width; x++) { + + //���o���ʂ��烉�x�����O���ꂽ�̈�𔭌� + if (BoneLabel[slice * ct.data.IMAGESIZE + j * ct.data.Width + x] != 0) { + int pointer = slice * ct.data.IMAGESIZE + j * ct.data.Width + x; + //�̐ς̏������ + CTV[BoneLabel[pointer]].Area++; + if (CTV[BoneLabel[pointer]].Area == 1) { + CTV[BoneLabel[pointer]].x = x; + CTV[BoneLabel[pointer]].y = j; + CTV[BoneLabel[pointer]].z = slice; + } + + CTV[BoneLabel[pointer]].CTValue += ct.image[pointer]; + if (ct.image[pointer] > 800) + high[BoneLabel[pointer]] ++; + if (ct.image[pointer] > 2200) + maxthre[BoneLabel[pointer]]++; + } + + } + + //������ + bool* BoneCheck = new bool[BoneCount + 1]; + std::vector> Points; + int bonedetect = 0; + + for (int i = 1; i < BoneCount + 1; i++) { + BoneCheck[i] = true; + + //������� + if (CTV[i].Area < 10000) { + //if (high[i] < 1500 || maxthre[i] > 50) { + //��������O�̃��x����false�Ƃ��Ă��� + BoneCheck[i] = false; + bonedetect++; + //} + } + else { + Points.push_back(mist::vector3(CTV[i].x, CTV[i].y, CTV[i].z)); + bonedetect++; + } + } + + auto BoneDelete = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + auto BoneOut = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + +#pragma omp parallel for + for (int i = 0; i < BoneLabel.size(); i++) { + if (BoneLabel[i] != 0) { + //bonecheck��false�̗̈�͏��� + if (!BoneCheck[BoneLabel[i]]) { + BoneOut[i] = 0; + BoneDelete[i] = 1; + } + else { + BoneOut[i] = 255; + gets[i] = 180; + } + } + + } + + if (type == LOWER) { + + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat original(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + + uchar* pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + if (BoneOut[slice * ct.data.IMAGESIZE + i] != 0) + *pm = BoneOut[slice * ct.data.IMAGESIZE + i]; + } + + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 6); + morphologyEx(original.clone(), original, cv::MORPH_DILATE, elemental, cv::Point(-1, -1), 2); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 2); + + //�O�֊s���o + std::vector> contours; + cv::Mat maskb(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + cv::findContours(original, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE); + for (int i = 0; i < contours.size(); i++) + cv::drawContours(maskb, contours, i, cv::Scalar(255), cv::FILLED); + + original = maskb; + + pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + BoneOut[slice * ct.data.IMAGESIZE + i] = *pm; + + if (ct.image[slice * ct.data.IMAGESIZE + i] > 600 || art[slice * ct.data.IMAGESIZE + i] != 0) { + BoneOut[slice * ct.data.IMAGESIZE + i] = 255; + } + } + + } + + mist::closing(BoneOut, 2, 0); + } + if (type == HEAD) { + mist::region_growing_utility::equal get(180); + mist::region_growing_utility::sphere component(1); + array3 Boneregion; + mist::region_growing(gets, Boneregion, Points, 140, component, get); + for (int i = 0; i < ct.data.VOXELS; i++) { + if (Boneregion[i] == 140) + BoneOut[i] = 255; + } + + + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat original(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + + uchar* pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + if (BoneOut[slice * ct.data.IMAGESIZE + i] != 0) + *pm = BoneOut[slice * ct.data.IMAGESIZE + i]; + } + + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 5); + morphologyEx(original.clone(), original, cv::MORPH_DILATE, elemental, cv::Point(-1, -1), 2); + + pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + BoneOut[slice * ct.data.IMAGESIZE + i] = *pm; + + if (ct.image[slice * ct.data.IMAGESIZE + i] > 600) { + BoneOut[slice * ct.data.IMAGESIZE + i] = 255; + } + } + + } + //mist::dilation(BoneOut, 2, 0); + } + if (type == UPPER) { + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat original(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + + uchar* pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + if (BoneOut[slice * ct.data.IMAGESIZE + i] != 0) + *pm = BoneOut[slice * ct.data.IMAGESIZE + i]; + } + + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 5); + morphologyEx(original.clone(), original, cv::MORPH_DILATE, elemental, cv::Point(-1, -1), 2); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 3); + + + //�O�֊s���o + std::vector> contours; + cv::Mat maskb(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + cv::findContours(original, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE); + for (int i = 0; i < contours.size(); i++) + cv::drawContours(maskb, contours, i, cv::Scalar(255), cv::FILLED); + + original = maskb; + + pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + BoneOut[slice * ct.data.IMAGESIZE + i] = *pm; + + //if (ct.image[slice * ct.data.IMAGESIZE + i] > 600 || art[slice * ct.data.IMAGESIZE + i] != 0) { + // BoneOut[slice * ct.data.IMAGESIZE + i] = 255; + //} + } + + } + mist::closing(BoneOut, 3, 0); + //mist::closing(BoneOut, 3, 0); + } + + delete[] maxthre; + delete[] high; + delete[] Bone; + delete[] BoneCheck; + return(BoneOut); +} +array3 Extravasation_detect::Bone_Extract(Raw_image& ct, BODYTYPE type, int threshold) +{ + + short* Bone = new short[ct.data.VOXELS]; + array3 gets(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + int exstaract = threshold; + + for (int i = 0; i < ct.data.VOXELS; i++) { + if (ct.image[i] > exstaract) + Bone[i] = ct.image[i]; + else + Bone[i] = 0; + + } + + array3 art = Artificial_Segmentation(ct); + + if (type != HEAD) { + for (int i = 0; i < ct.data.VOXELS; i++) { + if (art[i] != 0) { + Bone[i] = 0; + } + } + } + + auto BoneOnly = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + +#pragma omp parallel for + for (int i = 0; i < ct.data.VOXELS; i++) { + if (Bone[i] > 0) + BoneOnly[i] = 255; + else + BoneOnly[i] = 0; + } + + mist::closing(BoneOnly, 1); + + auto BoneLabel = mist::array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + size_t BoneCount = mist::labeling6(BoneOnly, BoneLabel); + + //���̑̐� + std::vector CTV = std::vector(BoneCount + 1); + + + for (int slice = 0; slice < ct.data.Depth; slice++) + for (int j = 0; j < ct.data.Height; j++) + for (int x = 0; x < ct.data.Width; x++) { + + //���o���ʂ��烉�x�����O���ꂽ�̈�𔭌� + if (BoneLabel[slice * ct.data.IMAGESIZE + j * ct.data.Width + x] != 0) { + int pointer = slice * ct.data.IMAGESIZE + j * ct.data.Width + x; + //�̐ς̏������ + CTV[BoneLabel[pointer]].Area++; + if (CTV[BoneLabel[pointer]].Area == 1) { + CTV[BoneLabel[pointer]].x = x; + CTV[BoneLabel[pointer]].y = j; + CTV[BoneLabel[pointer]].z = slice; + } + + } + + } + + //������ + bool* BoneCheck = new bool[BoneCount + 1]; + std::vector> Points; + int bonedetect = 0; + + for (int i = 1; i < BoneCount + 1; i++) { + BoneCheck[i] = true; + + //������� + if (CTV[i].Area < 7000) { + BoneCheck[i] = false; + bonedetect++; + } + else { + Points.push_back(mist::vector3(CTV[i].x, CTV[i].y, CTV[i].z)); + bonedetect++; + } + } + + auto BoneDelete = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + auto BoneOut = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + +#pragma omp parallel for + for (int i = 0; i < BoneLabel.size(); i++) { + if (BoneLabel[i] != 0) { + //bonecheck��false�̗̈�͏��� + if (!BoneCheck[BoneLabel[i]]) { + BoneOut[i] = 0; + BoneDelete[i] = 1; + } + else { + BoneOut[i] = 255; + gets[i] = 180; + } + } + + } + + if (type == LOWER) { + + BoneOut = Watershed(ct.image,BoneOut, threshold); + + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat original(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + + uchar* pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + if (BoneOut[slice * ct.data.IMAGESIZE + i] != 0) + *pm = BoneOut[slice * ct.data.IMAGESIZE + i]; + } + + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + //morphologyEx(original.clone(), original, MORPH_CLOSE, elemental, cv::Point(-1, -1), 6); + //morphologyEx(original.clone(), original, MORPH_DILATE, elemental, cv::Point(-1, -1), 2); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 2); + + //�O�֊s���o + std::vector> contours; + cv::Mat maskb(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + cv::findContours(original, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE); + for (int i = 0; i < contours.size(); i++) + cv::drawContours(maskb, contours, i, cv::Scalar(255), cv::FILLED); + + original = maskb; + + pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + BoneOut[slice * ct.data.IMAGESIZE + i] = *pm; + } + + } + + mist::closing(BoneOut, 1, 0); + } + if (type == HEAD) { + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat original(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + + uchar* pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + if (BoneOut[slice * ct.data.IMAGESIZE + i] != 0) + *pm = BoneOut[slice * ct.data.IMAGESIZE + i]; + } + + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 2); + morphologyEx(original.clone(), original, cv::MORPH_DILATE, elemental, cv::Point(-1, -1), 2); + + pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + BoneOut[slice * ct.data.IMAGESIZE + i] = *pm; + + } + + } + //mist::dilation(BoneOut, 2, 0); + } + if (type == UPPER) { + BoneOut = Watershed(ct.image, BoneOut, threshold); + + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat original(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + + uchar* pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + if (BoneOut[slice * ct.data.IMAGESIZE + i] != 0) + *pm = BoneOut[slice * ct.data.IMAGESIZE + i]; + } + + cv::Mat elemental(3, 3, CV_8U, cv::Scalar::all(255)); + //morphologyEx(original.clone(), original, MORPH_CLOSE, elemental, cv::Point(-1, -1), 5); + //morphologyEx(original.clone(), original, MORPH_DILATE, elemental, cv::Point(-1, -1), 1); + morphologyEx(original.clone(), original, cv::MORPH_CLOSE, elemental, cv::Point(-1, -1), 2); + + + //�O�֊s���o + std::vector> contours; + cv::Mat maskb(ct.data.Height, ct.data.Width, CV_8U, cv::Scalar(0)); + cv::findContours(original, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE); + for (int i = 0; i < contours.size(); i++) + cv::drawContours(maskb, contours, i, cv::Scalar(255), cv::FILLED); + + original = maskb; + + pm = original.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, pm++) { + BoneOut[slice * ct.data.IMAGESIZE + i] = *pm; + + //if (ct[slice * ct.data.IMAGESIZE + i] > 600 || art[slice * ct.data.IMAGESIZE + i] != 0) { + // BoneOut[slice * ct.data.IMAGESIZE + i] = 255; + //} + } + + } + mist::closing(BoneOut, 3, 0); + + } + + delete[] Bone; + delete[] BoneCheck; + return(BoneOut); +} + +//���̈撊�o�֘A +array3 Extravasation_detect::double_thresholding(Raw_image& ct, array3& bone, BODYTYPE type) +{ + array3 Return(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + + short* DataL = new short[ct.data.VOXELS]; + short* DataH = new short[ct.data.VOXELS]; + Parameter PL, PH; + if (type == LOWER) { + PL = Parameter{ 1000,160,100,15,90,3,2.0,10000 }; + PH = Parameter{ 1000,180,115,30,114,3,5,75000 }; + } + else if (type == HEAD) { + PL = Parameter{ 1000,155,105,20,102,3,2.0,10000 }; + PH = Parameter{ 1000,265,180,30,162,3,5,75000 }; + } + else { + PL = Parameter{ 1000,185,95,20,86,3,2.0,10000 }; + PH = Parameter{ 1000,300,140,30,136,3,5,75000 }; + } + + for (int i = 0; i < ct.data.VOXELS; i++) { + DataL[i] = DataH[i] = ct.image[i]; + } + + //CT�l���� + for (int i = 0; i < ct.data.VOXELS; i++) { + if (DataL[i] > PL.HighCut) + DataL[i] = 0; + + + if (bone[i] != 0) + DataL[i] = 0; + + if (DataL[i] < -500) + DataL[i] = -2000; + + if (DataH[i] > PH.HighCut) + DataH[i] = 0; + + + if (bone[i] != 0) + DataH[i] = 0; + + if (DataH[i] < -500) + DataH[i] = -2000; + } + + //���̈撊�o + array3 VenoOnlyL = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso); + array3 VenoOnlyH = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso); + array3 LV = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso); + array3 HV = array3(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso); + + + //�Ö������i�[�p�̃R���e�i + //�Ö��̃J�b�g�I�t + for (int slice = 0; slice < ct.data.Depth; slice++) { + cv::Mat blurL(ct.data.Height, ct.data.Width, CV_16S, (short*)(DataL + ct.data.Width * ct.data.Height * slice)); + cv::Mat blurH(ct.data.Height, ct.data.Width, CV_16S, (short*)(DataH + ct.data.Width * ct.data.Height * slice)); + cv::Mat getblurL, getblurH; + //�K�E�V�A���t�B���^ + GaussianBlur(blurL, getblurL, cv::Size(PL.ExtractKarnel, PL.ExtractKarnel), (float)PL.sigma); + GaussianBlur(blurH, getblurH, cv::Size(PH.ExtractKarnel, PH.ExtractKarnel), (float)PH.sigma); + short* p = getblurL.ptr(0); + short* p2 = getblurH.ptr(0); + for (int i = 0; i < ct.data.IMAGESIZE; i++, p++, p2++) { + + //���W + int point = i + slice * ct.data.IMAGESIZE; + if (*p < PL.CutVeno || bone[point] != 0) + LV[point] = VenoOnlyL[point] = 0; + else + LV[point] = VenoOnlyL[point] = 255; + + if (*p2 < PH.CutVeno || bone[point] != 0) + HV[point] = VenoOnlyH[point] = 0; + else + HV[point] = VenoOnlyH[point] = 255; + } + } + + //�댟�o���� + //���x�����O + array3 VenoLabelL; + array3 VenoLabelH; + size_t VenoCountL = mist::labeling6(VenoOnlyL, VenoLabelL); + size_t VenoCountH = labeling6(VenoOnlyH, VenoLabelH); + //�����ݒ� + auto VenoPointL = std::vector(VenoCountL + 1); + auto VenoPointH = std::vector(VenoCountH + 1); + + VenoPointL[0].Delete = VenoPointH[0].Delete = false; + VenoPointL[0].Detectcheck = VenoPointH[0].Detectcheck = false; + + for (int i = 0; i < ct.data.Depth; i++) + for (int j = 20; j < 500; j++) + for (int x = 0; x < ct.data.Height; x++) { + + //���W + int point = i * ct.data.IMAGESIZE + j * ct.data.Width + x; + if (VenoLabelL[point] != 0) { + VenoPointL[VenoLabelL[point]].Area++; + VenoPointL[VenoLabelL[point]].CTValue += DataL[point]; + VenoPointL[VenoLabelL[point]].xCount += x; + VenoPointL[VenoLabelL[point]].yCount += j; + VenoPointL[VenoLabelL[point]].zCount += i; + if (VenoPointL[VenoLabelL[point]].minz == 0 || VenoPointL[VenoLabelL[point]].minz > i) + VenoPointL[VenoLabelL[point]].minz = i; + if (VenoPointL[VenoLabelL[point]].maxz == 0 || VenoPointL[VenoLabelL[point]].maxz < i) + VenoPointL[VenoLabelL[point]].maxz = i; + } + if (VenoLabelH[point] != 0) { + VenoPointH[VenoLabelH[point]].Area++; + VenoPointH[VenoLabelH[point]].CTValue += DataH[point]; + VenoPointH[VenoLabelH[point]].xCount += x; + VenoPointH[VenoLabelH[point]].yCount += j; + VenoPointH[VenoLabelH[point]].zCount += i; + if (VenoPointH[VenoLabelH[point]].minz == 0 || VenoPointH[VenoLabelH[point]].minz > i) + VenoPointH[VenoLabelH[point]].minz = i; + if (VenoPointH[VenoLabelH[point]].maxz == 0 || VenoPointH[VenoLabelH[point]].maxz < i) + VenoPointH[VenoLabelH[point]].maxz = i; + } + } + + //���� + for (int i = 1; i < VenoCountL + 1; i++) { + //������� + if (VenoPointL[i].Area < PL.MINV || VenoPointL[i].Area > PL.MaxVolume || VenoPointL[i].CTValue / VenoPointL[i].Area > PL.MaxCTValue || VenoPointL[i].CTValue / VenoPointL[i].Area < PL.MinCTValue) { + VenoPointL[i].Delete = false; + VenoPointL[i].Detectcheck = false; + } + } + + for (int i = 1; i < VenoCountH + 1; i++) { + //������� + if (VenoPointH[i].Area < PH.MINV || VenoPointH[i].Area > PH.MaxVolume || VenoPointH[i].CTValue / VenoPointH[i].Area > PH.MaxCTValue || VenoPointH[i].CTValue / VenoPointH[i].Area < PH.MinCTValue) { + VenoPointH[i].Delete = false; + VenoPointH[i].Detectcheck = false; + } + } + + for (int i = 1; i < VenoCountL + 1; i++) { + if (VenoPointL[i].Delete) { + VenoPointL[i].x = (int)(VenoPointL[i].xCount / VenoPointL[i].Area); + VenoPointL[i].y = (int)(VenoPointL[i].yCount / VenoPointL[i].Area); + VenoPointL[i].z = (int)(VenoPointL[i].zCount / VenoPointL[i].Area); + + int width = VenoPointL[i].maxz - VenoPointL[i].minz; + + if (width == 0) { + VenoPointL[i].Delete = false; + VenoPointL[i].Detectcheck = false; + } + else { + int AveArea = VenoPointL[i].Area / width; + + if (VenoPointL[i].y > THRESHOLDCENTER || width > MAXW || width < MINW || AveArea > MINA) { + VenoPointL[i].Delete = false; + VenoPointL[i].Detectcheck = false; + } + + if (width > 120 && AveArea < 80) { + VenoPointL[i].Delete = false; + VenoPointL[i].Detectcheck = false; + } + } + + } + } + + for (int i = 1; i < VenoCountH + 1; i++) { + if (VenoPointH[i].Delete) { + VenoPointH[i].x = (int)(VenoPointH[i].xCount / VenoPointH[i].Area); + VenoPointH[i].y = (int)(VenoPointH[i].yCount / VenoPointH[i].Area); + VenoPointH[i].z = (int)(VenoPointH[i].zCount / VenoPointH[i].Area); + + int width = VenoPointH[i].maxz - VenoPointH[i].minz; + + if (width == 0) { + VenoPointH[i].Delete = false; + VenoPointH[i].Detectcheck = false; + } + else { + int AveArea = VenoPointH[i].Area / width; + + if (VenoPointH[i].y > THRESHOLDCENTER || width > MAXW || width < MINW || AveArea > MINA) { + VenoPointH[i].Delete = false; + VenoPointH[i].Detectcheck = false; + } + + if (width > 120 && AveArea < 80) { + VenoPointH[i].Delete = false; + VenoPointH[i].Detectcheck = false; + } + } + + } + } + +#pragma omp parallel for + for (int i = 0; i < VenoLabelL.size(); i++) { + if (VenoLabelL[i] != 0 || VenoLabelH[i] != 0) { + if (!VenoPointL[VenoLabelL[i]].Delete && !VenoPointH[VenoLabelH[i]].Delete) { + Return[i] = 0; + } + else { + Return[i] = 255; + } + } + else { + Return[i] = 0; + } + } + + delete[] DataL; + delete[] DataH; + return (Return); +} + +template +inline array3 Extravasation_detect::unshape(array3& input) +{ + array3 CTb, laplace; + mist::gaussian_filter(input, CTb); + array3 sub(input.width(), input.height(), input.depth(), 0.5, 0.5, 0.5, 0); + for (int i = 0; i < input.size(); i++) { + sub[i] = 2 * input[i] - CTb[i]; + } + + return sub; +} +std::vector> Extravasation_detect::reduction_point(array3 ct, std::vector>& point, float rate) +{ + { + float max = (float)point.size(); + int hist[600]; + std::vector> result; + for (int i = 0; i < 600; i++) hist[i] = 0; + + for (int i = 0; i < max; i++) { + int p = ct.width() * ct.height() * point[i].z + ct.width() * point[i].y + point[i].x; + if (ct[p] >= 600) hist[599]++; + else if (ct[p] < 0) hist[0] ++; + else hist[ct[p]]++; + } + + float count = 0; + int thre; + for (int i = 599; i > 0; i--) { + count += hist[i]; + + if (count / max > rate) { + thre = i; + break; + } + } + printf("peeklimit:%d\n", thre); + + for (int i = 0; i < max; i++) { + int p = ct.width() * ct.height() * point[i].z + ct.width() * point[i].y + point[i].x; + if (ct[p] >= thre) { + result.push_back(point[i]); + } + } + + return result; + } +} +array3 Extravasation_detect::dencity_vertex(Raw_image& ct, array3& bone) +{ + std::vector element; + for (int z = -1; z < 2; z++) + for (int y = -1; y < 2; y++) + for (int x = -1; x < 2; x++) + element.push_back(z * ct.data.IMAGESIZE + y * ct.data.Width + x); + + array3 sub(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, ct.image, ct.data.VOXELS); + auto CT = unshape(sub); + array3 CTb2, CTb, CTba, laplacea, laplaceb; + mist::gaussian_filter(CT, CTb, 0.5); + mist::laplacian_filter(sub, laplacea); + std::vector> Points; + float nearlaplace = 0; + array3 region(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + + for (int z = 1; z < ct.data.Depth - 1; z++) + for (int y = 1; y < ct.data.Height - 1; y++) + for (int x = 1; x < ct.data.Width - 1; x++) { + int point = z * ct.data.IMAGESIZE + y * ct.data.Width + x; + if (CTb[point] > 0) { + + bool check = true; + + for (int k = 0; k < 27; k++) { + if (CTb[point] < CTb[point + element[k]]) { + check = false; + break; + } + } + + if (check) { + + Points.push_back(mist::vector3(x, y, z)); + + } + } + } + + //return points; + auto rpoint = reduction_point(sub, Points, 0.1); + + mist::region_growing_utility::less condi(-50); + mist::region_growing_utility::sphere component(1); + mist::region_growing(laplacea, region, rpoint, 255, component, condi); + laplacea.clear(); + Points.clear(); + + array3 label; + auto count = mist::labeling6(region, label); + + //�����ݒ� + auto VenoPointL = std::vector(count + 1); + + VenoPointL[0].Delete = false; + VenoPointL[0].Detectcheck = false; + + for (int i = 0; i < ct.data.Depth; i++) + for (int j = 0; j < ct.data.Height; j++) + for (int x = 0; x < ct.data.Width; x++) { + + //���W + int point = i * ct.data.IMAGESIZE + j * ct.data.Width + x; + if (label[point] != 0) { + VenoPointL[label[point]].Area++; + VenoPointL[label[point]].CTValue += sub[point]; + VenoPointL[label[point]].xCount += x; + VenoPointL[label[point]].yCount += j; + VenoPointL[label[point]].zCount += i; + + if (VenoPointL[label[point]].minz == 0 || VenoPointL[label[point]].minz > i) + VenoPointL[label[point]].minz = i; + if (VenoPointL[label[point]].maxz == 0 || VenoPointL[label[point]].maxz < i) + VenoPointL[label[point]].maxz = i; + } + } + + //���� + for (int i = 1; i < count + 1; i++) { + //������� + if (VenoPointL[i].Area < 20 || VenoPointL[i].Area > 2000 || VenoPointL[i].CTValue / VenoPointL[i].Area < 40) { + VenoPointL[i].Delete = false; + VenoPointL[i].Detectcheck = false; + } + } + + for (int i = 0; i < ct.data.VOXELS; i++) { + if (region[i] != 0) + region[i] = VenoPointL[label[i]].Delete ? region[i] : 0; + } + + CTb.clear(); + label.clear(); + + return region; +} + +//�]�����@ +Extravasation_detect::Result Extravasation_detect::Evaluation(array3& InputLabel, size_t& VenoCount, + array3& Mask,int SIZE1, int SIZE2) +{ + array3 MaskLabel; + size_t TrueCount = mist::labeling6(Mask, MaskLabel); + auto TruePoint = std::vector(TrueCount + 1); + + struct TrueGET + { + float MIND = FAI; + int VenoNumber = 0; + }; + + std::vector Get(TrueCount + 1); + int ImageSize = Mask.width() * Mask.height(); + for (int i = 0; i < Mask.depth(); i++) + for (int j = 0; j < Mask.height(); j++) + for (int x = 0; x < Mask.width(); x++) { + int point = i * ImageSize + j * Mask.width() + x; + if (MaskLabel[point] != 0) { + TruePoint[MaskLabel[point]].Area++; + TruePoint[MaskLabel[point]].xCount += x; + TruePoint[MaskLabel[point]].yCount += j; + TruePoint[MaskLabel[point]].zCount += i; + } + } + + for (int i = 1; i < TrueCount + 1; i++) { + + if (TruePoint[i].Area < 30) { + TruePoint[i].Delete = false; + } + + if (TruePoint[i].Delete) { + TruePoint[i].x = (int)(TruePoint[i].xCount / TruePoint[i].Area); + TruePoint[i].y = (int)(TruePoint[i].yCount / TruePoint[i].Area); + TruePoint[i].z = (int)(TruePoint[i].zCount / TruePoint[i].Area); + } + } + + Result Return; + auto VenoPoint = std::vector(VenoCount + 1); + + for (int i = 0; i < VenoCount + 1; i++) { + VenoPoint[i].Delete = VenoPoint[i].Detectcheck = false; + } + + for (int i = 0; i < Mask.depth(); i++) + for (int j = 0; j < Mask.height(); j++) + for (int x = 0; x < Mask.width(); x++) { + int point = i * ImageSize + j * Mask.width() + x; + if (InputLabel[point] != 0) { + VenoPoint[InputLabel[point]].Area++; + VenoPoint[InputLabel[point]].xCount += x; + VenoPoint[InputLabel[point]].yCount += j; + VenoPoint[InputLabel[point]].zCount += i; + + if (MaskLabel[point] != 0) { + VenoPoint[InputLabel[point]].Delete = true; + VenoPoint[InputLabel[point]].Detectcheck = true; + } + } + } + + int dcount = 0; + + for (int i = 1; i < VenoCount + 1; i++) { + VenoPoint[i].x = (int)(VenoPoint[i].xCount / VenoPoint[i].Area); + VenoPoint[i].y = (int)(VenoPoint[i].yCount / VenoPoint[i].Area); + VenoPoint[i].z = (int)(VenoPoint[i].zCount / VenoPoint[i].Area); + if (VenoPoint[i].Detectcheck) + dcount++; + } + printf("candi:%d\n", dcount); + + + for (int i = 1; i < VenoCount + 1; i++) { + bool TPcheck = false; + for (int j = 1; j < TrueCount + 1; j++) { + + if (VenoPoint[i].Detectcheck) { + float dx = (TruePoint[j].x - VenoPoint[i].x) * Mask.reso1(); + float dy = (TruePoint[j].y - VenoPoint[i].y) * Mask.reso1(); + float dz = (TruePoint[j].z - VenoPoint[i].z) * Mask.reso1(); + + //�����`�F�b�N,����̈�Ƃ̏d�S�ԋ����𒲂ׂ� + if (TruePoint[j].Delete) + if (sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)) < 7.5) { + if (Get[j].MIND == FAI) { + TPcheck = true; + TruePoint[j].Detectcheck = false; + Get[j].VenoNumber = i; + Get[j].MIND = (float)sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)); + } + else if (Get[j].MIND > sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2))) { + Get[j].MIND = (float)sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2)); + Get[j].VenoNumber = i; + } + } + } + } + } + + for (int i = 1; i < VenoCount + 1; i++) { + if (VenoPoint[i].Detectcheck) { + bool TPcheck = false; + for (int j = 1; j < TrueCount + 1; j++) { + if (!TruePoint[j].Detectcheck) { + if (i == Get[j].VenoNumber) + TPcheck = true; + } + } + + //�ǂ̐������x���Ƃ��߂��Ȃ��ꍇ,FP�ɒlj� + if (!TPcheck) { + Return.FP++; + VenoPoint[i].Detectcheck = false; + if (VenoPoint[i].z < SIZE1) { + Return.HFP++; + } + else if (VenoPoint[i].z < SIZE2) { + Return.UFP++; + } + else { + Return.LFP++; + } + + } + } + else { + Return.FP++; + if (VenoPoint[i].z < SIZE1) { + Return.HFP++; + } + else if (VenoPoint[i].z < SIZE2) { + Return.UFP++; + } + else { + Return.LFP++; + } + } + } + + for (int i = 0; i < ImageSize * Mask.depth(); i++) { + + } + + for (int i = 1; i < VenoCount + 1; i++) + if (VenoPoint[i].Detectcheck) { + Return.TP++; + } + + + for (int i = 1; i < TrueCount + 1; i++) { + if (TruePoint[i].Delete) { + Return.OriginalP++; + if (TruePoint[i].z < SIZE1) { + Return.HEADP++; + } + else if (TruePoint[i].z < SIZE2) { + Return.UPPERP++; + } + else { + Return.LOWERP++; + } + + if (TruePoint[i].Detectcheck) { + Return.MissCount++; + if (TruePoint[i].z < SIZE1) { + Return.MissH++; + } + else if (TruePoint[i].z < SIZE2) { + Return.MissU++; + } + else { + Return.MissL++; + } + } + } + } + + Return.region = (int)VenoCount; + + //���ʂ̏o�� + Return.Evalu_Array = array3(Mask.width(), Mask.height(), Mask.depth(), Mask.reso1(), Mask.reso1(), Mask.reso1(), 0); + +#pragma omp parallel for + for (int i = 0; i < InputLabel.size(); i++) { + if (InputLabel[i] != 0) { + if (!VenoPoint[InputLabel[i]].Detectcheck) + Return.Evalu_Array[i] = 255; //FP + else + Return.Evalu_Array[i] = 140; //TP + + } + + if (MaskLabel[i] != 0) + if (TruePoint[MaskLabel[i]].Delete && TruePoint[MaskLabel[i]].Detectcheck) { + Return.Evalu_Array[i] = 0; //FN + } + } + + return(Return); +} + +//���������摜��l�‚Ȃ��ɖ߂� +array3 Extravasation_detect::ArrayLink(array3& FirstArray, array3& SecondArray) +{ + array3 Return(FirstArray.width(), FirstArray.height(), FirstArray.depth() + SecondArray.depth(), 0.5, 0.5, 0.5, 0); + + int VOXELS = (int)(FirstArray.width() * FirstArray.height() * (FirstArray.depth() + SecondArray.depth())); + + for (int i = 0; i < VOXELS; i++) { + if (i < FirstArray.size()) { + Return[i] = FirstArray[i]; + } + else { + Return[i] = SecondArray[i - FirstArray.size()]; + } + } + return(Return); +} + +//���ʂ̕\�� +void Extravasation_detect::Result::ToString(int casenum) +{ + if (casenum != 0) { + printf("�Ǘ�:%d�����o��:%d����̈�̐�:%dTP:%dFP:%d����:%d����:%d�����g:%d�����o�̈�:%d����:%d����:%d�����g:%d\n", casenum, + region, OriginalP, TP, FP, HFP, UFP, LFP, MissCount, MissH, MissU, MissL); + } + else { + printf("�����o��:%d����̈�̐�:%dTP:%dFP:%d����:%d����:%d�����g:%d�����o�̈�:%d����:%d����:%d�����g:%d\n", + region, OriginalP, TP, FP, HFP, UFP, LFP, MissCount, MissH, MissU, MissL); + } +} diff --git a/extravasation_detect.h b/extravasation_detect.h new file mode 100644 index 0000000..ab2dfa6 --- /dev/null +++ b/extravasation_detect.h @@ -0,0 +1,101 @@ +#pragma once +#include "raw_image_class.h" +//#include +//#include +#include "myOpenCV3.h" +#include "mymist.h" + +//�o�����o�Ɋւ���֐����܂Ƃ߂��N���X +//��ɍ����o��@�Əo�����̈撊�o��@ +static class Extravasation_detect +{ +public: + //���W + struct LabelData { + int x = 0; + int y = 0; + int z = 0; + int minz = 0; + int maxz = 0; + int minx = 0; + int maxx = 0; + int miny = 0; + int maxy = 0; + int xCount = 0; + int yCount = 0; + int zCount = 0; + int Area = 0; + int CTValue = 0; + bool Delete = true; + bool Detectcheck = true; + }; + //�ǂ̈ʒu�������Ă��邩��\���񋓌^ + enum BODYTYPE{HEAD,UPPER,LOWER}; + //�����o�֘A + //�l�H������ + static array3 Artificial_Segmentation(Raw_image& ct); + //�K���I臒l���� + static int Adaptive_Threshold(Raw_image& ct); + //watershed + static array3 Watershed(short* ct, array3 prebone, int thre); + + //�]���̍����o��@ + static array3 Old_BoneExtractMethod(Raw_image& ct, BODYTYPE type); + //threshold��ς��邱�ƂœK���I臒l�����ɂȂ� + static array3 Bone_Extract(Raw_image& ct, BODYTYPE type,int threshold = 350); + + //�o�����̈撊�o��@ + /*���̈撊�o�p�����[�^*/ + struct Parameter { + + //�����o + int HighCut; + + //����CT�l + int MaxCTValue; + int MinCTValue; + + //���̈撊�o + int MINV; + int CutVeno; + int ExtractKarnel; + float sigma; + int MaxVolume; + }; + //��i�K臒l���� + static array3 double_thresholding(Raw_image& ct, array3& bone, BODYTYPE type); + + //�s�[�N�_�x�[�X(dencity_vertex�����O�ɍ��̈揜�����K�v) + template + static array3 unshape(array3& input); + static std::vector> reduction_point(array3 ct, std::vector>& point, float rate); + static array3 dencity_vertex(Raw_image& ct,array3& bon); + + //�]�� + //���ʂ̍\���� + struct Result { + array3 Evalu_Array; + int region = 0; + int TP = 0; + int FP = 0; + int MissCount = 0; + int OriginalP = 0; + int HEADP = 0; + int LOWERP = 0; + int UPPERP = 0; + int MissH = 0; + int MissL = 0; + int MissU = 0; + int HFP = 0; + int UFP = 0; + int LFP = 0; + + public: + void ToString(int casenum = 0); + }; + static Result Evaluation(array3& InputLabel, size_t& VenoCount, + array3& Mask, int SIZE1, int SIZE2); + + //���̑� + static array3 ArrayLink(array3& FirstArray, array3& SecondArray); +}; diff --git a/function.h b/function.h new file mode 100644 index 0000000..0926deb --- /dev/null +++ b/function.h @@ -0,0 +1,219 @@ +#pragma once +#include "Pass_management.h" +#include "extravasation_detect.h" +#include "raw_image_class.h" +#include +void OutputName(char* Bone, char* Ev, char* PreEv, char* Res, int casenum) { + + sprintf_s(Bone, 256, DIR CASENUM OUTNAME THREBON RAW, casenum, casenum); + sprintf_s(Ev, 256, DIR CASENUM OUTNAME NAME_DOUBLE EVA RAW, casenum, casenum); + sprintf_s(PreEv, 256, DIR CASENUM OUTNAME NAME_DOUBLE PREEVA RAW, casenum, casenum); + sprintf_s(Res, 256, DIR CASENUM OUTNAME NAME_DOUBLE RES RAW, casenum, casenum); + +} + + +void EvansDetect(Raw_image& ct, array3& evans, array3& bone, int size1, int size2) { + Raw_image CTHEAD; + Raw_image CTLOWER; + Raw_image CTUPPER; + + CTHEAD.data.Depth = size1; + CTHEAD.data.VOXELS = CTHEAD.data.IMAGESIZE * size1; + CTHEAD.image = new short[CTHEAD.data.VOXELS]; + + CTUPPER.data.Depth = size2 - size1; + CTUPPER.data.VOXELS = CTUPPER.data.IMAGESIZE * CTUPPER.data.Depth; + CTUPPER.image = new short[CTUPPER.data.VOXELS]; + + CTLOWER.data.Depth = ct.data.Depth - size2; + CTLOWER.data.VOXELS = CTLOWER.data.IMAGESIZE * CTLOWER.data.Depth; + CTLOWER.image = new short[CTLOWER.data.VOXELS]; + + for (int i = 0; i < ct.data.VOXELS; i++) { + if (i < ct.data.IMAGESIZE * size1) { + CTHEAD.image[i] = ct.image[i]; + } + else if (i < ct.data.IMAGESIZE * size2) { + CTUPPER.image[i - (ct.data.IMAGESIZE * size1)] = ct.image[i]; + } + else { + CTLOWER.image[i - (ct.data.IMAGESIZE * size2)] = ct.image[i]; + } + } + + int thre = 350; + if (size2 != 0) { + thre = Extravasation_detect::Adaptive_Threshold(CTUPPER); + printf("thre:%d\n", thre); + } + + array3 BoneH, BoneU, BoneL; + array3 EvansH, EvansL, EvansU; + + if (CTHEAD.data.Depth != 0) { + auto th1 = std::thread([&BoneH, &CTHEAD,&EvansH, thre] + { + switch (BE) + { + case OLD: + BoneH =Extravasation_detect::Old_BoneExtractMethod(CTHEAD,Extravasation_detect::HEAD); + break; + case WATER: + BoneH =Extravasation_detect::Bone_Extract(CTHEAD,Extravasation_detect::HEAD); + break; + case ADAPTIVE: + BoneH =Extravasation_detect::Bone_Extract(CTHEAD,Extravasation_detect::HEAD,thre); + break; + default: + break; + } + + EvansH =Extravasation_detect::double_thresholding(CTHEAD, BoneH, Extravasation_detect::HEAD); + }); + + auto th2 = std::thread([&BoneU, &CTUPPER, &EvansU, thre] + { + switch (BE) + { + case OLD: + BoneU =Extravasation_detect::Old_BoneExtractMethod(CTUPPER,Extravasation_detect::UPPER); + break; + case WATER: + BoneU =Extravasation_detect::Bone_Extract(CTUPPER, Extravasation_detect::UPPER); + break; + case ADAPTIVE: + BoneU =Extravasation_detect::Bone_Extract(CTUPPER,Extravasation_detect::UPPER,thre); + break; + default: + break; + } + + + + EvansU =Extravasation_detect::double_thresholding(CTUPPER, BoneU,Extravasation_detect::UPPER); + }); + + auto th3 = std::thread([&BoneL, &CTLOWER,&EvansL, thre] + { + switch (BE) + { + case OLD: + BoneL =Extravasation_detect::Old_BoneExtractMethod(CTLOWER,Extravasation_detect::LOWER); + break; + case WATER: + BoneL =Extravasation_detect::Bone_Extract(CTLOWER,Extravasation_detect::LOWER); + break; + case ADAPTIVE: + BoneL =Extravasation_detect::Bone_Extract(CTLOWER, Extravasation_detect::LOWER,thre); + break; + default: + break; + } + + EvansL =Extravasation_detect::double_thresholding(CTLOWER, BoneL,Extravasation_detect::LOWER); + }); + + th1.join(); + th2.join(); + th3.join(); + + array3 tmp = Extravasation_detect::ArrayLink(BoneH, BoneU); + bone = Extravasation_detect::ArrayLink(tmp, BoneL); + + array3 tmpE = Extravasation_detect::ArrayLink(EvansH, EvansU); + evans = Extravasation_detect::ArrayLink(tmpE, EvansL); + + BoneH.clear(); + BoneU.clear(); + BoneL.clear(); + + EvansH.clear(); + EvansU.clear(); + EvansL.clear(); + + } + else if (CTUPPER.data.Depth != 0) { + + auto th2 = std::thread([&BoneU, &CTUPPER, &EvansU, thre] + { + switch (BE) + { + case OLD: + BoneU =Extravasation_detect::Old_BoneExtractMethod(CTUPPER,Extravasation_detect::UPPER); + break; + case WATER: + BoneU =Extravasation_detect::Bone_Extract(CTUPPER,Extravasation_detect::UPPER); + break; + case ADAPTIVE: + BoneU =Extravasation_detect::Bone_Extract(CTUPPER,Extravasation_detect::UPPER,thre); + break; + default: + break; + } + EvansU =Extravasation_detect::double_thresholding(CTUPPER, BoneU,Extravasation_detect::UPPER); + }); + + auto th3 = std::thread([&BoneL, &CTLOWER,&EvansL, thre] + { + switch (BE) + { + case OLD: + BoneL =Extravasation_detect::Old_BoneExtractMethod(CTLOWER,Extravasation_detect::LOWER); + break; + case WATER: + BoneL =Extravasation_detect::Bone_Extract(CTLOWER,Extravasation_detect::LOWER); + break; + case ADAPTIVE: + BoneL =Extravasation_detect::Bone_Extract(CTLOWER,Extravasation_detect::LOWER,thre); + break; + default: + break; + } + EvansL =Extravasation_detect::double_thresholding(CTLOWER, BoneL,Extravasation_detect::LOWER); + }); + + th2.join(); + th3.join(); + + bone = Extravasation_detect::ArrayLink(BoneU, BoneL); + evans = Extravasation_detect::ArrayLink(EvansU, EvansL); + + BoneU.clear(); + BoneL.clear(); + EvansU.clear(); + EvansL.clear(); + } + else { + auto th3 = std::thread([&BoneL, &CTLOWER,&EvansL] + { + switch (BE) + { + case OLD: + BoneL =Extravasation_detect::Old_BoneExtractMethod(CTLOWER,Extravasation_detect::LOWER); + break; + case WATER: + BoneL =Extravasation_detect::Bone_Extract(CTLOWER,Extravasation_detect::LOWER); + break; + case ADAPTIVE: + BoneL =Extravasation_detect::Bone_Extract(CTLOWER,Extravasation_detect::LOWER); + break; + default: + break; + } + EvansL =Extravasation_detect::double_thresholding(CTLOWER, BoneL,Extravasation_detect::LOWER); + }); + + th3.join(); + + bone = BoneL; + evans = EvansL; + BoneL.clear(); + EvansL.clear(); + + } + + delete[] CTHEAD.image; + delete[] CTUPPER.image; + delete[] CTLOWER.image; +} diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..3d9cbdc --- /dev/null +++ b/main.cpp @@ -0,0 +1,758 @@ +#include +#include +#include +#include "raw_image_class.h" +#include "Pass_management.h" +#include +#include +#include "slice_predict.h" +#include "function.h" +#include "Feature_extraction.h" +#define Leng 64 +#define USE_MODEL true +enum { RTREE, SVM }; + +//using namespace cv; +using namespace std; + +int main(int argc, char* argv[]) { + int mode_select; + printf("���o���[�h:1 �̈�����ʒ��o���[�h:2 �w�K���[�h:3 FROC��̓��[�h:4 BoneSegmentation:5"); + scanf("%d", &mode_select); + + if (mode_select == 1) { + int startcase; + printf("�J�n�Ǘ�:"); + scanf("%d", &startcase); + + int endcase; + printf("�I��:"); + scanf("%d", &endcase); + + cv::Mat RegionInfo; + + auto AStart = std::chrono::steady_clock::now(); + double EvaluationTime = 0; + + FILE* fpOut; + FILE* fpOutP; + + + fopen_s(&fpOut, DIR DOUBLEOUT OUTNAME CSV, "w"); + fopen_s(&fpOutP, DIR DOUBLEPRE OUTNAME CSV, "w"); + + if (fpOut == NULL) { + printf("�t�@�C�����J���܂���\n"); + getchar(); + + } + + fprintf(fpOut, "�Ǘ�ԍ�,���̈搔,�o�����ʐ�,����,�㔼�g,�����g,TP��,FP��,����,�㔼�g,�����g,�����o�̈搔,����,�㔼�g,�����g,���o����(s),�X���C�X��,\n"); + fprintf(fpOutP, "�Ǘ�ԍ�,���̈搔,�o�����ʐ�,����,�㔼�g,�����g,TP��,FP��,����,�㔼�g,�����g,�����o�̈搔,����,�㔼�g,�����g,���o����(s),�X���C�X��,\n"); + + //cv::Ptr r_tree = cv::ml::RTrees::load("D:\\YL\\Data\\BleedingData\\Ktest\\rtree_model.xml"); + cv::Ptr r_tree = cv::ml::RTrees::load(DIR MODEL XML); + cv::dnn::Net net = cv::dnn::readNet(DIR CLASS); + cv::dnn::Net snet = cv::dnn::readNet(DIR SHOULDER); + cv::dnn::Net lnet = cv::dnn::readNet(DIR WAIST); + + printf(DIR WAIST); + + + + bool Calib = false; + + for (int i = startcase; i < endcase; i++) { + //�o�͐�̐ݒ� + char Bon[256]; + char EV[256]; + char PRE[256]; + char RE[256]; + + OutputName(Bon, EV, PRE, RE, i); + char Name[256]; + sprintf_s(Name, 256, DIR CTI RAW, i, i); + Raw_image ct(Name); + + char MaskName[256]; + sprintf_s(MaskName, 256, DIR MAK RAW, i, i); + Raw_image mask(MaskName); + printf("Detection start"); + if (!ct.empty() && !mask.empty()) { + + printf("�Ǘ�:%d\n", i); + for (int sli = 0; sli < mask.data.VOXELS; sli++) { + if (mask.image[sli] != 0) + mask.image[sli] = 255; + else + mask.image[sli] = 0; + } + + FILE* fp; + char NaSIZE[256]; + sprintf_s(NaSIZE, 256, DIR SIZEFILE CSV, i, i); + fopen_s(&fp, NaSIZE, "r"); + + if (fp == NULL) { + printf("csv not exist---------------------------------------"); + return(false); + } + int SIZE1 = 0; + int SIZE2 = 0; + fscanf(fp, "%d,%d,%f", &SIZE1, &SIZE2, &ct.data.reso); + fclose(fp); + printf("%d,%d,%f\n", SIZE1, SIZE2, ct.data.reso); + + Slice_predict::predict(ct, SIZE1, SIZE2, net, snet, lnet); + printf("%d,%d,%f\n", SIZE1, SIZE2, ct.data.reso); + array3 DetectEvans; + array3 Bone; + + auto start = std::chrono::steady_clock::now(); + + EvansDetect(ct, DetectEvans, Bone, SIZE1, SIZE2); + + //�v���I�� + auto end = std::chrono::steady_clock::now(); + double time = (double)std::chrono::duration_cast(end - start).count(); + printf("���o����%f\n", (time) / 1000); + + array3 MArray(ct.data.Height, ct.data.Width, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + for (int sli = 0; sli < ct.data.VOXELS; sli++) + MArray[sli] = mask.image[sli]; + + //���ʂ̃��x�����O + array3 Label; + size_t ResCount = mist::labeling6(DetectEvans, Label); + + auto evastart = std::chrono::steady_clock::now(); + + Extravasation_detect::Result Res = Extravasation_detect::Evaluation(Label, ResCount, MArray, SIZE1, SIZE2); + Res.ToString(i); + + auto evaend = std::chrono::steady_clock::now(); + double evatime = (double)std::chrono::duration_cast(evaend - evastart).count(); + printf("�]������:%f\n", evatime / 1000); + EvaluationTime += evatime; + + fprintf(fpOut, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%0.3f,%d\n", i, Res.region, Res.OriginalP, + Res.HEADP, Res.UPPERP, Res.LOWERP, Res.TP, Res.FP, Res.HFP, Res.UFP, Res.LFP, Res.MissCount, + Res.MissH, Res.MissU, Res.MissL, time / 1000, ct.data.Depth); + + mist::write_raw(Bone, Bon); + mist::write_raw(DetectEvans, RE); + mist::write_raw(Res.Evalu_Array, EV); + + if (USE_MODEL) { + auto prestart = std::chrono::steady_clock::now(); + cv::Mat FestureMat = Feature_extraction::All_feature_exstract(ct, Label, ResCount, Leng); + array3 Pre = Feature_extraction::Predict(r_tree, Label, ResCount, FestureMat); + auto preend = std::chrono::steady_clock::now(); + double pretime = (double)std::chrono::duration_cast(preend - prestart).count(); + printf("���o����%f\n", (time + pretime) / 1000); + + //���ʂ̃��x�����O + array3 PreLabel; + size_t PreCount = mist::labeling6(Pre, PreLabel); + + evastart = std::chrono::steady_clock::now(); + Extravasation_detect::Result PredictRes = Extravasation_detect::Evaluation(PreLabel, PreCount, MArray, SIZE1, SIZE2); + evaend = std::chrono::steady_clock::now(); + evatime = (double)std::chrono::duration_cast(evaend - evastart).count(); + printf("�]������:%f\n", evatime / 1000); + EvaluationTime += evatime; + mist::write_raw(PredictRes.Evalu_Array, PRE); + + PredictRes.ToString(i); + + fprintf(fpOutP, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%0.3f,%d\n", i, PredictRes.region, PredictRes.OriginalP, + PredictRes.HEADP, PredictRes.UPPERP, PredictRes.LOWERP, PredictRes.TP, PredictRes.FP, PredictRes.HFP, PredictRes.UFP, PredictRes.LFP, PredictRes.MissCount, + PredictRes.MissH, PredictRes.MissU, PredictRes.MissL, (pretime + time) / 1000, ct.data.Depth); + } + + delete[] ct.image; + delete[] mask.image; + printf("-----------------------------------------------------------------------------------------------------------------------------\n"); + + } + + } + + auto Aend = std::chrono::steady_clock::now(); + double Atime = (double)std::chrono::duration_cast(Aend - AStart).count(); + + printf("���v����:%f[s](�����]������%f[s])", Atime / 1000, EvaluationTime / 1000); + + fclose(fpOut); + fclose(fpOutP); + + } + else if (mode_select == 2) { + printf("�̈�����ʒ��o�J�n...\n"); + auto AStart = std::chrono::steady_clock::now(); + cv::Mat TrueSet; + cv::Mat DataSet; + cv::Mat TrueRegion; + bool Calib = false; + bool trueCalib = false; + + int TP = 0; + int FP = 0; + int TTP = 0; + int MTP = 0; + + for (int i = 0; i < 109; i++) { + char Name[256]; + sprintf_s(Name, 256, DIR CTI RAW, i, i); + Raw_image ct(Name); + + char MaskName[256]; + sprintf_s(MaskName, 256, DIR MAK RAW, i, i); + Raw_image mask(MaskName); + + if (!ct.empty() && !mask.empty()) { + for (int i = 0; i < ct.data.VOXELS; i++) { + if (mask.image[i] != 0) + mask.image[i] = 255; + else + mask.image[i] = 0; + } + + printf("�Ǘ�%d...\n", i); + + array3 MArray(ct.data.Height, ct.data.Width, ct.data.Depth, 0.5, 0.5, 0.5, 0); + + for (int i = 0; i < ct.data.VOXELS; i++) { + MArray[i] = mask.image[i]; + } + + array3 MaskLabel; + size_t TrueCount = mist::labeling6(MArray, MaskLabel); + auto TruePoint = std::vector(TrueCount + 1); + + for (int i = 0; i < MArray.depth(); i++) + for (int j = 0; j < MArray.height(); j++) + for (int x = 0; x < MArray.width(); x++) { + int point = i * MArray.width() * MArray.height() + j * MArray.width() + x; + if (MaskLabel[point] != 0) { + TruePoint[MaskLabel[point]].Area++; + } + } + + for (int i = 0; i < MArray.size(); i++) { + + if (TruePoint[MaskLabel[i]].Area < 20) { + MArray[i] = 0; + } + } + + + MaskLabel.clear(); + array3 TrueHE = Feature_extraction::Histogram_equalization(ct, Leng); + cv::Mat TL; + + for (int i = -1; i <= 1; i++) { + array3 cMarray(ct.data.Height, ct.data.Width, ct.data.Depth, 0.5, 0.5, 0.5, 0); + for (int j = 0; j < MArray.size(); j++) { + cMarray[j] = MArray[j]; + } + + if (i < 0) { + mist::erosion(cMarray, -1 * i, 0); + printf("Erosion..."); + } + if (i > 0) { + mist::dilation(cMarray, i, 0); + printf("Dilation..."); + } + + array3 Label; + auto ResCount = mist::labeling6(cMarray, Label); + + if (ResCount > 0) { + cv::Mat TFesture = Feature_extraction::All_feature_exstract(ct, Label, ResCount, Leng); + + + TL = Feature_extraction::LabelMat(TFesture, 1); + + if (!Calib) { + DataSet = TL.clone(); + TrueRegion = TL.clone(); + Calib = true; + } + else { + DataSet = Feature_extraction::Addsample(DataSet, TL); + TrueRegion = Feature_extraction::Addsample(TrueRegion, TL); + } + printf("%d,", DataSet.rows); + MTP += TL.rows; + } + } + + printf("row:%d\n", MTP); + + char ResultName[256]; + sprintf_s(ResultName, 256, DIR CASENUM "DOUBLEFeature%d" RAW, i, i); + Raw_image Resa(ResultName); + + if (!Resa.empty()) { + array3 RFArray(ct.data.Height, ct.data.Width, ct.data.Depth, 0.5, 0.5, 0.5, 0); + array3 RTArray(ct.data.Height, ct.data.Width, ct.data.Depth, 0.5, 0.5, 0.5, 0); + + for (int i = 0; i < ct.data.VOXELS; i++) { + if (Resa.image[i] > 10) { + if (Resa.image[i] == 255) { + RFArray[i] = 255; + } + else { + RTArray[i] = 255; + } + } + else { + RFArray[i] = 0; + RTArray[i] = 0; + } + } + + cv::Mat RFL, RTL; + + auto th2 = std::thread([&ct, &RFArray, &TrueHE, &RFL] { + array3 Label; + auto ResCount = mist::labeling6(RFArray, Label); + cv::Mat FalseFesture = Feature_extraction::All_feature_exstract(ct, Label, ResCount, Leng); + RFL = Feature_extraction::LabelMat(FalseFesture, 0); }); + + auto th3 = std::thread([&ct, &RTArray, &TrueHE, &RTL] { + + array3 Label; + auto ResCount = mist::labeling6(RTArray, Label); + if (ResCount > 0) { + cv::Mat TFesture = Feature_extraction::All_feature_exstract(ct, Label, ResCount, Leng); + RTL = Feature_extraction::LabelMat(TFesture, 1); + }}); + + th2.join(); + th3.join(); + + + if (RFL.rows > 0) { + DataSet = Feature_extraction::Addsample(DataSet, RFL); + } + if (RTL.rows > 0) { + DataSet = Feature_extraction::Addsample(DataSet, RTL); + } + + if (!trueCalib) { + TrueSet = Feature_extraction::Addsample(RFL, RTL); + trueCalib = true; + } + else { + if (RFL.rows > 0) { + TrueSet = Feature_extraction::Addsample(TrueSet, RFL); + } + if (RTL.rows > 0) { + TrueSet = Feature_extraction::Addsample(TrueSet, RTL); + } + } + + TP += RTL.rows; + FP += RFL.rows; + TTP += TL.rows; + printf("�Ǘ�%d��͊��� �o���̈搔:%d TP:%d FP:%d\n", i, TL.rows, RTL.rows, RFL.rows); + printf("Dataset col:%d rows:%d\n", DataSet.cols, DataSet.rows); + Calib = true; + trueCalib = true; + delete[] Resa.image; + } + + delete[] ct.image; + delete[] mask.image; + + printf("-----------------------------------------------------------------------------------\n"); + + + } + } + + auto Aend = std::chrono::steady_clock::now(); + double Atime = (double)std::chrono::duration_cast(Aend - AStart).count(); + + printf("���v����:%f[s]", Atime / 1000); + printf("TP:%d FP:%d\n", TP, FP); + printf("TrueSet�F%d\n", TrueSet.rows); + printf("TrueRegion�F%d\n", TrueRegion.rows); + cv::Mat Re = Feature_extraction::Feature_Cut(DataSet, true, 1); + + char csvfile[256]; + sprintf_s(csvfile, 256, DIR DST CSV); + Feature_extraction::writeCSV(csvfile, Re); + + char Rcsv[256]; + sprintf_s(Rcsv, 256, DIR RCSV CSV); + Feature_extraction::writeCSV(Rcsv, TrueSet); + + } + else if (mode_select == 3) { + int random; + printf("�����_���t�H���X�g:1 SVM:2"); + scanf("%d", &random); + if (random == 1) { + printf("�����_���t�H���X�g�w�K�J�n\n"); + printf(DIR DST CSV); + cv::Ptr data = cv::ml::TrainData::loadFromCSV(DIR DST CSV, 0, 0, -1); + cv::Ptr r_tree = cv::ml::RTrees::create(); + printf("Traning DataSet:%s\n", DST CSV); + //data->setTrainTestSplitRatio(0.95, true); + + r_tree->setMaxDepth(23); + r_tree->setMinSampleCount(2); + r_tree->setCalculateVarImportance(true); + r_tree->setTermCriteria(cv::TermCriteria(cv::TermCriteria::MAX_ITER + cv::TermCriteria::EPS, 1000, 0.001)); + r_tree->train(data); + + printf("train error: %f\n", r_tree->calcError(data, false, cv::noArray())); + /* printf("test error: %f\n\n", r_tree->calcError(data, true, noArray()));*/ + cv::Mat s = r_tree->getVarImportance(); + + printf("row:%d,col:%d\n", s.rows, s.cols); + + for (int i = 0; i < s.rows; i++) { + for (int k = 0; k < s.cols; k++) { + printf("%f", s.at(i, k)); + if (k == s.cols - 1) { + printf("\n"); + } + else + printf(","); + } + } + + r_tree->save(DIR MODEL XML); + + } + else if (random == 2) { + cv::Ptr data = cv::ml::TrainData::loadFromCSV(DIR DST CSV, 0, 0, -1); + cv::Ptr testdata = cv::ml::TrainData::loadFromCSV(DIR RCSV CSV, 0, 0, -1); + //data->setTrainTestSplitRatio(0.90, true); + + printf("Traning DataSet:%s\n", DST CSV); + int adaC; + double adaGanma; + double maxScore = 100; + for (int C = 1; C < 10; C++) + for (double G = 1e-5; G < 10000; G *= 10) { + cv::Ptr svm = cv::ml::SVM::create(); + svm->setType(cv::ml::SVM::C_SVC); + svm->setKernel(cv::ml::SVM::RBF); + svm->setC(C); + svm->setGamma(G); + svm->setTermCriteria(cv::TermCriteria(cv::TermCriteria::MAX_ITER + cv::TermCriteria::EPS, 1000, 0.001)); + svm->train(data); + + float score = svm->calcError(testdata, false, cv::noArray()); + + //printf("Score:%f C:%f Gamma:%f\n", score, C, G); + if (maxScore > score) { + printf("Updata....Score:%f C:%d Gamma:%f\n", score, C, G); + maxScore = score; + adaC = C; + adaGanma = G; + } + } + + printf("final....Score:%f C:%d Gamma:%f\n", maxScore, adaC, adaGanma); + + cv::Ptr svm = cv::ml::SVM::create(); + svm->setType(cv::ml::SVM::C_SVC); + svm->setKernel(cv::ml::SVM::RBF); + svm->setC(adaC); + svm->setGamma(adaGanma); + svm->setTermCriteria(cv::TermCriteria(cv::TermCriteria::MAX_ITER + cv::TermCriteria::EPS, 1000, 0.001)); + svm->train(data); + + svm->save(DIR "\\SVMIso" XML); + } + printf("\a\a\a"); + printf("�w�K����\n"); + } + else if (mode_select == 4) { + printf("Rtree:0 SVM:1\n"); + int mode; + scanf("%d", &mode); + if (mode == RTREE) { + cv::Ptr r_tree = cv::ml::RTrees::load(DIR MODEL XML); + printf("model_name:%s\n", MODEL XML); + printf("ROCCSV:%s\n", RCSV CSV); + char name[256]; + sprintf_s(name, 256, DIR RCSV CSV); + Feature_extraction::CSVData s = Feature_extraction::readCSV(name); + float* pre = new float[s.Labels.rows]; + + FILE* fpOut; + fopen_s(&fpOut, DIR FROC CSV, "w"); + fprintf(fpOut, "Rate,TP,TN,FP,FN,���x,����FP\n"); + + for (float Rate = 0.6; Rate > 0; Rate -= (float)0.01) { + int TP = 0; + int TN = 0; + int FP = 0; + int FN = 0; + + for (int i = 0; i < s.Labels.rows; i++) { + cv::Mat sam = s.samples.row(i).clone(); + pre[i] = r_tree->predict(sam, cv::noArray(), cv::ml::StatModel::RAW_OUTPUT); + pre[i] = pre[i] / 1000; + + if (pre[i] > Rate) + pre[i] = 1; + else + pre[i] = 0; + + if (pre[i] == 1) { + if (pre[i] == s.Labels.at(i, 0)) { + TP++; + } + else { + FP++; + } + } + else { + if (pre[i] == s.Labels.at(i, 0)) { + FN++; + } + else { + TN++; + } + } + + } + printf("Rate:%f TP:%d TN:%d FP:%d FN:%d ���x:%f ����FP:%f\n", Rate, TP, TN, FP, FN, (float)TP / (float)(TP + TN), (float)FP / 28); + fprintf(fpOut, "%f,%d,%d,%d,%d,%f,%f\n", Rate, TP, TN, FP, FN, (float)TP / (float)(TP + TN), (float)FP / 35); + } + fclose(fpOut); + + fopen_s(&fpOut, DIR ROC TXT, "w"); + + for (int i = 0; i < s.Labels.rows; i++) { + cv::Mat sam = s.samples.row(i).clone(); + pre[i] = r_tree->predict(sam, cv::noArray(), cv::ml::StatModel::RAW_OUTPUT); + pre[i] = pre[i] / 1000; + + fprintf(fpOut, "%f,%d\n", pre[i], (int)s.Labels.at(i, 0)); + } + + fclose(fpOut); + } + else if (mode == SVM) { + cv::Ptr r_tree = cv::ml::SVM::load(DIR "\\SVMIso" XML); + printf("model_name:%s\n", "\\SVMIso" XML); + printf("ROCCSV:%s\n", RCSV CSV); + char name[256]; + sprintf_s(name, 256, DIR RCSV CSV); + Feature_extraction::CSVData s = Feature_extraction::readCSV(name); + double* pre = new double[s.Labels.rows]; + + + FILE* fpOut; + fopen_s(&fpOut, DIR FROC CSV, "w"); + fprintf(fpOut, "Rate,TP,TN,FP,FN,���x,����FP\n"); + + + for (float j = 0; j <= 1; j += 0.001) { + + int TP = 0; + int TN = 0; + int FP = 0; + int FN = 0; + for (int i = 0; i < s.Labels.rows; i++) { + cv::Mat sam = s.samples.row(i).clone(); + pre[i] = r_tree->predict(sam); + //printf("%f\n", pre[i]); + if (pre[i] <= j) { + if (1 == s.Labels.at(i, 0)) { + TP++; + } + else { + FP++; + } + } + else { + if (0 == s.Labels.at(i, 0)) { + FN++; + } + else { + TN++; + } + } + } + + printf("%f TP:%d TN:%d FP:%d FN:%d ���x:%f ����FP:%f\n", j, TP, TN, FP, FN, (float)TP / (float)(TP + TN), (float)FP / 35); + fprintf(fpOut, "%f,%d,%d,%d,%d,%f,%f\n", j, TP, TN, FP, FN, (float)TP / (float)(TP + TN), (float)FP / 35); + } + fclose(fpOut); + fopen_s(&fpOut, DIR ROC TXT, "w"); + + for (int i = 0; i < s.Labels.rows; i++) { + cv::Mat sam = s.samples.row(i).clone(); + pre[i] = r_tree->predict(sam, cv::noArray()); + + fprintf(fpOut, "%f,%d\n", pre[i], (int)s.Labels.at(i, 0)); + } + + fclose(fpOut); + } + + + } + else if (mode_select == 5) + { + int startcase; + printf("�J�n�Ǘ�:"); + scanf("%d", &startcase); + + int endcase; + printf("�I��:"); + scanf("%d", &endcase); + + cv::Mat RegionInfo; + cv::Ptr r_tree = cv::ml::RTrees::load(DIR MODEL XML); + cv::dnn::Net net = cv::dnn::readNet(DIR CLASS); + cv::dnn::Net snet = cv::dnn::readNet(DIR SHOULDER); + cv::dnn::Net lnet = cv::dnn::readNet(DIR WAIST); + auto AStart = std::chrono::steady_clock::now(); + //double EvaluationTime = 0; + + bool Calib = false; + + for (int i = startcase; i < endcase; i++) { + //�o�͐�̐ݒ� + char Bon[256]; + char EV[256]; + char PRE[256]; + char RE[256]; + + OutputName(Bon, EV, PRE, RE, i); + char Name[256]; + sprintf_s(Name, 256, DIR CTI RAW, i, i); + Raw_image ct(Name); + + /*char MaskName[256]; + sprintf_s(MaskName, 256, DIR MAK RAW, i, i); + Raw_image mask(MaskName);*/ + printf("Segmentation start"); + //if (!ct.empty() && !mask.empty()) { + if (!ct.empty()) { + printf("�Ǘ�:%d\n", i); + /*for (int sli = 0; sli < mask.data.VOXELS; sli++) { + if (mask.image[sli] != 0) + mask.image[sli] = 255; + else + mask.image[sli] = 0; + }*/ + + //FILE* fp; + //char NaSIZE[256]; + //sprintf_s(NaSIZE, 256, DIR SIZEFILE CSV, i, i); + //fopen_s(&fp, NaSIZE, "r"); + + //if (fp == NULL) { + // printf("csv not exist---------------------------------------"); + // return(false); + //} + int SIZE1 = 0; + int SIZE2 = 0; + /*fscanf(fp, "%d,%d,%f", &SIZE1, &SIZE2, &ct.data.reso); + fclose(fp); + printf("%d,%d,%f\n", SIZE1, SIZE2, ct.data.reso);*/ + + Slice_predict::predict(ct, SIZE1, SIZE2, net, snet, lnet); + //printf("%d,%d,%f\n", SIZE1, SIZE2, ct.data.reso); + array3 DetectEvans; + array3 Bone; + + auto start = std::chrono::steady_clock::now(); + + EvansDetect(ct, DetectEvans, Bone, SIZE1, SIZE2); + + //�v���I�� + /*auto end = std::chrono::steady_clock::now(); + double time = (double)std::chrono::duration_cast(end - start).count(); + printf("���o����%f\n", (time) / 1000); + + array3 MArray(ct.data.Height, ct.data.Width, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0); + for (int sli = 0; sli < ct.data.VOXELS; sli++) + MArray[sli] = mask.image[sli];*/ + + ////���ʂ̃��x�����O + //array3 Label; + //size_t ResCount = mist::labeling6(DetectEvans, Label); + + //auto evastart = std::chrono::steady_clock::now(); + + //Extravasation_detect::Result Res = Extravasation_detect::Evaluation(Label, ResCount, MArray, SIZE1, SIZE2); + //Res.ToString(i); + + //auto evaend = std::chrono::steady_clock::now(); + //double evatime = (double)std::chrono::duration_cast(evaend - evastart).count(); + //printf("�]������:%f\n", evatime / 1000); + //EvaluationTime += evatime; + + /*fprintf(fpOut, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%0.3f,%d\n", i, Res.region, Res.OriginalP, + Res.HEADP, Res.UPPERP, Res.LOWERP, Res.TP, Res.FP, Res.HFP, Res.UFP, Res.LFP, Res.MissCount, + Res.MissH, Res.MissU, Res.MissL, time / 1000, ct.data.Depth);*/ + + mist::write_raw(Bone, Bon); + /*mist::write_raw(DetectEvans, RE); + mist::write_raw(Res.Evalu_Array, EV);*/ + + //if (USE_MODEL) { + // auto prestart = std::chrono::steady_clock::now(); + // cv::Mat FestureMat = Feature_extraction::All_feature_exstract(ct, Label, ResCount, Leng); + // array3 Pre = Feature_extraction::Predict(r_tree, Label, ResCount, FestureMat); + // auto preend = std::chrono::steady_clock::now(); + // double pretime = (double)std::chrono::duration_cast(preend - prestart).count(); + // printf("���o����%f\n", (time + pretime) / 1000); + + // //���ʂ̃��x�����O + // array3 PreLabel; + // size_t PreCount = mist::labeling6(Pre, PreLabel); + + // evastart = std::chrono::steady_clock::now(); + // Extravasation_detect::Result PredictRes = Extravasation_detect::Evaluation(PreLabel, PreCount, MArray, SIZE1, SIZE2); + // evaend = std::chrono::steady_clock::now(); + // evatime = (double)std::chrono::duration_cast(evaend - evastart).count(); + // printf("�]������:%f\n", evatime / 1000); + // EvaluationTime += evatime; + // mist::write_raw(PredictRes.Evalu_Array, PRE); + + // PredictRes.ToString(i); + + // /*fprintf(fpOutP, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%0.3f,%d\n", i, PredictRes.region, PredictRes.OriginalP, + // PredictRes.HEADP, PredictRes.UPPERP, PredictRes.LOWERP, PredictRes.TP, PredictRes.FP, PredictRes.HFP, PredictRes.UFP, PredictRes.LFP, PredictRes.MissCount, + // PredictRes.MissH, PredictRes.MissU, PredictRes.MissL, (pretime + time) / 1000, ct.data.Depth);*/ + //} + + delete[] ct.image; + /*delete[] mask.image;*/ + printf("-----------------------------------------------------------------------------------------------------------------------------\n"); + + } + + } + + //auto Aend = std::chrono::steady_clock::now(); + //double Atime = (double)std::chrono::duration_cast(Aend - AStart).count(); + + //printf("���v����:%f[s](�����]������%f[s])", Atime / 1000, EvaluationTime / 1000); + + //fclose(fpOut); + //fclose(fpOutP); + } + printf("����"); + + printf("\a"); + getchar(); + + + return(0); +} diff --git a/myOpenCV3.h b/myOpenCV3.h new file mode 100644 index 0000000..b1232e1 --- /dev/null +++ b/myOpenCV3.h @@ -0,0 +1,109 @@ + +// OpenCV 3�n ���ʃw�b�_�[�t�@�C�� +// T.Nakaguchi + +#pragma once + +// �w�b�_�[�t�@�C�� +#pragma warning(disable: 4819) +#include +#pragma warning(default: 4819) + +// �o�[�W�����擾 +#define CV_VERSION_STR CVAUX_STR(CV_MAJOR_VERSION) CVAUX_STR(CV_MINOR_VERSION) CVAUX_STR(CV_SUBMINOR_VERSION) + +// �r���h���[�h +#ifdef _DEBUG +#define CV_EXT_STR "d.lib" +#else +#define CV_EXT_STR ".lib" +#endif + +// ���C�u�����̃����N�i�s�v�ȕ��̓R�����g�A�E�g�j +#define PRE_COMPILE 0 // �C���X�g�[���łȂ� 1 �•ʃ��C�u�����g�p���� 0 +#define PREHEAD "opencv_" + +#if PRE_COMPILE +// OpenCV3.0 �C���X�g�[���� +#pragma comment(lib, PREHEAD "world" CV_VERSION_STR CV_EXT_STR) // �S�� +#pragma comment(lib, PREHEAD "ts" CV_VERSION_STR CV_EXT_STR) // �f���֘A + +#else +// �•ʂ̃��C�u�����w�� �iCmake��Static�Ŏg�p���Ȃǁj +// ��{���W���[�� +#pragma comment(lib, PREHEAD "world" CV_VERSION_STR CV_EXT_STR) // ��{�@�\ +//#pragma comment(lib, PREHEAD "core" CV_VERSION_STR CV_EXT_STR) // ��{�@�\ +//#pragma comment(lib, PREHEAD "imgproc" CV_VERSION_STR CV_EXT_STR) // �摜���� +//#pragma comment(lib, PREHEAD "imgcodecs" CV_VERSION_STR CV_EXT_STR) // �摜�t�@�C�����o�� +//#pragma comment(lib, PREHEAD "videoio" CV_VERSION_STR CV_EXT_STR) // ����t�@�C�����o�� +//#pragma comment(lib, PREHEAD "highgui" CV_VERSION_STR CV_EXT_STR) // ���@�\GUI +//#pragma comment(lib, PREHEAD "video" CV_VERSION_STR CV_EXT_STR) // ���摜��� +//#pragma comment(lib, PREHEAD "calib3d" CV_VERSION_STR CV_EXT_STR) // �J�����Z���ƎO�����č\�z +//#pragma comment(lib, PREHEAD "features2d" CV_VERSION_STR CV_EXT_STR) // �摜������� +//#pragma comment(lib, PREHEAD "objdetect" CV_VERSION_STR CV_EXT_STR) // ���̌��o +//#pragma comment(lib, PREHEAD "ml" CV_VERSION_STR CV_EXT_STR) // �@�B�w�K +//#pragma comment(lib, PREHEAD "flann" CV_VERSION_STR CV_EXT_STR) // �������N���X�^�����O�ƌ��� +//#pragma comment(lib, PREHEAD "photo" CV_VERSION_STR CV_EXT_STR) // �v�Z�@�ʐ^ +//#pragma comment(lib, PREHEAD "stitching" CV_VERSION_STR CV_EXT_STR) // �摜�ڑ� +//#pragma comment(lib, PREHEAD "hal" CV_VERSION_STR CV_EXT_STR) // �n�[�h�E�F�A������ +//#pragma comment(lib, PREHEAD "shape" CV_VERSION_STR CV_EXT_STR) // �`���v���o +//#pragma comment(lib, PREHEAD "superres" CV_VERSION_STR CV_EXT_STR) // ���� +//#pragma comment(lib, PREHEAD "videostab" CV_VERSION_STR CV_EXT_STR) // ���摜���艻 +//#pragma comment(lib, PREHEAD "vis" CV_VERSION_STR CV_EXT_STR) // 3�����Ž��� +// �g�����W���[�� +//#pragma comment(lib, PREHEAD "adas" CV_VERSION_STR CV_EXT_STR) // ��i�I�f�o�C�X�T�|�[�g +//#pragma comment(lib, PREHEAD "aruco" CV_VERSION_STR CV_EXT_STR) // AR�p�}�[�J�[ +//#pragma comment(lib, PREHEAD "bgsegm" CV_VERSION_STR CV_EXT_STR) // ���nj^�w�i�E�O�i���� +//#pragma comment(lib, PREHEAD "bioinspired" CV_VERSION_STR CV_EXT_STR) // ���̂Ɋ�Â����o�I���� +//#pragma comment(lib, PREHEAD "ccalib" CV_VERSION_STR CV_EXT_STR) // �J�X�^���p�^�[���ɂ��J�����Z���ƎO�����č\�� +//#pragma comment(lib, PREHEAD "cvv" CV_VERSION_STR CV_EXT_STR) // �Θb�I���o�I�f�o�b�OGUI +//#pragma comment(lib, PREHEAD "datasets" CV_VERSION_STR CV_EXT_STR) // ����f�[�^�Z�b�g��舵���t���[�����[�N +//#pragma comment(lib, PREHEAD "dnn" CV_VERSION_STR CV_EXT_STR) // DNN�T�|�[�g +//#pragma comment(lib, PREHEAD "face" CV_VERSION_STR CV_EXT_STR) // ��F�� +//#pragma comment(lib, PREHEAD "latentsvm" CV_VERSION_STR CV_EXT_STR) // Latent-SVM +//#pragma comment(lib, PREHEAD "line_descriptor" CV_VERSION_STR CV_EXT_STR) // �����o�̃o�C�i���\�� +//#pragma comment(lib, PREHEAD "matlab" CV_VERSION_STR CV_EXT_STR) // MATLAB�u���b�W +//#pragma comment(lib, PREHEAD "optflow" CV_VERSION_STR CV_EXT_STR) // �I�v�e�B�J���t���[ +//#pragma comment(lib, PREHEAD "reg" CV_VERSION_STR CV_EXT_STR) // �摜�ʒu���킹 +//#pragma comment(lib, PREHEAD "rgbd" CV_VERSION_STR CV_EXT_STR) // RGB-�[�x�J���� +//#pragma comment(lib, PREHEAD "saliency" CV_VERSION_STR CV_EXT_STR) // �摜 ������ API +//#pragma comment(lib, PREHEAD "surface_matching" CV_VERSION_STR CV_EXT_STR) // �\�ʃ��f����v���o +//#pragma comment(lib, PREHEAD "text" CV_VERSION_STR CV_EXT_STR) // �V�[���������o�ƔF�� +//#pragma comment(lib, PREHEAD "tracking" CV_VERSION_STR CV_EXT_STR) // �ǐ� +//#pragma comment(lib, PREHEAD "xfeatures2d" CV_VERSION_STR CV_EXT_STR) // �g���� �摜������� +//#pragma comment(lib, PREHEAD "ximgproc" CV_VERSION_STR CV_EXT_STR) // �g���� �摜���� +//#pragma comment(lib, PREHEAD "xobjdetect" CV_VERSION_STR CV_EXT_STR) // �g���� ���̌��o +//#pragma comment(lib, PREHEAD "xphoto" CV_VERSION_STR CV_EXT_STR) // �g���� �v�Z�@�ʐ^ +#endif + +//using namespace cv; + + +// �����摜��A���\�� +//void imShowMulti(String winname, std::vector& imgs, unsigned int cols, unsigned int rows, Size imgsize, unsigned int border) +//{ +// if (imgs.size() < 1 || cols < 1 || rows < 1) return; +// +// unsigned int w = imgsize.width + border, h = imgsize.height + border; +// Mat board(h * rows + border, w * cols + border, CV_8UC3, CV_RGB(128, 128, 128)); +// for (unsigned int r = 0, i = 0; r < rows; r ++) { +// for(unsigned int c = 0; c < cols; c ++, i ++) { +// Rect roi_rect = Rect(c * w + border, r * h + border, imgsize.width, imgsize.height); +// Mat roi(board, roi_rect); +// if (i < imgs.size()) { +// if (imgs[i].type() == CV_8UC3) { +// resize(imgs[i], roi, imgsize); +// } else if (imgs[i].type() == CV_8UC1) { +// Mat c3; +// cvtColor(imgs[i], c3, COLOR_GRAY2BGR); +// resize(c3, roi, imgsize); +// } else { +// putText(roi, "Color mode not matched", Point(20,20), CV_FONT_HERSHEY_COMPLEX, 0.5, CV_RGB(0,0,0)); +// } +// } else { +// putText(roi, "No image", Point(20,20), CV_FONT_HERSHEY_COMPLEX, 0.5, CV_RGB(0,0,0)); +// } +// } +// } +// imshow(winname, board); +//} diff --git a/mymist.h b/mymist.h new file mode 100644 index 0000000..2c74cd4 --- /dev/null +++ b/mymist.h @@ -0,0 +1,31 @@ +#ifndef __MYMIST_H +#define __MYMIST_H + +#include "mist/mist.h" +#include "mist/io/raw.h" +#include "mist/timer.h" +#include "mist/filter/morphology.h" +#include "mist/filter/linear.h" +#include "mist/filter/thinning.h" +#include "mist\config\mist_conf.h" +#include "mist\filter\mode.h" +#include "mist/filter/region_growing.h" +#include "mist/filter/median.h" +#include "mist/filter/labeling.h" +//#include +//#include +//#include +//#include "mist/io/dicom.h" +//#include "mist/io/dicom_info.h" +//#include +//#include +#include "mist/filter/boundary.h" +//#include +//#include +//#include +//#include +//#include "mist/vector.h" +//#include +#endif //__MYMIST_H + +using namespace mist; diff --git a/raw_image_class.cpp b/raw_image_class.cpp new file mode 100644 index 0000000..50052ea --- /dev/null +++ b/raw_image_class.cpp @@ -0,0 +1,67 @@ +#include "raw_image_class.h" +template +inline Raw_image::Raw_image(char* filename) +{ + + //�摜�ǂݍ��� + FILE* file; + fopen_s(&file, filename, "rb"); + if (file == NULL) { + image = 0; + return; + } + + // �t�@�C���T�C�Y�m�F + fseek(file, 0, SEEK_END); + int fileSize = ftell(file); + int depth = fileSize / (data.Width * data.Height * sizeof(Type)); + + if (depth != data.Depth) { + if (data.Depth == 0) { + data.Depth = depth; + data.VOXELS = data.Width * data.Height * data.Depth; + } + else { + data.Depth = depth > data.Depth ? data.Depth : depth; + data.VOXELS = data.Width * data.Height * data.Depth; + } + } + else { + data.Depth = depth; + data.VOXELS = data.Width * data.Height * data.Depth; + } + + image = new Type[data.VOXELS]; + + // �f�[�^�ǂݍ��� + fseek(file, 0, SEEK_SET); + fread(image, sizeof(Type), data.VOXELS, file); + fclose(file); + + return; +} + +template +Raw_image::Raw_image() +{ +} + +template +Raw_image::~Raw_image() +{ + if (image != 0) { + delete[] image; + } + return; +} + +template +bool Raw_image::empty() +{ + if (image == 0) { + return false; + } + else { + return true; + } +} diff --git a/raw_image_class.h b/raw_image_class.h new file mode 100644 index 0000000..7c50738 --- /dev/null +++ b/raw_image_class.h @@ -0,0 +1,91 @@ +#pragma once +#include +//#include +//#include +#include "myOpenCV3.h" +#include "mymist.h" + +//Raw�摜�̓ǂݍ��݂��s���N���X +//�Ȃ񂩃f�X�g���N�^�ŃG���[���o���̂ō��͒��ډ�����Ă���(�����炭function.h�̏o�����o�̊֐�������) +template +class Raw_image +{ +public: + //Raw�摜��3�������̕ۑ��\���� + struct DataInfo { + int Width = 512; + int Height = 512; + int IMAGESIZE = Width * Height; + int Depth = 0; + int VOXELS = 0; + float reso = 0.5; + }; + DataInfo data; + + //�摜�z�� + Type* image = NULL; + + //�R���X�g���N�^�ƃf�X�g���N�^ + Raw_image(char* filename) + { + //�摜�ǂݍ��� + FILE* file; + fopen_s(&file, filename, "rb"); + if (file == NULL) { + image = NULL; + return; + } + + // �t�@�C���T�C�Y�m�F + fseek(file, 0, SEEK_END); + int fileSize = ftell(file); + int depth = fileSize / (data.Width * data.Height * sizeof(Type)); + + if (depth != data.Depth) { + if (data.Depth == 0) { + data.Depth = depth; + data.VOXELS = data.Width * data.Height * data.Depth; + } + else { + data.Depth = depth > data.Depth ? data.Depth : depth; + data.VOXELS = data.Width * data.Height * data.Depth; + } + } + else { + data.Depth = depth; + data.VOXELS = data.Width * data.Height * data.Depth; + } + + image = new Type[data.VOXELS]; + + // �f�[�^�ǂݍ��� + fseek(file, 0, SEEK_SET); + fread(image, sizeof(Type), data.VOXELS, file); + fclose(file); + + return; + } + Raw_image() + { + } + ~Raw_image() + { + if (image != NULL) { + //delete[] image; + } + return; + } + + //�ǂݍ��߂Ă��邩�̊m�F + bool empty() + { + if (image == NULL) { + return true; + } + else { + return false; + } + } +}; + + diff --git a/slice_predict.cpp b/slice_predict.cpp new file mode 100644 index 0000000..d7c902b --- /dev/null +++ b/slice_predict.cpp @@ -0,0 +1,175 @@ +#include "slice_predict.h" + +//MIP�摜�̐��� +cv::Mat Slice_predict::createMIP(Raw_image& images) +{ + cv::Mat MIP(images.data.Depth, images.data.Width, CV_16S, cv::Scalar::all(0)); + + short* pt = MIP.ptr(0); + for (int i = 0; i < images.data.Width * images.data.Depth; i++) { + short max = 0; + for (int j = 0; j < images.data.Height; j++) { + int point = (i / images.data.Width) * images.data.IMAGESIZE + j * images.data.Width + (i % images.data.Width); + if (max < images.image[point]) + max = images.image[point]; + } + pt[i] = max; + } + double grad = 255.0 / 1700.0; + cv::Mat dispMIP; + MIP.convertTo(dispMIP, CV_8U, grad, 127.5 - (1000.0 * grad)); + cvtColor(dispMIP.clone(), dispMIP, cv::COLOR_GRAY2BGR); + return dispMIP; +} + +//���֌W���̌v�Z +std::vector Slice_predict::Correlation(std::vector& input1, std::vector& input2) +{ + int range = input2.size() / 2; + std::vector result; + float ave2 = 0; + //printf("range:%d\n", range); + for (int i = range; i < input1.size() - range; i++) { + //input1 Average; + float ave1 = 0; + for (int j = i - range; j <= i + range; j++) { + ave1 += input1[j]; + } + ave1 = ave1 / input2.size(); + float Sxy = 0; + float Sx = 0, Sy = 0; + for (int j = i - range, k = 0; j < i + range; j++, k++) { + Sxy += (input1[j] - ave1) * (input2[k] - ave2); + Sx += (input1[j] - ave1) * (input1[j] - ave1); + Sy += (input2[k] - ave2) * (input2[k] - ave2); + } + Sxy = Sxy / input2.size(); + Sx = sqrt(Sx / input2.size()); + Sy = sqrt(Sy / input2.size()); + result.push_back(Sxy / (Sx * Sy)); + } + return result; +} + +//�w�K�������f����ǂݍ���ŃX���C�X�ʒu�𐄒� +void Slice_predict::predict(Raw_image& image, int& shoulder, int& low, cv::dnn::Net& class_model, cv::dnn::Net& shoulde_model, cv::dnn::Net& lower_model) +{ + + cv::Mat im = createMIP(image); + cv::Mat imr; + resize(im, imr, cv::Size(224, 224)); + cv::Mat blob; + class_model.setInput(cv::dnn::blobFromImage(imr, 1.0 / 255.0, cv::Size(224, 224), cv::Scalar(), true)); + cv::Mat prob = class_model.forward(); + auto p = prob.ptr(0); + float max = 0; + int arg = 0; + for (int j = 0; j < 4; j++) + if (max < p[j]) { + max = p[j]; + arg = j; + } + + if (arg == 0) { + cv::Mat ims(im.rows + 512, im.cols, CV_8UC3, cv::Scalar::all(0)); + cv::Mat q(ims, cv::Rect(0, 256, im.cols, im.rows)); + im.copyTo(q); + + std::vector liner = { -60, -30, 0 , 30 , 60 }; + int range = liner.size() / 2; + std::vector imgs; + for (int j = 0; j < im.rows - 1; j += 30) { + cv::Mat input(ims, cv::Rect(0, j, im.cols, im.cols)); + imgs.push_back(input); + } + + std::vector predS(imgs.size()); + + shoulde_model.setInput(cv::dnn::blobFromImages(imgs, 1.0 / 255.0, cv::Size(224, 224), cv::Scalar(), true)); + cv::Mat pred = shoulde_model.forward(); + for (int j = 0; j < pred.rows; j++) { + auto p = pred.ptr(0); + predS[j] = p[j]; + } + + auto correration = Correlation(predS, liner); + float min = 9999; + int place = 0; + + for (int j = 0; j < correration.size(); j++) { + auto lengthabs = abs(predS[j + range]); + if (correration[j] > 0.98 && lengthabs < 15) { + if (min > lengthabs) { + min = lengthabs; + place = j + range; + } + } + } + printf("%d,%d,%d\n", correration.size(), predS.size(), place); + int splace; + if (place != 0) { + splace = (place * 30) - predS[place]; + shoulder = splace; + } + else { + splace = shoulder = 0; + } + + imgs = std::vector(); + std::vector predL; + for (int j = splace; j < im.rows - 1; j += 30) { + cv::Mat input(ims, cv::Rect(0, j, im.cols, im.cols)); + imgs.push_back(input); + } + lower_model.setInput(cv::dnn::blobFromImages(imgs, 1.0 / 255.0, cv::Size(224, 224), cv::Scalar(), true)); + cv::Mat predLs = lower_model.forward(); + auto p = predLs.ptr(0); + for (int j = 0; j < predLs.rows; j++) { + predL.push_back(p[j]); + } + + auto cL = Correlation(predL, liner); + min = 9999; + place = im.rows; + for (int j = 0; j < cL.size(); j++) { + auto lengthabs = abs(predL[j + range]); + if (cL[j] > 0.98 && lengthabs < 15) { + if (min > lengthabs) { + min = lengthabs; + place = j + range; + } + } + } + + int lplace; + if (place != im.rows) { + lplace = (place * 30) - predL[place] + splace; + low = lplace; + } + else { + lplace = im.rows; + low = im.rows; + } + printf("Soulder:%d Lower:%d\n", splace, lplace); + return; + } + else { + switch (arg) + { + case 1: + shoulder = im.rows; + low = im.rows; + break; + case 2: + shoulder = 0; + low = 0; + break; + case 3: + shoulder = 0; + low = im.rows; + default: + break; + } + } + +} diff --git a/slice_predict.h b/slice_predict.h new file mode 100644 index 0000000..44641a0 --- /dev/null +++ b/slice_predict.h @@ -0,0 +1,27 @@ +#pragma once +#include "raw_image_class.h" +#include +//#include +//#include +#include "myOpenCV3.h" +#include "mymist.h" +//using namespace cv; + +//MIP�摜����̌��y�э��̃X���C�X�ʒu����Ɋւ���N���X +//DNN���W���[����p����keras�Ŋw�K�������f����ǂݍ��ݐ���(.pb�t�@�C���ւ̕ϊ����K�v) +static class Slice_predict +{ + //MIP�摜���� + static cv::Mat createMIP(Raw_image& images); + + //���֌W���̌v�Z + static std::vector Correlation(std::vector& input1, std::vector& input2); + +public: + //�X���C�X�ʒu���� + static void predict(Raw_image& image, int& shoulder, int& low, + cv::dnn::Net& class_model, cv::dnn::Net& shoulde_model, + cv::dnn::Net& lower_model); + +}; +