Custom subarray of a multidimensional array

I have an array A which has a size [N1 N2 N3 3 3]. Array A can be considered as a rectangular box in three dimensions that for each of its element, an array of size [3 3] is defined, and here I call it B. In other words, B with size [3 3] is defined for each point in the domain.

I am trying to subdivide the rectangular box into small cubes of size d×d×d, where d is an even number. Also, there should be 50% overlap between adjacent cubes and . For example, let's say:

    A; % size(A) = [14 16 20 3 3];
    d = 4;

With the above example, the sampling will be as follows. The cubes along with with first dimension of A, i.e., 1:14 will include the following elements of A:

    cube1: 1, 2, 3, 4
    cube2: 3, 4, 5, 6
    cube3: 5, 6, 7, 8
    cube4: 7, 8, 9, 10
    cube5: 9, 10,11,12
    cube6: 11,12,13,14

along the second dimension of A, i.e., 1:16, we have:

    cube1: 1, 2, 3, 4
    cube2: 3, 4, 5, 6
    cube3: 5, 6, 7, 8
    cube4: 7, 8, 9, 10
    cube5: 9, 10,11,12
    cube6: 11,12,13,14
    cube7: 13,14,15,16

And finally along the third dimension of A, i.e., 1:20:

    cube1: 1, 2, 3, 4
    cube2: 3, 4, 5, 6
    cube3: 5, 6, 7, 8
    cube4: 7, 8, 9, 10
    cube5: 9, 10,11,12
    cube6: 11,12,13,14
    cube7: 13,14,15,16
    cube8: 15,16,17,18
    cube9: 17,18,19,20

Now Anew will have a size [6 7 9 3 3] and the value at each of its element will be the summation of each element of B over all the element of A that are in each cube. For example, Anew(1, 4, 9, 1, 1) is actually:

    A(elements of cube1, elements of cube4, elements of cube9, 1, 1)

or in other words, A(1:4, 7:10, 17:20). Therefore:

    Anew(1, 4, 9, 1, 1) = 0;
    for i = 1:4 % cube1 along the 1st dimension uses 1, 2, 3, 4
        for j = 7:10 % cube4 along the 2nd dimension uses 7, 8, 9, 10
            for k = 17:20 % cube9 along the 3rd dimension uses 17, 18, 19, 20
                Anew(1, 4, 9, 1, 1) = Anew(1, 4, 9, 1, 1) + A(i, j, k, 1, 1);
            end
        end
    end

In the above for loop, A(i, j, k, 1, 1) is B(1, 1) at point (i, j, k) in the original domain. Could someone help me how I should approach in a vectorize fashion ?

Answer 1

[Matt J]

I recommend the following, which uses the two FEX submissions interpMatrix and KronProd

    F=full(interpMatrix([1 1 1 1],1,11,2)).';

    P=F(1:6,1:14);
    Q=F(1:7,1:16);
    R=F(1:9,1:20);

    K=KronProd({P,Q,R,1},[1 2 3 4 4],[nan,nan,nan,3,3]);

    Anew=K*A;