For starters, put semi-colons at the end of the delta(etc) = etc lines so that intermediate stuff doesn't print to the screen. Doing just that cuts the execution time to under 1 sec on my computer. E.g.,
delta(k,i,j)=1; % <-- appended semi-colon
else delta(k,i,j)=0; <-- appended semi-colon
For more speed improvements, could work on vectorizing the loops. This should get you all the speed improvements needed without resorting to some type of GPU conversion. E.g., here is one way to vectorize the code with simple brute force application of the bsxfun function to each of your operations:
pc1_1 = pc1(:); % Convert to column vector in 1st dimension
pc2_1 = pc2(:); % Convert to column vector in 1st dimension
x1_2 = reshape(x1,1,40,1); % Convert to "vector" in 2nd dimension
x2_3 = reshape(x2,1,1,40); % Convert to "vector" in 3rd dimension
arg1 = bsxfun(@le,x1_2-c1,pc1_1);
arg2 = bsxfun(@lt,pc1_1,x1_2+c1);
arg3 = bsxfun(@le,x2_3-c2,pc2_1);
arg4 = bsxfun(@lt,pc2_1,x2_3+c2);
arg12 = bsxfun(@and,arg1,arg2);
arg34 = bsxfun(@and,arg3,arg4);
delta = bsxfun(@and,arg12,arg34);
delta(:,40,:) = bsxfun(@eq,pc1_1,x1_2+c1);
delta(:,:,40) = bsxfun(@eq,pc2_1,x2_3+c2);
There may be a way to simplify this even further, but the above code runs in about 0.02 sec on my machine so I stopped working on it.
