matlab上3D数据的椭球拟合 [英] Ellipsoid fitting for 3D data on matlab
问题描述
我正在处理 CT 肺部图像的 3D 体积，为了检测结节，我需要为每个疑似结节拟合一个椭球模型，我该如何制作一个代码???结节是疑似肿瘤的对象，我的算法需要检查每个对象，并将其近似为一个椭球，我们从椭球参数中计算出8个特征来构建一个分类器，通过训练和测试来检测它是否是一个结节数据，所以我需要拟合这样的椭球
这里是一张体积CT肺图像
这里是另一个相同体积的切片，但它包含一个结节(黄色圆圈中有一个结节)所以我需要我的代码来检查每个形状以确定它是否是一个结节
由于我们没有 3D 数据集可供使用，我从 2D 开始.
所以首先我们需要选择肺，这样我们就不会计算内部的任何其他对象.由于这是灰度，我们首先需要以某种方式将其二值化.我使用我自己的类图片作为DIP，这将大量使用我的growthfill，所以我强烈建议先阅读:
二值化图像
只需根据直方图中的**绿色*强度对图像进行阈值处理.因此，如果
i0,i1是直方图中左右两侧的局部最大强度，则针对(i0+i1)/2的阈值.结果如下:除肺外的所有东西都去除
这很简单，只需将黑色从外部填充到某些预定义的背景颜色即可.然后以同样的方式将所有白色的东西相邻的背景颜色.这将去除人体表面、骨骼、器官和 CT 机器，只留下肺部.现在用一些预定义的肺颜色重新着色剩余的黑色.
不应该留下黑色，剩下的白色是可能的结节.
处理所有剩余的白色像素
因此只需循环遍历图像，并在第一个白色像素命中时使用预定义的结节颜色或不同的对象索引填充它以供以后使用.我还区分了表面(浅绿色)和内部(洋红色).结果如下:
现在您可以计算每个结节的特征.如果您为此编写自定义
floodfill代码，那么您可以直接从中获取以下内容:- 体积
[像素] [pixels] 中的表面积- 边界框
- 位置(相对于肺)
- 重心或质心
所有这些都可以用作您的特征变量并有助于拟合.
- 体积
拟合找到的表面点
有很多方法可以做到这一点，但我会尽可能地简化它以提高性能和准确性.例如，您可以使用质心作为椭圆体中心.然后找到距离它的
min和max远点并将它们用作半轴(+/- 一些正交校正).然后只搜索这些初始值.有关更多信息，请参阅:您会在链接的问答中找到使用示例.
[注释]
所有项目符号都适用于 3D.在构建自定义 floodfill 时要注意递归尾.太多的信息会很快溢出你的堆栈，也会大大减慢速度.这是我如何使用我使用的几个自定义返回参数 + growthfill 处理它的小例子:
//---------------------------------------------------------------------无效增长填充(双字c0，双字c1，双字c2)；//用 c2 生长/淹没 c0 邻近 c1void floodfill(int x,int y,DWORD c);//从 (x,y) 填充颜色 cDWORD _floodfill_c0,_floodfill_c1;//递归填充颜色和填充颜色int _floodfill_x0,_floodfill_x1,_floodfill_n;//递归边界框和填充像素数int _floodfill_y0,_floodfill_y1;void _floodfill(int x,int y);//洪水填充的递归//---------------------------------------------------------------------------无效图片::growfill(DWORD c0,DWORD c1,DWORD c2){整数 x,y,e;对于 (e=1;e;)对于 (e=0,y=1;yx)_floodfill_x0=x；如果(_floodfill_y0>y)_floodfill_y0=y；如果 (_floodfill_x1 0) _floodfill(x-1,y);如果 (x 0)_floodfill(x，y-1)；如果 (y=xs)||(y<0)||(y>=ys)) 返回；_floodfill_c0=p[y][x].dd;_floodfill_c1=c;_floodfill_n=0;_floodfill_x0=x;_floodfill_y0=y;_floodfill_x1=x;_floodfill_y1=y;_floodfill(x,y);}//--------------------------------------------------------------------------- 这里的 C++ 代码是我使用的示例图像:
picture pic0,pic1;//pic0 - 源 img//pic1 - 输出 imgint x0,y0,x1,y1,x,y,i,hist[766];颜色 c;//复制源图像到输出pic1=pic0;pic1.pixel_format(_pf_u);//灰度 <0,765>//0xAARRGGBBconst DWORD col_backg=0x00202020;//灰色的const DWORD col_lungs=0x00000040;//蓝色const DWORD col_out =0x0000FFFF;//水色结节表面const DWORD col_in =0x00800080;//里面有品红色结节const DWORD col_test =0x00008040;//绿色的独特颜色只是为了安全重新着色//[去除白色背景]//找到白色背景区域(通过将光线从边框中间投射到图像中心直到命中非白色像素)对于 (x0=0 ,y=pic1.ys>>1;x01;x1>0;x1--) 如果 (pic1.p[y][x1].dd<700) 中断;for (y0=0 ,x=pic1.xs>>1;y0<pic1.ys;y0++) if (pic1.p[y0][x].dd<700);对于 (y1=pic1.ys-1,x=pic1.xs>1;y1>0;y1--) 如果 (pic1.p[y1][x].dd<700) 中断;//剪掉它pic1.bmp->Canvas->Draw(-x0,-y0,pic1.bmp);pic1.resize(x1-x0+1,y1-y0+1);//[准备数据]//原始直方图对于 (i=0;i<766;i++) hist[i]=0;对于 (y=0;y>1;对于 (i=765;i>0;i--) hist[i]=(hist[i]+hist[i-1])>>1;}//在直方图中找到峰值(用于阈值)对于 (x=0,x0=x,y0=hist[x];x<766;x++){y=历史[x];如果 (y0=0;x--){y=历史[x];如果 (y1=i) pic1.p[y][x].dd=0x00FFFFFF;否则 pic1.p[y][x].dd=0x00000000;pic1.save("out0.png");//重新着色外部背景对于 (x=0;x>1>6)&&(y>1].dd=0x00FF0000；对于 (x=x1,y=0;(y<=100)&&(y>1].dd=0x00FF0000；对于 (x=(x0+x1)>>1,y=0;(y<=100)&&(y I am working on a 3D volume of CT lung images, In order to detect nodules I need to fit an ellipsoid model for each suspected nodule, How can I make a code for that ??? Nodule is the suspected object to be a tumor, my algorithm needs to check every object, and approximate it to an ellipsoid, and from the ellipsoid parameters we calculate 8 features to build a classifier which detects whether it a nodule or not through training and testing data, so I need to fit such ellipsoid
here one slice of the volume CT lung image
here another slice of the same volume but it contains a nodule (the yellow circle there is a nodule) so I need my code to check every shape to determine whether it is a nodule or not
As we do not have 3D dataset at disposal I start with 2D.
So first we need to select the lungs so we do not count any other objects then those inside. As this is gray scale we first need to binarize it somehow. I use my own class picture for DIP and this will heavily use my growthfill so I strongly recommend to first read this:
Where you will find all explanations you need for this. Now for your task I would:
turn RGBA to grayscale
<0,765>I just compute intensity
i=R+G+Bas for 24 bit image the channels are 8 bit the result is up to3*255 = 765. As the input image was compressed by JPEG there are distortions in color and also noise in the image so do not forget about that.crop out the white border
Just cast rays (scan lines) from middle of each outer line of the border towards middle and stop if non white-ish pixel hit. I treshold with
700instead of765to compensate the noise in image. Now you got the bounding box of usable image so crop out the rest.compute histogram
To compensate distortions in image smooth the histogram enough to remove all unwanted noise and gaps. Then find the local maximum from left and from right (red). This will be used for binarisation treshold (middle between them green) This is my final result:
binarise image
Just treshold image against the **green* intensity from histogram. So if
i0,i1are the local maximum intensities from left and right in the histogram then treshold against(i0+i1)/2. This is the result:remove everything except lungs
That is easy just fill the black from outside to some predefined background color. Then the same way all the white stuff neighboring background color. That will remove human surface, skeleton, organs and the CT machine leaving just the lungs. Now recolor the remaining black with some predefined Lungs color.
There should be no black color left and the remaining white are the possible nodules.
process all remaining white pixels
So just loop through the image and on first white pixel hit flood fill it with predefined nodule color or distinct object index for latter use. I also distinct the surface (aqua) and the inside (magenta). This is the result:
Now you can compute your features per nodule. If you code your custom
floodfillfor this then you can obtain from it directly things like:- Volume in
[pixels] - Surface area in
[pixels] - Bounding box
- Position (relative to Lungs)
- Center of gravity or centroid
Which all can be used as your feature variables and also to help with fitting.
- Volume in
fit the found surface points
There are many methods for this but I would ease up it as much as I could to improve performance and accuracy. For example you can use centroid as your ellipsoid center. Then find the
minandmaxdistant points from it and use them as semi-axises (+/- some orthogonality corrections). And then just search around these initial values. For more info see:You will find examples of use in linked Q/As there.
[Notes]
All of the bullets are applicable in 3D. While constructing custom floodfill be careful with the recursion tail. Too much info will really fast overflow your stack and also slow things down considerably. Here small example of how I deal with it with few custom return parameters + growthfill I used:
//---------------------------------------------------------------------------
void growfill(DWORD c0,DWORD c1,DWORD c2); // grow/flood fill c0 neigbouring c1 with c2
void floodfill(int x,int y,DWORD c); // flood fill from (x,y) with color c
DWORD _floodfill_c0,_floodfill_c1; // recursion filled color and fill color
int _floodfill_x0,_floodfill_x1,_floodfill_n; // recursion bounding box and filled pixel count
int _floodfill_y0,_floodfill_y1;
void _floodfill(int x,int y); // recursion for floodfill
//---------------------------------------------------------------------------
void picture::growfill(DWORD c0,DWORD c1,DWORD c2)
{
int x,y,e;
for (e=1;e;)
for (e=0,y=1;y<ys-1;y++)
for ( x=1;x<xs-1;x++)
if (p[y][x].dd==c0)
if ((p[y-1][x].dd==c1)
||(p[y+1][x].dd==c1)
||(p[y][x-1].dd==c1)
||(p[y][x+1].dd==c1)) { e=1; p[y][x].dd=c2; }
}
//---------------------------------------------------------------------------
void picture::_floodfill(int x,int y)
{
if (p[y][x].dd!=_floodfill_c0) return;
p[y][x].dd=_floodfill_c1;
_floodfill_n++;
if (_floodfill_x0>x) _floodfill_x0=x;
if (_floodfill_y0>y) _floodfill_y0=y;
if (_floodfill_x1<x) _floodfill_x1=x;
if (_floodfill_y1<y) _floodfill_y1=y;
if (x> 0) _floodfill(x-1,y);
if (x<xs-1) _floodfill(x+1,y);
if (y> 0) _floodfill(x,y-1);
if (y<ys-1) _floodfill(x,y+1);
}
void picture::floodfill(int x,int y,DWORD c)
{
if ((x<0)||(x>=xs)||(y<0)||(y>=ys)) return;
_floodfill_c0=p[y][x].dd;
_floodfill_c1=c;
_floodfill_n=0;
_floodfill_x0=x;
_floodfill_y0=y;
_floodfill_x1=x;
_floodfill_y1=y;
_floodfill(x,y);
}
//---------------------------------------------------------------------------
And here C++ code I did the example images with:
picture pic0,pic1;
// pic0 - source img
// pic1 - output img
int x0,y0,x1,y1,x,y,i,hist[766];
color c;
// copy source image to output
pic1=pic0;
pic1.pixel_format(_pf_u); // grayscale <0,765>
// 0xAARRGGBB
const DWORD col_backg=0x00202020; // gray
const DWORD col_lungs=0x00000040; // blue
const DWORD col_out =0x0000FFFF; // aqua nodule surface
const DWORD col_in =0x00800080; // magenta nodule inside
const DWORD col_test =0x00008040; // green-ish distinct color just for safe recoloring
// [remove white background]
// find white background area (by casting rays from middle of border into center of image till non white pixel hit)
for (x0=0 ,y=pic1.ys>>1;x0<pic1.xs;x0++) if (pic1.p[y][x0].dd<700) break;
for (x1=pic1.xs-1,y=pic1.ys>>1;x1> 0;x1--) if (pic1.p[y][x1].dd<700) break;
for (y0=0 ,x=pic1.xs>>1;y0<pic1.ys;y0++) if (pic1.p[y0][x].dd<700) break;
for (y1=pic1.ys-1,x=pic1.xs>>1;y1> 0;y1--) if (pic1.p[y1][x].dd<700) break;
// crop it away
pic1.bmp->Canvas->Draw(-x0,-y0,pic1.bmp);
pic1.resize(x1-x0+1,y1-y0+1);
// [prepare data]
// raw histogram
for (i=0;i<766;i++) hist[i]=0;
for (y=0;y<pic1.ys;y++)
for (x=0;x<pic1.xs;x++)
hist[pic1.p[y][x].dd]++;
// smooth histogram a lot (remove noise and fill gaps due to compression and scanning nature of the image)
for (x=0;x<100;x++)
{
for (i=0;i<765;i++) hist[i]=(hist[i]+hist[i+1])>>1;
for (i=765;i>0;i--) hist[i]=(hist[i]+hist[i-1])>>1;
}
// find peaks in histogram (for tresholding)
for (x=0,x0=x,y0=hist[x];x<766;x++)
{
y=hist[x];
if (y0<y) { x0=x; y0=y; }
if (y<y0>>1) break;
}
for (x=765,x1=x,y1=hist[x];x>=0;x--)
{
y=hist[x];
if (y1<y) { x1=x; y1=y; }
if (y<y1>>1) break;
}
// binarize image (tresholding)
i=(x0+x1)>>1; // treshold with middle intensity between peeks
pic1.pf=_pf_rgba; // result will be RGBA
for (y=0;y<pic1.ys;y++)
for (x=0;x<pic1.xs;x++)
if (pic1.p[y][x].dd>=i) pic1.p[y][x].dd=0x00FFFFFF;
else pic1.p[y][x].dd=0x00000000;
pic1.save("out0.png");
// recolor outer background
for (x=0;x<pic1.xs;x++) pic1.p[ 0][x].dd=col_backg; // render rectangle along outer border so the filling starts from there
for (x=0;x<pic1.xs;x++) pic1.p[pic1.ys-1][x].dd=col_backg;
for (y=0;y<pic1.ys;y++) pic1.p[y][ 0].dd=col_backg;
for (y=0;y<pic1.ys;y++) pic1.p[y][pic1.xs-1].dd=col_backg;
pic1.growfill(0x00000000,col_backg,col_backg); // fill its black content outside in
// recolor white human surface and CT machine
pic1.growfill(0x00FFFFFF,col_backg,col_backg);
// recolor Lungs dark matter
pic1.growfill(0x00000000,col_backg,col_lungs); // outer border
pic1.growfill(0x00000000,col_lungs,col_lungs); // fill its black content outside in
pic1.save("out1.png");
// find/recolor individual nodules
for (y=0;y<pic1.ys;y++)
for (x=0;x<pic1.xs;x++)
if (pic1.p[y][x].dd==0x00FFFFFF)
{
pic1.floodfill(x,y,col_test);
pic1.growfill(col_lungs,col_test,col_out);
pic1.floodfill(x,y,col_in);
}
pic1.save("out2.png");
// render histogram
for (x=0;(x<766)&&(x>>1<pic1.xs);x++)
for (y=0;(y<=hist[x]>>6)&&(y<pic1.ys);y++)
pic1.p[pic1.ys-1-y][x>>1].dd=0x000040FF;
for (x=x0 ,y=0;(y<=100)&&(y<pic1.ys);y++) pic1.p[pic1.ys-1-y][x>>1].dd=0x00FF0000;
for (x=x1 ,y=0;(y<=100)&&(y<pic1.ys);y++) pic1.p[pic1.ys-1-y][x>>1].dd=0x00FF0000;
for (x=(x0+x1)>>1,y=0;(y<=100)&&(y<pic1.ys);y++) pic1.p[pic1.ys-1-y][x>>1].dd=0x0000FF00;
这篇关于matlab上3D数据的椭球拟合的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
