function d2 = dist2(data, centres)

[ndata, ddim] = size(data);
[ncentre, cdim] = size(centres);


d2 = (ones(ncentre,1) * sum((data.^2)',1))' + ...
ones(ndata,1) * sum((centres.^2)',1) - ...
2.*(data*(centres'));