#include "extravasation_detect.h"
//using namespace cv;
#define MAXW 400
#define MINW 2
#define MINA 600
#define THRESHOLDCENTER 425
#define FAI -999
//骨抽出関連
array3<uchar> Extravasation_detect::Artificial_Segmentation(Raw_image<short>& ct)
{
auto Artficial = mist::array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0);
auto Artf = mist::array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0);
auto result = mist::array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0);
std::vector<mist::vector3<int>> 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<int>(x, y, z));
}
}
}
mist::region_growing_utility::equal<int> equal(180);
mist::region_growing_utility::sphere component(1);
array3<uchar> 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<short>& ct)
{
array3<uchar> 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<short> Upper = array3<short>(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<uchar> Extravasation_detect::Watershed(short* ct, array3<uchar> prebone, int thre)
{
array3<uchar> 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<uchar> thinf(512, 512, prebone.reso1(), prebone.reso2(), (uchar*)front.ptr<uchar>(0), sizeof(uchar) * 512 * 512);
mist::hilditch::thinning(thinf, thinf);
uchar* pf = front.ptr<uchar>(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<short>(0);
uchar* imwp = image_w.ptr<uchar>(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<cv::Vec3b> 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<int>(i);
cv::Vec3b* pix = Dst.ptr<cv::Vec3b>(i);
uchar* rix = DstU.ptr<uchar>(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<uchar>(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;
}
array3<uchar>Extravasation_detect::Old_BoneExtractMethod(Raw_image<short>& ct, BODYTYPE type)
{
short* Bone = new short[ct.data.VOXELS];
array3<uchar> 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<uchar> 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<uchar>(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<unsigned int>(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<LabelData> CTV = std::vector<LabelData>(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++) {
//検出結果からラベリングされた領域を発見
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<mist::vector3<int>> 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) {
//判定条件外のラベルはfalseとしておく
BoneCheck[i] = false;
bonedetect++;
//}
}
else {
Points.push_back(mist::vector3<int>(CTV[i].x, CTV[i].y, CTV[i].z));
bonedetect++;
}
}
auto BoneDelete = array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0);
auto BoneOut = array3<uchar>(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<uchar>(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);
//外輪郭抽出
std::vector<std::vector<cv::Point>> 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<uchar>(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<int> get(180);
mist::region_growing_utility::sphere component(1);
array3<uchar> 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<uchar>(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<uchar>(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<uchar>(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);
//外輪郭抽出
std::vector<std::vector<cv::Point>> 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<uchar>(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<uchar> Extravasation_detect::Bone_Extract(Raw_image<short>& ct, BODYTYPE type, int threshold)
{
short* Bone = new short[ct.data.VOXELS];
array3<uchar> 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<uchar> 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<uchar>(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<unsigned int>(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<LabelData> CTV = std::vector<LabelData>(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++) {
//検出結果からラベリングされた領域を発見
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<mist::vector3<int>> 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<int>(CTV[i].x, CTV[i].y, CTV[i].z));
bonedetect++;
}
}
auto BoneDelete = array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso, 0);
auto BoneOut = array3<uchar>(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<uchar>(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);
//外輪郭抽出
std::vector<std::vector<cv::Point>> 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<uchar>(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<uchar>(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<uchar>(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<uchar>(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);
//外輪郭抽出
std::vector<std::vector<cv::Point>> 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<uchar>(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);
}
//候補領域抽出関連
array3<uchar> Extravasation_detect::double_thresholding(Raw_image<short>& ct, array3<uchar>& bone, BODYTYPE type)
{
array3<uchar> 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値調整
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;
}
//候補領域抽出
array3<uchar> VenoOnlyL = array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso);
array3<uchar> VenoOnlyH = array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso);
array3<uchar> LV = array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso);
array3<uchar> HV = array3<uchar>(ct.data.Width, ct.data.Height, ct.data.Depth, ct.data.reso, ct.data.reso, ct.data.reso);
//静脈時相格納用のコンテナ
//静脈のカットオフ
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;
//ガウシアンフィルタ
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<short>(0);
short* p2 = getblurH.ptr<short>(0);
for (int i = 0; i < ct.data.IMAGESIZE; i++, p++, p2++) {
//座標
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;
}
}
//誤検出除去
//ラベリング
array3<unsigned int> VenoLabelL;
array3<unsigned int> VenoLabelH;
size_t VenoCountL = mist::labeling6(VenoOnlyL, VenoLabelL);
size_t VenoCountH = labeling6(VenoOnlyH, VenoLabelH);
//初期設定
auto VenoPointL = std::vector<LabelData>(VenoCountL + 1);
auto VenoPointH = std::vector<LabelData>(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++) {
//座標
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<typename Type>
inline array3<float> Extravasation_detect::unshape(array3<Type>& input)
{
array3<float> CTb, laplace;
mist::gaussian_filter(input, CTb);
array3<float> 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<mist::vector3<int>> Extravasation_detect::reduction_point(array3<short> ct, std::vector<mist::vector3<int>>& point, float rate)
{
{
float max = (float)point.size();
int hist[600];
std::vector<mist::vector3<int>> 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<uchar> Extravasation_detect::dencity_vertex(Raw_image<short>& ct, array3<uchar>& bone)
{
std::vector<int> 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<short> 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<float> CTb2, CTb, CTba, laplacea, laplaceb;
mist::gaussian_filter(CT, CTb, 0.5);
mist::laplacian_filter(sub, laplacea);
std::vector<mist::vector3<int>> Points;
float nearlaplace = 0;
array3<uchar> 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<int>(x, y, z));
}
}
}
//return points;
auto rpoint = reduction_point(sub, Points, 0.1);
mist::region_growing_utility::less<float> condi(-50);
mist::region_growing_utility::sphere component(1);
mist::region_growing(laplacea, region, rpoint, 255, component, condi);
laplacea.clear();
Points.clear();
array3<unsigned int> label;
auto count = mist::labeling6(region, label);
//初期設定
auto VenoPointL = std::vector<LabelData>(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++) {
//座標
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<unsigned int>& InputLabel, size_t& VenoCount,
array3<uchar>& Mask,int SIZE1, int SIZE2)
{
array3<unsigned int> MaskLabel;
size_t TrueCount = mist::labeling6(Mask, MaskLabel);
auto TruePoint = std::vector<LabelData>(TrueCount + 1);
struct TrueGET
{
float MIND = FAI;
int VenoNumber = 0;
};
std::vector<TrueGET> 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<LabelData>(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();
//距離チェック,正解領域との重心間距離を調べる
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;
}
}
//どの正解ラベルとも近くない場合,FPに追加
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;
//結果の出力
Return.Evalu_Array = array3<uchar>(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);
}
//分割した画像を人つなぎに戻す
array3<uchar> Extravasation_detect::ArrayLink(array3<uchar>& FirstArray, array3<uchar>& SecondArray)
{
array3<uchar> 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総検出数:%d正解領域の数:%dTP:%dFP:%d頭部:%d胴体:%d下半身:%d未検出領域:%d頭部:%d胴体:%d下半身:%d\n", casenum,
region, OriginalP, TP, FP, HFP, UFP, LFP, MissCount, MissH, MissU, MissL);
}
else {
printf("総検出数:%d正解領域の数:%dTP:%dFP:%d頭部:%d胴体:%d下半身:%d未検出領域:%d頭部:%d胴体:%d下半身:%d\n",
region, OriginalP, TP, FP, HFP, UFP, LFP, MissCount, MissH, MissU, MissL);
}
}
