- •1. Программные системы и методы 3D-реконструкции биомедицинских данных
- •1.1. Обзор программных комплексов
- •1.2. Обзор основных подходов
- •1.2.1. Структуры визуализируемых объёмов
- •1.2.2. Сегментация и идентификация объектов
- •1.2.3. Построение поверхностей объектов
- •2. Модели, методы и алгоритмы, положенные в основу сегментации и поиска объектов
- •2.1. Формальная постановка задачи и основные этапы её решения
- •2.2. Сегментация данных компьютерной томографии и электронной микроскопии
- •2.3. Повышение производительности вычислений в методе k-средних
- •2.4. Метод иерархической декомпозиции сегментов
- •2.5. Морфологическая фильтрация сегментов
- •2.6. Построение признакового описания сегмента
- •2.7. Поиск объектов методом динамического программирования
- •3. Описание реализации программной системы
- •3.1. Основные классы, структуры и методы
- •3.2. Организация пользовательского интерфейса
- •Заключение
- •Список литературы
- •Приложение А. Примеры результатов сегментации и идентификации объектов
- •Приложение Б. Классы, реализующие семейство методов k-средних
- •Приложение В. Исходные коды базовых компонентов программной системы
- •В.1. Основные классы и структуры для хранения и обработки данных
- •В.2. Методы вычисления атрибутов сегментов
- •В.3. Поиск объектов по найденным сегментам
В.2. Методы вычисления атрибутов сегментов
#include "Structures.h"
using namespace cv; using namespace std;
/* Data – исходные данные для сегментации, Segment – сегмент из индексов вокселов, принадлежащих исходным данным, LayerIndex – индекс рассматриваемого слоя данных, ParentSegmentIndex и SegmentIndex однозначно идентифицируют рассматриваемый сегмент в иерархической структуре, SegmentType – тип сегмента, в зависимости от рассматриваемого уровня иерархии */
Листинг B.3: Средние значения плотности и координат вокселов
float GetMeanDensity(TVoxelsData* Data, vector <size_t> Segment)
{
float result = 0.0f;
vector |
<size_t>::iterator iter = Segment.begin(); |
vector |
<size_t>::iterator end = Segment.end(); |
for (; |
iter != end; ++iter) result += Data->Density[*iter]; |
return |
result/Segment.size(); |
} |
|
Point2f |
GetMeanCoords_2D(TVoxelsData* Data, vector <size_t> Segment) |
{ |
|
Point2f result(0.0f, 0.0f);
vector |
<size_t>::iterator |
iter = Segment.begin(); |
vector |
<size_t>::iterator |
end = Segment.end(); |
for (; |
iter != end; ++iter) |
{
size_t x = *iter%Data->sizeX,
y = ((*iter-x)/Data->sizeX)%Data->sizeY; result += Point2f((float)x, (float)y);
}
return (1.0f/Segment.size())*result;
}
Point3f GetMeanCoords_3D(TVoxelsData* Data, vector <size_t> Segment)
{
Point3f result(0.0f, 0.0f, 0.0f);
vector <size_t>::iterator iter = Segment.begin(); vector <size_t>::iterator end = Segment.end(); for (; iter != end; ++iter)
{
size_t x = *iter%Data->sizeX,
y= ((*iter-x)/Data->sizeX)%Data->sizeY,
z= (*iter-Data->sizeX*y)/(Data->sizeX*Data->sizeY); result += Point3f((float)x, (float)y, (float)z);
}
return (1.0f/Segment.size())*result;
}
Point3f GetMeanDensityCoords_2D(TVoxelsData* Data, vector <size_t> Segment)
{
Point3f result(0.0f, 0.0f, 0.0f);
vector <size_t>::iterator iter = Segment.begin(); vector <size_t>::iterator end = Segment.end();
68
for (; iter != end; ++iter)
{
size_t x = *iter%Data->sizeX,
y = ((*iter-x)/Data->sizeX)%Data->sizeY;
result += Point3f((float)Data->Density[*iter], (float)x, (float)y);
}
return (1.0f/Segment.size())*result;
}
Vec4f GetMeanDensityCoords_3D(TVoxelsData* Data, vector <size_t> Segment)
{
Vec4f result(0.0f, 0.0f, 0.0f, 0.0f);
vector <size_t>::iterator iter = Segment.begin(); vector <size_t>::iterator end = Segment.end(); for (; iter != end; ++iter)
{
size_t x = *iter%Data->sizeX,
y= ((*iter-x)/Data->sizeX)%Data->sizeY,
z= (*iter-x-Data->sizeX*y)/(Data->sizeX*Data->sizeY);
result += Vec4f((float)Data->Density[*iter], (float)x, (float)y, (float)z);
}
return (1.0f/Segment.size())*result;
}
Листинг B.4: Минимальное и максимальное значения плотности вокселов
Point_<short> MinMaxDensityOfSegment_2D(TVoxelsData* Data, size_t LayerIndex, Vec2i ParentSegmentIndex, int SegmentIndex, int SegmentType)
{
short result_min = Data->MaxDensity, result_max = Data->MinDensity;
for (size_t i = 0; i < Data->sizeX*Data->sizeY; ++i)
{
bool value;
switch (SegmentType)
{
case DENSITY:
{
value = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == SegmentIndex;
}
break;
case SPATIAL:
{
value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) && (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == SegmentIndex);
} break;
case CONNECTED:
{
value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) &&
((Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == ParentSegmentIndex[1]) || (ParentSegmentIndex[1] < 0)) &&
69
(Data->VoxelSegments[i+LayerIndex*Data->sizeX* Data->sizeY].ComponentIndex_2D == SegmentIndex);
} break;
default: break;
}
if (value)
{
result_min = min<short>(Data->Density[i+LayerIndex*Data->sizeX*Data->sizeY], result_min);
result_max = max<short>(Data->Density[i+LayerIndex*Data->sizeX*Data->sizeY], result_max);
}
}
return Point_<short>(result_min, result_max);
}
void CalcMinMaxDensityOfSegments_2D(TVoxelsData* Data, size_t LayerIndex, Vec2i ParentSegmentIndex, int SegmentType, size_t SegmentsCount, vector < Point_<short> > &MinMaxDensity)
{
MinMaxDensity.clear();
for (size_t k = 0; k < SegmentsCount; ++k) MinMaxDensity.push_back(Point_<short>(Data->MaxDensity, Data->MinDensity));
for (size_t i = 0; i < Data->sizeX*Data->sizeY; ++i)
{
int index; bool value;
switch (SegmentType)
{
case DENSITY:
{
index = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D; value = true;
} break;
case SPATIAL:
{
index = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D; value = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0];
} break;
case CONNECTED:
{
index = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].ComponentIndex_2D; value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) && ((ParentSegmentIndex[1] >= 0) && (Data->VoxelSegments[i+
LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == ParentSegmentIndex[1]) || (ParentSegmentIndex[1] < 0));
} break;
default: break;
70
}
if (index >= 0 && value)
{
MinMaxDensity[index].x = min<short>(Data->Density[i+ LayerIndex*Data->sizeX*Data->sizeY], MinMaxDensity[index].x); MinMaxDensity[index].y = max<short>(Data->Density[i+ LayerIndex*Data->sizeX*Data->sizeY], MinMaxDensity[index].y);
}
}
}
Point_<short> MinMaxDensityOfSegment_3D(TVoxelsData* Data, size_t StartLayerIndex, int SegmentIndex, int SegmentType)
{
short result_min = Data->MaxDensity, result_max = Data->MinDensity;
for (size_t i = 0; i < Data->TotalSize – StartLayerIndex*Data->sizeX*Data->sizeY; ++i)
{
bool value;
switch (SegmentType)
{
case DENSITY: value = Data->VoxelSegments[i+ StartLayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_3D == SegmentIndex; break;
case SPATIAL: value = Data->VoxelSegments[i+ StartLayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_3D == SegmentIndex; break;
default: break;
}
if (value)
{
result_min = min<short>(Data->Density[i+ StartLayerIndex*Data->sizeX*Data->sizeY], result_min); result_max = max<short>(Data->Density[i+ StartLayerIndex*Data->sizeX*Data->sizeY], result_max);
}
}
return Point_<short>(result_min, result_max);
}
void CalcMinMaxDensityOfSegments_3D(TVoxelsData* Data, size_t StartLayerIndex, int SegmentType, size_t SegmentsCount, vector <Point_<short> > &MinMaxDensity)
{
MinMaxDensity.clear();
for (size_t k = 0; k < SegmentsCount; ++k) MinMaxDensity.push_back(Point_<short>(Data->MaxDensity, Data->MinDensity));
for (size_t i = 0; i < Data->TotalSize – StartLayerIndex*Data->sizeX*Data->sizeY; ++i)
{
int index;
switch (SegmentType)
{
case DENSITY: index = Data->VoxelSegments[i+
71
StartLayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_3D; break;
case SPATIAL: index = Data->VoxelSegments[i+ StartLayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_3D; break;
default: break;
}
if (index >= 0)
{
MinMaxDensity[index].x = min<short>(Data->Density[i+ StartLayerIndex*Data->sizeX*Data->sizeY], MinMaxDensity[index].x); MinMaxDensity[index].y = max<short>(Data->Density[i+ StartLayerIndex*Data->sizeX*Data->sizeY], MinMaxDensity[index].y);
}
}
}
Листинг B.5: Среднее значение и среднеквадратическое отклонение плотности
Point2f DensityMeanDevOfSegment_2D(TVoxelsData* Data, size_t LayerIndex, Vec2i ParentSegmentIndex, int SegmentIndex, int SegmentType)
{
float sqr = 0.0f, mean = 0.0f; int N = 0;
for (size_t i = 0; i < Data->sizeX*Data->sizeY; ++i)
{
bool value;
switch (SegmentType)
{
case DENSITY:
{
value = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == SegmentIndex;
}
break;
case SPATIAL:
{
value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) && (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == SegmentIndex);
} break;
case CONNECTED:
{
value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) && ((Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == ParentSegmentIndex[1]) || (ParentSegmentIndex[1] < 0)) && (Data->VoxelSegments[i+LayerIndex*Data->sizeX*Data->sizeY].ComponentIndex_2D
== SegmentIndex);
}break;
default: break;
}
if (value)
{
short density = Data->Density[i+LayerIndex*Data->sizeX*Data->sizeY]; sqr += density*density; mean += density;
72
N++;
}
}
return Point2f(mean/N, sqrt(sqr/N-(mean/N)*(mean/N)));
}
void CalcDensityMeanDevOfSegments_2D(TVoxelsData* Data, size_t LayerIndex, Vec2i ParentSegmentIndex, int SegmentType, size_t SegmentsCount, vector <Point2f> &DensityMeanDev)
{
int * N = new int [SegmentsCount]; DensityMeanDev.clear();
for (size_t k = 0; k < SegmentsCount; ++k)
{
DensityMeanDev.push_back(Point2f(0.0f, 0.0f)); N[k] = 0;
}
for (size_t i = 0; i < Data->sizeX*Data->sizeY; ++i)
{
int index; bool value;
switch (SegmentType)
{
case DENSITY:
{
index = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D; value = true;
} break;
case SPATIAL:
{
index = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D; value = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0];
} break;
case CONNECTED:
{
index = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].ComponentIndex_2D; value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0])
&& ((Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == ParentSegmentIndex[1]) || (ParentSegmentIndex[1] < 0));
} break;
default: break;
}
if (index >= 0 && value)
{
short density = Data->Density[i+LayerIndex*Data->sizeX*Data->sizeY]; DensityMeanDev[index].x += density;
DensityMeanDev[index].y += density*density; N[index]++;
}
73
}
for (size_t k = 0; k < SegmentsCount; ++k)
{
DensityMeanDev[k].x /= N[k];
DensityMeanDev[k].y = sqrt(DensityMeanDev[k].y/N[k]- DensityMeanDev[k].x*DensityMeanDev[k].x);
}
delete [] N;
}
Листинг B.6: Средние значения и среднеквадратические отклонения координат X и Y
Vec4f MeanDevXYOfSegment_2D(TVoxelsData* Data, size_t LayerIndex, Vec2i ParentSegmentIndex, int SegmentIndex, int SegmentType)
{
float sqr_x = 0.0f, sqr_y = 0.0f, mean_x = 0.0f, mean_y = 0.0f; int N = 0;
for (size_t i = 0; i < Data->sizeX; ++i) for (size_t j = 0; j < Data->sizeY; ++j)
{
bool value;
switch (SegmentType)
{
case DENSITY:
{
value = Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == SegmentIndex;
}
break;
case SPATIAL:
{
value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) && (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == SegmentIndex);
} break;
case CONNECTED:
{
value = (Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SegmentIndex_2D == ParentSegmentIndex[0]) && ((Data->VoxelSegments[i+ LayerIndex*Data->sizeX*Data->sizeY].SpatialSegmentIndex_2D == ParentSegmentIndex[1]) || (ParentSegmentIndex[1] < 0)) && (Data->VoxelSegments[i+LayerIndex*Data->sizeX*Data->sizeY].ComponentIndex_2D
== SegmentIndex);
}break;
default: break;
}
if (value)
{
sqr_x += i*i; sqr_y += j*j; mean_x += i; mean_y += j; N++;
}
}
74
return Vec4f(mean_x/N, mean_y/N,
sqrt(sqr_x/N-(mean_x/N)*(mean_x/N)), sqrt(sqr_y/N-(mean_y/N)*(mean_y/N)));
}
void CalcMeanDevXYOfSegments_2D(TVoxelsData* Data, size_t LayerIndex, Vec2i ParentSegmentIndex, int SegmentType, size_t SegmentsCount, vector <Vec4f> &MeanDevXY)
{
int * N = new int [SegmentsCount]; MeanDevXY.clear();
for (size_t k = 0; k < SegmentsCount; ++k)
{
MeanDevXY.push_back(Vec4f(0.0f, 0.0f, 0.0f, 0.0f)); N[k] = 0;
}
for (size_t i = 0; i < Data->sizeX; ++i) for (size_t j = 0; j < Data->sizeY; ++j)
{
int index; bool value;
switch (SegmentType)
{
case DENSITY:
{
index = Data->VoxelSegments[Data->ReducedIndex(i, j, LayerIndex)].SegmentIndex_2D;
value = true; } break;
case SPATIAL:
{
index = Data->VoxelSegments[Data->ReducedIndex(i, j, LayerIndex)].SpatialSegmentIndex_2D;
value = Data->VoxelSegments[Data->ReducedIndex(i, j, LayerIndex)].SegmentIndex_2D == ParentSegmentIndex[0];
} break;
case CONNECTED:
{
index = Data->VoxelSegments[Data->ReducedIndex(i, j, LayerIndex)].ComponentIndex_2D;
value = (Data->VoxelSegments[Data->ReducedIndex(i, j, LayerIndex)].SegmentIndex_2D == ParentSegmentIndex[0]) && ((Data->VoxelSegments[Data->ReducedIndex(i, j, LayerIndex)].SpatialSegmentIndex_2D == ParentSegmentIndex[1]) || (ParentSegmentIndex[1] < 0));
} break;
default: break;
}
if (index >= 0 && value)
{
MeanDevXY[index][0] += i; MeanDevXY[index][1] += j; MeanDevXY[index][2] += i*i; MeanDevXY[index][3] += j*j; N[index]++;
}
}
for (size_t k = 0; k < SegmentsCount; ++k)
{
MeanDevXY[k][0] /= N[k]; MeanDevXY[k][1] /= N[k];
75