// neighborhood
std::array<int, 6> elment6 = { 4,10,12,14,16,22 };
std::array<int, 18> elment18 = { 1,3,4,5,7,9,10,11,12,14,15,16,17,19,21,22,23,25 };
std::array<int, 18> elment6_18 = { 1,3,5,7,9,11,15,17,19,21,23,25 };
std::array<int, 18> elment18_26 = { 0,2,6,8,18,20,24,26 };
itk::Neighborhood<float, 3> neighborhood6, neighborhood18, neighborhood26;
neighborhood6.SetRadius(1);
neighborhood18.SetRadius(1);
neighborhood26.SetRadius(1);
for (unsigned int i = 0; i < neighborhood26.GetSize()[0] * neighborhood26.GetSize()[1] * neighborhood26.GetSize()[2]; ++i)
{
neighborhood18[i] = 0;
neighborhood6[i] = 0;
neighborhood26[i] = 1;
}
neighborhood26[13] = 0;
for (auto index : elment6)
{
neighborhood6[index] = 1;
}
for (auto index : elment18)
{
neighborhood18[index] = 1;
}
itk::Neighborhood<float, 3> neighborhood_6_18, neighborhood_18_26;
neighborhood_6_18.SetRadius(1);
neighborhood_18_26.SetRadius(1);
for (unsigned int i = 0; i < neighborhood26.GetSize()[0] * neighborhood26.GetSize()[1] * neighborhood26.GetSize()[2]; ++i)
{
neighborhood_6_18[i] = neighborhood18[i] - neighborhood6[i];
neighborhood_18_26[i] = neighborhood26[i] - neighborhood18[i];
}
using NeighborhoodOperatorImageFunctionType = itk::NeighborhoodOperatorImageFunction<InputImageType, float>;
auto neighborhoodOperatorImageFunction6 = NeighborhoodOperatorImageFunctionType::New();
neighborhoodOperatorImageFunction6->SetOperator(neighborhood6);
neighborhoodOperatorImageFunction6->SetInputImage(fix);
auto neighborhoodOperatorImageFunction18 = NeighborhoodOperatorImageFunctionType::New();
neighborhoodOperatorImageFunction18->SetOperator(neighborhood18);
neighborhoodOperatorImageFunction18->SetInputImage(fix);
auto neighborhoodOperatorImageFunction6_18 = NeighborhoodOperatorImageFunctionType::New();
neighborhoodOperatorImageFunction6_18->SetOperator(neighborhood_6_18);
neighborhoodOperatorImageFunction6_18->SetInputImage(fix);
auto neighborhoodOperatorImageFunction18_26 = NeighborhoodOperatorImageFunctionType::New();
neighborhoodOperatorImageFunction18_26->SetOperator(neighborhood_18_26);
neighborhoodOperatorImageFunction18_26->SetInputImage(fix);
typedef itk::ImageRegionIterator<InputImageType> FieldIterator;
FieldIterator fieldIterFix(fix, fix->GetLargestPossibleRegion());
fieldIterFix.GoToBegin();
InputImageType::ValueType value;
itk::Index<3> centerIndex;
int pointId = 0;
itk::Point<float, 3> temp;
while (!fieldIterFix.IsAtEnd()) {
value = fieldIterFix.Get();
if (value > 0) {
centerIndex = fieldIterFix.GetIndex();
temp[0] = centerIndex[0] * spacefix[0] + orifix[0];
temp[1] = centerIndex[1] * spacefix[1] + orifix[1];
temp[2] = centerIndex[2] * spacefix[2] + orifix[2];
fixedPointContainer->InsertElement(pointId, temp);
fixedPoints->InsertPoint(pointId, temp[0], temp[1], temp[2]);
pointId++;
}
++fieldIterFix;
}
pointId = 0;
while (!fieldIterMov.IsAtEnd()) {
value = fieldIterMov.Get();
if (value > 0) {
centerIndex = fieldIterMov.GetIndex();
temp[0] = centerIndex[0] * spacemov[0] + orimov[0];
temp[1] = centerIndex[1] * spacemov[1] + orimov[1];
temp[2] = centerIndex[2] * spacemov[2] + orimov[2];
movingPointContainer->InsertElement(pointId, temp);
movingPoints->InsertPoint(pointId, temp[0], temp[1], temp[2]);
pointId++;
}
++fieldIterMov;
}
pointId = 0;
fieldIterFix.GoToBegin();
/*
while (!fieldIterFix.IsAtEnd())
{
int number = 0;
value = fieldIterFix.Get();
if (value > 0) {
centerIndex = fieldIterFix.GetIndex();
std::vector< itk::Index<3>> p_6_18, p_18_26;
if (neighborhoodOperatorImageFunction18_26->EvaluateAtIndex(centerIndex) > 0)
{
for (auto i_18_26 : elment18_26)
{
p_18_26.push_back(centerIndex + neighborhood_18_26.GetOffset(i_18_26));
}
number++;
}
if (neighborhoodOperatorImageFunction6_18->EvaluateAtIndex(centerIndex) > 0)
{
for (auto i_6_18 : elment6_18)
{
p_6_18.push_back(centerIndex + neighborhood_6_18.GetOffset(i_6_18));
}
for (auto elment : p_6_18)
{
if (neighborhoodOperatorImageFunction6->EvaluateAtIndex(elment) == 0)
{
number++;
break;
}
else {
std::vector< itk::Index<3>> x_6;
for (auto i_6 : elment6)
{
x_6.push_back(elment + neighborhood6.GetOffset(i_6));
}
if (!test(x_6, p_18_26))
{
number++;
break;
}
}
}
}
if (neighborhoodOperatorImageFunction6->EvaluateAtIndex(centerIndex) > 0)
{
for (auto i_6 : elment6)
{
auto x = centerIndex + neighborhood6.GetOffset(i_6);
if (neighborhoodOperatorImageFunction18->EvaluateAtIndex(x) > 0)
{
std::vector< itk::Index<3>> x_18;
// x 18
for (auto i_18 : elment18)
{
x_18.push_back(x + neighborhood18.GetOffset(i_18));
}
if (!test(x_18, p_6_18) && !test(x_18, p_18_26))
{
number++;
break;
}
}
else
{
number++;
break;
}
}
}
if (number > 1)
{
fieldIterFix.Set(1);
temp[0] = centerIndex[0] * spacefix[0] + orifix[0];
temp[1] = centerIndex[1] * spacefix[1] + orifix[1];
temp[2] = centerIndex[2] * spacefix[2] + orifix[2];
fixedPointContainer->InsertElement(pointId, temp);
pointId++;
}
else
fieldIterFix.Set(0);
}
++fieldIterFix;
number++;
}