我正在使用 ITK 进行一些预处理,并且想使用快速行进滤波器和测地线事件轮廓滤波器进行测试。 我遵循 ITK software guide 中描述的算法,第 9.3.3 节。 但是,我没有得到预期的结果。我正在处理 3D 图像。
这是我的代码:
AnisotropicDiffusionFilter::Pointer anisotropic_filter = AnisotropicDiffusionFilter::New();
anisotropic_filter->SetInput(itk_image_in);
anisotropic_filter->SetTimeStep(0.0625);
anisotropic_filter->SetNumberOfIterations(5);
anisotropic_filter->SetConductanceParameter(3.0);
anisotropic_filter->Update();
GradientFilter::Pointer gradient_filter = GradientFilter::New();
gradient_filter->SetInput(anisotropic_filter->GetOutput());
gradient_filter->SetSigma(0.5);
gradient_filter->Update();
SigmoidFilter::Pointer sigmoid_filter = SigmoidFilter::New();
sigmoid_filter->SetInput(gradient_filter->GetOutput());
sigmoid_filter->SetOutputMinimum(0.0);
sigmoid_filter->SetOutputMaximum(1.0);
sigmoid_filter->SetAlpha(-1.5);
sigmoid_filter->SetBeta(4.0);
sigmoid_filter->Update();
FastMarchingFilter::Pointer fast_marching = FastMarchingFilter::New();
NodeContainer::Pointer seeds = NodeContainer::New();
Node node;
const double seedValue = -50.0;
node.SetValue(seedValue);
seeds->Initialize();
vector<GeoVec3s>::iterator it = m_clicks_.begin();
int i=0;
for(; it != m_clicks_.end(); it++)
{
itkIndex index;
index[0] = (*it)[0];
index[1] = (*it)[1];
index[2] = (*it)[2];
node.SetIndex(index);
seeds->InsertElement(i++, node);
}
fast_marching->SetTrialPoints(seeds);
fast_marching->SetSpeedConstant(1.0);
fast_marching->SetStoppingValue(100);
//fast_marching->SetInput(sigmoid_filter->GetOutput());
fast_marching->SetOutputSize(sigmoid_filter->GetOutput()->GetBufferedRegion().GetSize());
fast_marching->Update();
GeodesicFilter::Pointer geodesic_filter = GeodesicFilter::New();
geodesic_filter->SetInput(fast_marching->GetOutput());
geodesic_filter->SetFeatureImage(sigmoid_filter->GetOutput());
geodesic_filter->SetPropagationScaling(0.5);
geodesic_filter->SetCurvatureScaling(5.0);
geodesic_filter->SetAdvectionScaling(1.0);
geodesic_filter->SetMaximumRMSError( 0.02 );
geodesic_filter->Update();
BinaryThresholdFilter::Pointer thresholder = BinaryThresholdFilter::New();
thresholder->SetLowerThreshold(-1000);
thresholder->SetUpperThreshold(0);
thresholder->SetOutsideValue(0);
thresholder->SetInsideValue(255);
thresholder->SetInput( geodesic_filter->GetOutput() );
我正在使用 this paper 中描述的指标哪个目标和我的一样。
我有几个问题:
- 快速行进过滤器应该输出距离图,对吗?相反,当我将音量输出为一系列 png(值在 0 到 4095 之间)时,我有一个二进制图像(像素为 0 或 4095)。我想我应该得到一个灰度体积,指示从种子获得每个像素所需的时间。
- 按照铃木描述的过程,我成功地使算法或多或少地工作,但是我改变了测地线滤波器的参数值。我不记得确切的值,但它与论文中描述的值并不接近。当我们使用在 0 和 1 之间标准化的 sigmoid 输入时,发生了什么?
- 我应该对快速行进滤波器或 sigmoid 图像使用恒速函数吗?什么时候应该首选任一方法?
- 我正在使用重新缩放器来输出 float 图像(从过滤器输出)。这可能是我看到的不一致的原因吗?
- 对我可能做错的事情有什么建议吗?
谢谢。
最佳答案
好的,我发现了我的问题。快速行进过滤器确实输出时间交叉图(距离图),但当我在算法中指定停止值时,所有未访问的像素都具有较高的值(1.7e+38,因为它是最大值的一半)用于输出图像的类型,在我的例子中是 float 的 3.4e+38)。因此,当我使用重新缩放滤镜时,它会压缩所有图像动态,结果是二进制图像。 我认为使用 sigmoid 图像作为快速行进滤波器的输入可以获得更好的结果。 感谢@nav 的建议。
关于c++ - ITK快速行进输出,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24488370/