Take a look at my peak_nii toolbox in NITRC. It will generate a map of clusters where each cluster has a unique value. Then, you can do something similar to the following:
hdr=spm_vol(clusterimage);
img=spm_read_vols(hdr);
mimg=spm_read_vols(spm_vol(AALimage));
if ~all(size(img)==size(mimg))
disp('Images are not the same size')
error()
end
for ii=1:max(img(:))
for jj=1:max(mimg(:))
oimg=zeros(size(img));
oimg(img==ii && mimg==jj)=1;
hdr.fname=['Cluster_' num2str(ii) '_AALRegion_' num2str(jj) '.nii'];
spm_write_vol(hdr,img)
end
end
Now you have a file for each cluster split by anatomical region.