Simple Image Processing Optimisation

I have a simple image processing application that is working fine, but I think it could be vectorised and sped up substantially. I've had a shot but am having trouble with it.

The idea is to take a series of images and get

- timex: time exposure (the average of each pixel over time)

- immax: each pixel with the maximum intensity over time

- immix: each pixel with the minimum intensity over time

Working code:

im_stack=[timestep,x_pixel,y_pixel,[red,green,blue,intensity]]

timex = nan(x,y,3); immax = timex; immin = timex; %preallocate
for i = 1:x
for j = 1:y
  for k = 1:z
    timex(i,j,k)= mean(im_stack(:,i,j,k));
  end
  [~,a] = max(im_stack(:,i,j,4));%from intensity values
  immax(i,j,1:3)= im_stack(a,i,j,1:3);
  [~,a] = min(im_stack(:,i,j,4));
  immin(i,j,1:3)= im_stack(a,i,j,1:3);
end
end
timex=uint8(timex);
immax=uint8(immax);
immin=uint8(immin);

[Walter Roberson]

Could you clarify what is special about im_stack(:,:,:,4) compared to im_stack(:,:,:,1:3) ? And what values is z ?

[Ben Modra]

Hi Walter,

x= number of x_pixels

y = number of y_pixels

z = number of colour bands (3)

im_stack(:,:,:,4) is the intensity values of each pixel generated using rgb2grey. im_stack(:,:,:,1:3) is the three colour bands. The selection of the brightest pixel is based on the intensity, but I need to return the pixel colour not the greyscale.

[Ben Modra]

So I've progressed somewhat but am stuck on one line. This is a simpler example but if I can get the indexing right then I can roll out to higher dimensions.

The same idea in this example, just with two images. I need to use the index rather than just max, because the index will used to recall the equivalent colour values.

   tmpint = cat(3, rgb2grey(im1), rgb2grey(im2)); %stack the intensity of images together
[~,I] = max(tmpint,[],3); %get indices of highest intensity
% I has indicies of the image stack, but how to be used in the line below??
max_intensity = reshape(tmpint(**I**),[x_pix,y_pix]); % new BW image of highest intensity pixels
%
tmpcol = cat(4,im1,im2); %stack colour images
J=repmat(I,[1,1,3,1]);
% J needs to now have the indicies of the pixel colours
max_col = reshape(tmpcol(J),[x_pix,y_pix,3]);
%
image(max_intensity)
image(max_col)

Answer 1

[Walter Roberson]

As a start:

With no loops needed,

timex = squeeze( mean(im_stack(:,:,:,1:z)) );

And if z is the same as size(im_stack,4) then it would simplify to

timex = squeeze( mean(im_stack) );

[Ben Modra]

Thanks! I knew there was an elegant solution, I was just not seeing how it would average along one dimension.

The max& min will be different though. Do I need to reshape to get the indexes somehow?