GPU is the slowest on my PC. I knew for-loop is now away fast, but I'm still amazed by how fast it is. There is really a remarkable progress made by TMW on Excution Engine last few releases.
- For-loop time = 0.010799
- Vectorized CPU time = 0.012376
- Vectorized GPU time = 0.069613
Code:
N = 1000;
Image = peaks(N);
Xmin = 1; Xmax = N;
Ymin = 1; Ymax = N;
CentX = 0;
CentY = 0;
Resolution = 1/N;
BinDensity = 100;
MaxTargetBin = 100;
tic
Profile = zeros(MaxTargetBin,1);
SumBinSignal = zeros(MaxTargetBin,1);
for X=Xmin:1:Xmax
for Y=Ymin:1:Ymax
if ~isnan(Image(Y,X)) %NaN pixels are just ignored
Dist = Resolution * ( (X-CentX)^2 + (Y-CentY)^2 )^0.5; %distance of point to centre measured in mm
TargetBin= ceil(Dist * BinDensity); %convert distance from mm to "bins/rings"
if (TargetBin>=1) && (TargetBin<=MaxTargetBin)
Profile(TargetBin)= Profile(TargetBin)+1; %record number of pixels in bin/ring
SumBinSignal(TargetBin)= SumBinSignal(TargetBin) + Image(Y, X); % record total signal in ring/bin
end
end
end
end
fprintf('For-loop time = %1.6f\n', toc)
% Vectorized code, CPU
tic
X = Xmin:Xmax;
Y = (Ymin:Ymax)';
Dist = Resolution * sqrt( (X-CentX).^2 + (Y-CentY).^2 );
TargetBin= ceil(Dist * BinDensity);
Im = Image;
keep = (TargetBin>=1) & (TargetBin<=MaxTargetBin) & ~isnan(Im);
TargetBin = TargetBin(keep);
Profile2 = accumarray(TargetBin, 1, [MaxTargetBin,1]);
if ~isequal([length(Y) length(X)], size(Im))
Im = Image(Y,X); % indexing is costly
end
SumBinSignal2 = accumarray(TargetBin, Im(keep), [MaxTargetBin,1]);
fprintf('Vectorized CPU time = %1.6f\n', toc)
if gpuDeviceCount("available")
% Vectorized code, GPU
tic
X = gpuArray(Xmin:Xmax);
Y = gpuArray((Ymin:Ymax))';
Dist = Resolution * sqrt( (X-CentX).^2 + (Y-CentY).^2 );
TargetBin= ceil(Dist * BinDensity);
Im = gpuArray(Image);
keep = (TargetBin>=1) & (TargetBin<=MaxTargetBin) & ~isnan(Im);
TargetBin = TargetBin(keep);
Profile3 = gather(accumarray(TargetBin, 1, [MaxTargetBin,1]));
if ~isequal([length(Y) length(X)], size(Im))
Im = Im(Y,X); % indexing is costly
end
SumBinSignal3 = gather(accumarray(TargetBin, Im(keep), [MaxTargetBin,1]));
fprintf('Vectorized GPU time = %1.6f\n', toc)
end
% Compare results
%norm(Profile-Profile2)
norm(SumBinSignal-SumBinSignal2)
if gpuDeviceCount("available")
%norm(Profile-Profile3)
norm(SumBinSignal-SumBinSignal3)
end
