Datacube indexing speed up

I use some weighted probability routines and the following code is to slow. Is there a way to optimize it? I already preallocated everything.

function normalizedFrame = normalizer (currentFrame, normalizedFrame, oneLevelFrame, dataCube, maxADU)   %#codegen
%
for i = maxADU:-1:1
    oneLevelFrame(currentFrame==i) = 1;
    dataCube(:,:,i) = oneLevelFrame;
    oneLevelFrame   = zeros(31, 31);
    currentFrame(currentFrame==i) = currentFrame(currentFrame==i)-1;
end
%
possitionsIn3dCube = find(dataCube);
[row, col, ~] = ind2sub(size(dataCube), possitionsIn3dCube);
index = randi(length(row));
photonPossition = [row(index), col(index)];
normalizedFrame(photonPossition(1), photonPossition(2)) = 1;

Best regards, Alex

Answer 1

[Walter Roberson]

A little more efficient would be

for i = maxADU:-1:1
  idx == currentFrame == i;
  oneLevelFrame(idx) = 1;
  dataCube(:,:,i) = oneLevelFrame;
  oneLevelFrame   = zeros(31, 31);
  currentFrame(idx) = currentFrame(idx)-1;
end

Could you confirm that there is no further code in the routine than what you show here? In particular, that there is no further code that uses currentFrame? Because if that is true, there is a more efficient version that does not use loops.

[Alex Kurek]

I confirm. There is no more code using currentFrame, it's just a temporary variable and used only in the function here

[Walter Roberson]

The logic of the routine is impeded by the need to match the size and datatypes and initial values of the input arrays. If the input conditions had been documented some of this might not be needed.

function normalizedFrame = normalizer (currentFrame, normalizedFrame, oneLevelFrame, dataCube, maxADU)   %#codegen

nottoobig = repmat(currentFrame <= maxADU, 1, 1, maxADU);
bigenough = bsxfun(@ge, currentFrame, reshape(1:maxADU, 1, 1, []));
local_dataCube = bigenough & nottoobig;   %logical

%the input dataCube might not be logical. And that matters
%because the initial value of oneLevelFrame is copied in
%in the positions where there is no match, and the datatype of
%oneLevelFrame is unknown and might not match dataCube. For
%example, dataCube might be uint8 but oneLevelFrame might be
%double and have fractional values. We have to ensure that the
%datatype conversions are done properly to match the original code

local_dataCube = cast(local_dataCube, class(dataCube));
idx = oneLevelFrame ~= maxADU;
dataCube_at_max = local_dataCube(:,:,maxADU);
dataCube_at_max(idx) = oneLevelFrame(idx);
local_dataCube(:,:,maxADU) = dataCube_at_max;

if size(dataCube,3) > maxADU
  dataCube(:,:,1:maxADU) = local_dataCube;
else
  dataCube = local_dataCube;
end

positionsIn3dCube = find(dataCube);
[row, col, ~] = ind2sub(size(dataCube), positionsIn3dCube);
index = randi(length(row));
photonPosition = [row(index), col(index)];
normalizedFrame(photonPosition(1), photonPosition(2)) = 1;

I suspect there might be a more efficient way of doing the calculation, but it is going to depend upon the initial conditions:

• is the localFrame that is passed in always type logical() ?
• if it is not, then is localFrame always a double precision array of 0?
• if it is not, then is localFrame double precision and contains only the values 0 and 1?
• is the dataCube that is passed in always type logical() ?
• if it is not, then is dataCube always a double precision array of 0?
• if it is not, then is dataCube double precision and contains only the values 0 and 1?
• is size(dataCube,3) ever greater than maxADU, so that the calculation is over a subset of dataCube but the entire original dataCube needs to be considered for the normalization phase?
• can currentFrame contain any values that are greater than maxADU ?
• can currentFrame contain any values that are not integers but might be between 0 and maxADU?
• is maxADU guaranteed to be a non-negative integer?
• is it correct that you want the photon position probability to increase according to the value in currentFrame?

[Alex Kurek]

Thank you for all this. I just finished testing everything. I also learned a lot!

1. the "A little more efficient" version is little less than 30% faster
2. The long version is less than 50% faster from "A little more efficient" beginning from maxADU = 50. Otherwise it is slower, so I implemented it using if / else.

• localFrame and dataCube are double.
• is size(dataCube,3) ever greater than maxADU - NO
• can currentFrame contain any values that are greater than maxADU - NO
• can currentFrame contain any values that are not integers but might be between 0 and maxADU? - NO, always integers
• maxADU is between 1 and 255, but usually <5
• you want the photon position probability to increase according to the value in currentFrame? YES, its true

I changed two things only: This:

idx == currentFrame == i;

to this:

idx = currentFrame == i;

and this:

local_dataCube = cast(local_dataCube, class(dataCube));

to this:

local_dataCube = cast(local_dataCube, 'like', dataCube);

in the latter case Matlab told mi to do it