I am trying to implement convex hull algorithm calling hit-or-miss several times with 90 degrees rotated structural elements.
It is not working properly.
using ShapedNeighborhoodIteratorType_Heart_Hull =itk::ConstShapedNeighborhoodIterator<InternalImageType>;
using IteratorType_Heart_Hull = itk::ImageRegionIterator<InternalImageType>;
using FaceCalculatorType_Heart_Hull =
itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InternalImageType>;
using MaximumImageFilterType_Heart_Hull = itk::MaximumImageFilter<InternalImageType>;
unsigned int element_radius_heart_hull = 1;
ShapedNeighborhoodIteratorType_Heart_Hull::RadiusType radius_heart_hull;
radius_heart_hull.Fill(element_radius_heart_hull);
// const auto rad_convex_hull = static_cast<float>(element_radius_heart_hull);
////
InternalImageType::Pointer x0 = InternalImageType::New();
InternalImageType::Pointer x1 = InternalImageType::New();
InternalImageType::Pointer x2 = InternalImageType::New();
InternalImageType::Pointer x3 = InternalImageType::New();
InternalImageType::Pointer x4 = InternalImageType::New();
InternalImageType::Pointer x1c = InternalImageType::New();
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x0);
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x1);
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x2);
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x3);
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x4);
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x1c);
x1->FillBuffer(0.0);
x2->FillBuffer(0.0);
x3->FillBuffer(0.0);
x4->FillBuffer(0.0);
x1c->FillBuffer(0.0);
x0 -> SetSpacing(spacing);
x0 -> SetOrigin(newOrigin);
x0 -> SetDirection(direction);
x1 -> SetSpacing(spacing);
x1 -> SetOrigin(newOrigin);
x1 -> SetDirection(direction);
x2 -> SetSpacing(spacing);
x2 -> SetOrigin(newOrigin);
x2 -> SetDirection(direction);
x3 -> SetSpacing(spacing);
x3 -> SetOrigin(newOrigin);
x3 -> SetDirection(direction);
x4 -> SetSpacing(spacing);
x4 -> SetOrigin(newOrigin);
x4 -> SetDirection(direction);
x1c -> SetSpacing(spacing);
x1c -> SetOrigin(newOrigin);
x1c -> SetDirection(direction);
bool se_hull_ready = false;
for(int convex_hull_counter = 1; convex_hull_counter < 5 ; convex_hull_counter++)
{
se_hull_ready = false;
DeepCopy<InternalImageType>(Segmented_Lungs_Heart_for_close_float, x0);
int while_counter_convex_hull = 0;
while(se_hull_ready== false){
IteratorType_Heart_Hull out_heart_hull;
IteratorType_Heart_Hull out_heart_hull_2;
constexpr InternalPixelType background_value = 0.0;
constexpr InternalPixelType foreground_value = 255.0;
FaceCalculatorType_Heart_Hull faceCalculator_heart_hull;
FaceCalculatorType_Heart_Hull::FaceListType faceList_heart_hull;
FaceCalculatorType_Heart_Hull faceCalculator_heart_hull_c;
FaceCalculatorType_Heart_Hull::FaceListType faceList_heart_hull_c;
auto invertIntensityFilter_heart_hull = InvertIntensityImageFilterType_Heart_btw::New();
invertIntensityFilter_heart_hull->SetInput(x0);
invertIntensityFilter_heart_hull->SetMaximum(255);
invertIntensityFilter_heart_hull->Update();
auto x01c = invertIntensityFilter_heart_hull->GetOutput();
x1 -> FillBuffer(0.0);
x1c -> FillBuffer(0.0);
faceList_heart_hull = faceCalculator_heart_hull(x0,
x1->GetRequestedRegion(),
radius_heart_hull);
faceList_heart_hull_c = faceCalculator_heart_hull_c(x01c,
x1c->GetRequestedRegion(),
radius_heart_hull);
FaceCalculatorType_Heart_Hull::FaceListType::iterator fit_heart_hull;
FaceCalculatorType_Heart_Hull::FaceListType::iterator fit_heart_hull_2;
std::cerr << "define ok " << std::endl;
for( fit_heart_hull = faceList_heart_hull.begin(); fit_heart_hull != faceList_heart_hull.end(); ++fit_heart_hull)
{
// set iterators
ShapedNeighborhoodIteratorType_Heart_Hull shape_it(radius_heart_hull, x0, *fit_heart_hull);
out_heart_hull = IteratorType_Heart_Hull( x1, *fit_heart_hull );
std::cerr << "define it " << std::endl;
float half = element_radius_heart_hull/2;
for (float y = -element_radius_heart_hull; y <= element_radius_heart_hull ; y++)
{
for (float x = -element_radius_heart_hull; x <= element_radius_heart_hull ; x++)
{
for (float z = -element_radius_heart_hull; z <= element_radius_heart_hull ; z++)
{
ShapedNeighborhoodIteratorType_Heart_Hull::OffsetType off;
if(convex_hull_counter == 1){
if((x <= 0.0)&&(z = 0.0)){
off[0] = static_cast<int>(255);
off[1] = static_cast<int>(255);
off[2] = static_cast<int>(255);
shape_it.ActivateOffset(off);}
if(((x == 0.0)&&(abs(y)<half))&&(z = 0.0)){
off[0] = static_cast<int>(0);
off[1] = static_cast<int>(0);
off[2] = static_cast<int>(0);
shape_it.ActivateOffset(off);
}
}
else if(convex_hull_counter == 2){
if((x >= 0.0)&&(z = 0.0)){
off[0] = static_cast<int>(255);
off[1] = static_cast<int>(255);
off[2] = static_cast<int>(255);
shape_it.ActivateOffset(off);}
if(((x == 0.0)&&(abs(y)<half))&&(z = 0.0)){
off[0] = static_cast<int>(0);
off[1] = static_cast<int>(0);
off[2] = static_cast<int>(0);
shape_it.ActivateOffset(off);
}
}
else if(convex_hull_counter == 3){
if((y <= 0.0)&&(z = 0.0)){
off[0] = static_cast<int>(255);
off[1] = static_cast<int>(255);
off[2] = static_cast<int>(255);
shape_it.ActivateOffset(off);}
if(((y == 0.0)&&(abs(x)<half))&&(z = 0.0)){
off[0] = static_cast<int>(0);
off[1] = static_cast<int>(0);
off[2] = static_cast<int>(0);
shape_it.ActivateOffset(off);
}
}
else if(convex_hull_counter == 4){
if((y >= 0.0)&&(z = 0.0)){
off[0] = static_cast<int>(255);
off[1] = static_cast<int>(255);
off[2] = static_cast<int>(255);
shape_it.ActivateOffset(off);}
if(((y == 0.0)&&(abs(x)<half))&&(z = 0.0)){
off[0] = static_cast<int>(0);
off[1] = static_cast<int>(0);
off[2] = static_cast<int>(0);
shape_it.ActivateOffset(off);
}
}
}
}
}
// Implements erosion
for (shape_it.GoToBegin(), out_heart_hull.GoToBegin(); !shape_it.IsAtEnd(); ++shape_it, ++out_heart_hull){
ShapedNeighborhoodIteratorType_Heart_Hull::ConstIterator ci;
bool flag = true;
for (ci = shape_it.Begin(); ci != shape_it.End(); ci++)
{
if (ci.Get() == background_value)
{
flag = false;
break;
}
}
if (flag == true)
{
out_heart_hull.Set(foreground_value);
}
else
{
out_heart_hull.Set(background_value);
}
}
}
std::cerr << "erosion ok " << std::endl;
// COMPLEMENT
for ( fit_heart_hull_2 = faceList_heart_hull_c.begin(); fit_heart_hull_2 != faceList_heart_hull_c.end(); ++fit_heart_hull_2)
{
// set iterators
ShapedNeighborhoodIteratorType_Heart_Hull shape_it_2(radius_heart_hull, x01c, *fit_heart_hull_2);
out_heart_hull_2 = IteratorType_Heart_Hull( x1c, *fit_heart_hull_2 );
float half = element_radius_heart_hull/2;
for (float y = -element_radius_heart_hull; y <= element_radius_heart_hull + 1; y++)
{
for (float x = -element_radius_heart_hull; x <= element_radius_heart_hull + 1; x++)
{
for (float z = -element_radius_heart_hull; z <= element_radius_heart_hull + 1; z++) {
ShapedNeighborhoodIteratorType_Heart_Hull::OffsetType off_2;
if(((y == 0.0)&&(x == 0.0))&&(z == 0)){
off_2[0] = static_cast<int>(255);
off_2[1] = static_cast<int>(255);
off_2[2] = static_cast<int>(255);
shape_it_2.ActivateOffset(off_2);}
}
}
}
// Implements erosion
for (shape_it_2.GoToBegin(), out_heart_hull_2.GoToBegin(); !shape_it_2.IsAtEnd(); ++shape_it_2, ++out_heart_hull_2){
ShapedNeighborhoodIteratorType_Heart_Hull::ConstIterator ci;
bool flag = true;
for (ci = shape_it_2.Begin(); ci != shape_it_2.End(); ci++)
{
if (ci.Get() == background_value)
{
flag = false;
break;
}
}
if (flag == true)
{
out_heart_hull_2.Set(foreground_value);
}
else
{
out_heart_hull_2.Set(background_value);
}
}
}
std::cerr << "complement ok " << std::endl;
auto maskFilter_heart_hull = MaskFilterType_Heart_btw::New();
maskFilter_heart_hull->SetInput(x1);
maskFilter_heart_hull->SetMaskImage(x1c);
maskFilter_heart_hull->Update();
auto intersection_image_hit_or_miss = maskFilter_heart_hull->GetOutput();
intersection_image_hit_or_miss -> SetSpacing(spacing);
intersection_image_hit_or_miss -> SetOrigin(newOrigin);
intersection_image_hit_or_miss -> SetDirection(direction);
auto maximumImageFilter_heart_hull = MaximumImageFilterType_Heart_Hull::New();
maximumImageFilter_heart_hull->SetInput(0, intersection_image_hit_or_miss);
maximumImageFilter_heart_hull->SetInput(1, Segmented_Lungs_Heart_for_close_float);
maximumImageFilter_heart_hull->Update();
InternalImageType::Pointer trial_image = maximumImageFilter_heart_hull->GetOutput();
InternalImageType::Pointer sum_trial_image = maximumImageFilter_heart_hull->GetOutput();
sum_trial_image->FillBuffer(0.0);
IteratorType_Heart_Hull it_trial(trial_image, trial_image->GetRequestedRegion());
IteratorType_Heart_Hull it_trial_c(x0, x0->GetRequestedRegion());
IteratorType_Heart_Hull it_trial_sum(sum_trial_image, sum_trial_image->GetRequestedRegion());
if(while_counter_convex_hull>0){
for (it_trial.GoToBegin(), it_trial_c.GoToBegin(), it_trial_sum.GoToBegin(); !it_trial.IsAtEnd(); ++it_trial, ++it_trial_c, ++it_trial_sum)
{
it_trial_sum.Set(it_trial.Get()-it_trial_c.Get());
}
int trial_counter =0;
for (it_trial_sum.GoToBegin(); !it_trial_sum.IsAtEnd(); ++it_trial_sum)
{
if(it_trial_sum.Get()!=0)trial_counter++;
}
if(trial_counter == 0){
if(convex_hull_counter==1)x1 = maximumImageFilter_heart_hull->GetOutput();
else if(convex_hull_counter==2)x2 = maximumImageFilter_heart_hull->GetOutput();
else if(convex_hull_counter==3)x3 = maximumImageFilter_heart_hull->GetOutput();
else if(convex_hull_counter==4)x4 = maximumImageFilter_heart_hull->GetOutput();
se_hull_ready = true;
}
}
while_counter_convex_hull++;
DeepCopy<InternalImageType>(trial_image, x0);
}
}
std::cerr << "iteration ok " << std::endl;
auto maximumImageFilter_heart_hull_2 = MaximumImageFilterType_Heart_Hull::New();
maximumImageFilter_heart_hull_2->SetInput(0, x1);
maximumImageFilter_heart_hull_2->SetInput(1, x2);
maximumImageFilter_heart_hull_2->Update();
auto maximumImageFilter_heart_hull_3 = MaximumImageFilterType_Heart_Hull::New();
maximumImageFilter_heart_hull_3->SetInput(0, x3);
maximumImageFilter_heart_hull_3->SetInput(1, x4);
maximumImageFilter_heart_hull_3->Update();
///
auto maximumImageFilter_heart_hull_4 = MaximumImageFilterType_Heart_Hull::New();
maximumImageFilter_heart_hull_4->SetInput(0, maximumImageFilter_heart_hull_2->GetOutput());
maximumImageFilter_heart_hull_4->SetInput(1, maximumImageFilter_heart_hull_3->GetOutput());
maximumImageFilter_heart_hull_4->Update();
auto Convex_Hull_Lungs = maximumImageFilter_heart_hull_4->GetOutput();
auto Writer_Heart_hull = WriterType_Heart_btw::New();
Writer_Heart_hull->SetFileName("Convex Hull Lungs.tiff");
Writer_Heart_hull->SetInput(Convex_Hull_Lungs);
Writer_Heart_hull->Update();