the following code is easier and faster:
X_hist=sum(image,1);
Y_hist=sum(image,2);
X=1:N; Y=1:M;
centX=sum(X.*X_hist)/sum(X_hist);
centY=sum(Y'.*Y_hist)/sum(Y_hist);
Besides, you may use "meshgrid" command instead of defining x variable in your code (i.e. x=meshgrid(1:columns,1:rows);)
Comment only