I’m wanting to implement a filter using itkDenseFiniteDifferenceImageFilter 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 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 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 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 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?