Move all counts of the image in random direction

I have the following problem. I have 8-bit image. It represents photon counts, so each pixel is capable to register up to 255 photons. I need to introduce some random error - each photons trajectory is disturbed, so that finaly instead of e.g. (150, 100) it is regiesterd at (150+random1, 100+random2). Random1 and random2 both have Poisson distirbution and are small, ~2 pixels.

Notice, that is is not enough to displace each pixel in random direction, since each count from each pixel has to be displaced separately.

What is the easiest and fastest way to code it?

Best regards, Alex

Answer 1

[Image Analyst]

You need to scan your original image, and if the value is not zero (meaning that there is a photon in the original CCD well or whatever the image represents) then you need to calculate the new location (row,col) by adding a deltax and deltay taken from a Poisson distribution. Then add the photon to the output location and (optionally) remove it from the input location.

outputImage = zeros(rows, columns);
for col = 1 : columns
    for row = 1 : rows
        if originalImage(row, col) > 1
            newRow = round(row + random('Poisson', pd));
            newCol = round(col + random('Poisson', pd));
            % Increment the number of pixels at this output image location.
            outputImage(newRow, newCol) = outputImage(newRow, newCol) + 1;
            % Decrement the number of photons at this output image location (only if you want to)
            originalImage(row, col) = originalImage(row, col) - 1;
        end
    end
end

See documentation for random() in the Statistics and Machine Learning Toolbox for info on how to get pd, the probability distribution object.

[Alex Kurek]

There was one small issue (should be while instead of if), but it is working:

tic
for col = 1 : columns
    for row = 1 : rows
        while imageTestTemp(row, col) > 0
            newRow = round(row + randRow = random('Poisson', pd));
            newCol = round(col + randCol = random('Poisson', pd));
            if newRow>0 && newRow<columns+1    &&    newCol>0 && newCol<rows+1
                outputImage(newRow, newCol) = outputImage(newRow, newCol) + 1;
            end
            imageTestTemp(row, col) = imageTestTemp(row, col) - 1;
        end
    end
    disp([num2str(col) ' / ' num2str(columns)]);
end
toc

The problem is, that it is working VERY slow. Is there some way to speed it up, e.g. vectorize it?

[Image Analyst]

You could perhaps get all the Poisson distances in 256 arrays in advance instead of getting them one number at a time.

[Alex Kurek]

Yes, it worked. Like this it is >300x faster:

countNo = sum(sum(imageTest));
randContainer = round(normrnd(0, pd, [1 2*countNo]));
i = 1;

tic
for col = 1 : columns
    for row = 1 : rows
        while imageTestTemp(row, col) > 0
            newRow = row + randContainer(i);
            newCol = col + randContainer(i+1);
            i = i+2;
            if newRow>0 && newRow<columns+1    &&    newCol>0 && newCol<rows+1
                outputImage(newRow, newCol) = outputImage(newRow, newCol) + 1;
            end
            imageTestTemp(row, col) = imageTestTemp(row, col) - 1;
        end
    end
end
toc