You start with what appears to be a RGB image and then you split it into 4 RGB images. If the input is color, the output should still be color. Is the question simply how to split it into 16 instead of 4? I'm going to assume that.
You can split into 16 images just the same as you've split into 4, though you're counting on the image height and width being both integer-divisible by 4; otherwise, you'll tend to wind up with errors as you'll have subscripts that are either fractional or out of bounds. You will have to take that into consideration if you want your subscript range calulations to be robust.
If you want a simpler method for splitting, the MIMT (look on the file exchange) has a function imdetile() that can split flat images into 4D arrays (which can then be separated into individual images if desired). If the image needs to be reassembled after some processing, that can be done with imtile() from the same toolbox. For example:
inpict=imread('sources/RGB16Million.png');
wstack=imdetile(inpict,[4 4]);
[imsize(inpict); imsize(wstack)]
This reveals that we've divided a flat image into multiframe image.
ans =
4096 4096 3 1
1024 1024 3 16
We can work on this multiframe image however we choose -- as a single array, or by splitting it up into 16 individual flat images. When we're done, we can rearrange it back into a 4D array and run imtile() to reassemble it into a single flat image.
outpict=imtile(wstack,[4 4]);
imsize(outpict)
And our image is the same size it used to be.
Of course, our image was integer-divisible by 4 to begin with. What if it wasn't? Let's divide it into 25 images.
ans =
4096 4096 3 1
819 819 3 25
ans =
4095 4095 3 1
By default, imdetile() either pads or crops edges (whichever is least) to make the image integer-divisible. As a consequence of using defaults, the output image is not the same size. This can be remedied by explicitly choosing an appropriate option and then adjusting the output accordingly.
inpict=imread('sources/RGB16Million.png');
wstack=imdetile(inpict,[5 5],'fittype','scale');
outpict=imresize(imtile(wstack,[5 5]),imsize(inpict,2));
In this case, we used a scaling method to make it integer-divisible and then we rescaled the output to match the input geometry. If you don't want to risk the interpolation error, you can use the 'grow' fit type and then crop off the excess edge vectors using cropborder() instead.
Of course, that's all using MIMT tools. Is there a purely Matlab/IPT way to do it? Probably, but you'll have to work out all the tilings manually.