Optimizing voxel merge algorithm
Try to determine the complexity of the algorithms you are implementing. The code you wrote above has a complexity of O(N⁹), where N is equal to grid.dimension. This is because you have 6 nested for loops, and then voxels.all{} which itself also has to loop over all three dimensions of the box voxels.

Grouping voxels by their int value is a good idea. You correctly worry about the case of voxels with the same value being spread out. The solution to that is to divide the box into connected voxels with the same int value. For example, if you had only one dimension, then if you have the values [0, 1, 2, 2, 1, 2], then you would divide this into 5 groups: [0], [1], [2, 2], [1] and [2]. It should be possible to do this in O(N³) time.

Now that you have a set of smaller shapes, the function that decomposes those shapes into boxes should run a lot faster. The question is now whether you want to have that function as fast as possible, or as accurate as possible. The function as you wrote does not necessarily result in the least amount of boxes. Apart from the fact that it doesn't actually try out boxes in a strict order from largest to smallest, even trying the largest boxes first won't necessarily result in the least amount of boxes being produced. Take for example this 2D array:

1 0 0 1
1 1 1 1
1 0 0 1


The voxels with value 1 form a H-shape. If you remove the largest possible box, you end up with:

1 0 0 1
x x x x
1 0 0 1


Now you are left with 4 1x1 boxes, so in total this will result in 5 boxes for the H-shape. If instead you would start by taking out the left and right side, you are left with the middle bar, resulting in just 3 boxes. Finding the solution with the least amount of boxes is actually an NP-complete problem, and an algorithm that would produce that solution would take even longer to run that yours. Given that your algorithm already does not give acceptable performance, you should instead focus on implementing the fastest algorithm that produces acceptable results.

The common approach is to start by finding strips of adjacent voxels in the x-direction. Then try to find adjacent strips with the same size and x-position in the y-direction, and merge them into squares. Finally do the same in the z-direction and merge adjacent squares into boxes. This should be O(N³) as well. If you know that some shapes are more common than others, you could tweak this algorithm to reduce the number of boxes.