Itk:: neighborhood get 6 neighborhood, 18 neighborhood, 26 neighborhood, 18/6 neighborhood, 26/18 neighborhood

   // 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++;

    }