# Implementing Non-Square Boundary Conditions
**URL:** https://discourse.itk.org/t/implementing-non-square-boundary-conditions/1690
**Category:** Engineering
**Created:** [March 20, 2019, 9:48pm UTC](https://discourse.itk.org/t/implementing-non-square-boundary-conditions/1690 "2019-03-20T21:48:30Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@Besler](https://discourse.itk.org/u/Besler)
#### Post date: [March 20, 2019, 9:48pm UTC](https://discourse.itk.org/t/implementing-non-square-boundary-conditions/1690/1 "2019-03-20T21:48:30Z")
I’m wanting to implement a filter using [itkDenseFiniteDifferenceImageFilter](https://itk.org/Doxygen/html/classitk_1_1DenseFiniteDifferenceImageFilter.html) with non-square boundary conditions.
The problem is defined as:
I have many PDE problems I would like to solve with similar boundary conditions.
I only want to solve the problem inside the mask, but I want to implement boundary conditions across the mask boundary. I am OK on the finite difference approximations and other technical aspect. What I would like is some feedback on implementing the boundary condition described above.
I have read through the [itkDenseFiniteDifferenceImageFilter](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/FiniteDifference/include/itkDenseFiniteDifferenceImageFilter.hxx) code many times.
**My Thoughts on Implementation**
There are two major hurtles to overcome: the boundary condition and iterating over the domain.
**_Boundary Condition_**
For boundary condition, I was going to modify the [itkZeroFluxNeumannBoundaryCondition](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/Common/include/itkZeroFluxNeumannBoundaryCondition.hxx) to operate on an image and the mask. A quick prototype would be
`GetPixel( const IndexType & index, const TInputImage * image, const TMaskImage * mask )`
In the body, the boundary checking would be the same but with an additional check for the mask being non-zero at the index. Of course, this would assume the image and the mask have the same region, etc.
**_Iterating over the Domain_**
For iterating, the current implementation uses [ImageBoundaryFacesCalculator](https://github.com/InsightSoftwareConsortium/ITK/blob/c95b99852a9109a8fc176a06f88b2aa6d442154b/Modules/Core/Common/include/itkNeighborhoodAlgorithm.hxx) to calculate a set of image regions to iterator over. The first region is guaranteed to not overlap the boundary (given some radius) so the boundary conditions don’t need to be checked. Then, the regions which will overlap the boundary are processed in sequence. This allows a computationally fast and memory efficient filter.
I was going to implement a similar filter that returned a set of indicies (instead of regions) to be iterated over. Of course, this is crazy memory expensive. Alternatively, I could use a region iterator on the [bounding box of the mask](https://itk.org/Doxygen/html/classitk_1_1LabelGeometryImageFilter.html) and check if every index is in the boundary and then check the boundary conditions for every neighbourhood. But, this is computationally expensive.
**_Scale_**
For a sense of scale, the size of datasets I am working with are in the range of 700x700x168 and 512x512x512 voxels. You can assume the masks defining the domain are relativly blob like (low surface area to volume ratio and dense - most of the voxels inside the bounding box of the mask are foreground).
**Summary**
With that, I would appreciate the insight of the ITK community. I am leaning towards the modified zero flux Neumann boundary condition combined with the bounding box iteration. Any recommendations for this implementation?
---
### Author:  [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [March 21, 2019, 2:23pm UTC](https://discourse.itk.org/t/implementing-non-square-boundary-conditions/1690/2 "2019-03-21T14:23:40Z")
Your approach sounds fine to me. To speed up evaluation of boundary condition, you might want to calculate [DistanceMap](https://itk.org/Doxygen/html/classitk_1_1DanielssonDistanceMapImageFilter.html) inside the bounding box. You can then use the **vector map** to get to the closest boundary pixel in one step.
---
### Author:  [@Besler](https://discourse.itk.org/u/Besler)
#### Post date: [March 22, 2019, 8:17pm UTC](https://discourse.itk.org/t/implementing-non-square-boundary-conditions/1690/3 "2019-03-22T20:17:37Z")
Having a vector map of offsets is really brilliant! I could remove all branching by expanding the image by one layer and doing an offset lookup on every index. It’s more memory intensive, but it should be fast. It doesn’t extend well to Dirichlet or non-zero flux boundary conditions, but would be fast for this implementation.
Thanks @dzenanz!