code: def _get_blocks(self, mask):
          
          blocks = []
          outerblk = self.max_blk_edge + 2*self.sl_rad
          for i in range(0, mask.shape[0], self.max_blk_edge):
              for j in range(0, mask.shape[1], self.max_blk_edge):
                  for k in range(0, mask.shape[2], self.max_blk_edge):
                      block_shape = mask[i:i+outerblk,
                                         j:j+outerblk,
                                         k:k+outerblk
                                         ].shape
                      if np.any(
                          mask[i+self.sl_rad:i+block_shape[0]-self.sl_rad,
                               j+self.sl_rad:j+block_shape[1]-self.sl_rad,
                               k+self.sl_rad:k+block_shape[2]-self.sl_rad]):
                          blocks.append(((i, j, k), block_shape))
          return blocks
docs: Divide the volume into a set of blocks
  
          Ignore blocks that have no active voxels in the mask
  
          Parameters
          ----------
  
          mask: a boolean 3D array which is true at every active voxel
  
          Returns
          -------
  
          list of tuples containing block information:
             - a triple containing top left point of the block and
             - a triple containing the size in voxels of the block