Image manipulation with Matlab - intensity and tif image

Question

I have to analyse a set of images, and these are the operations I need to perform:

1. sum another set of images (called open beam in the code), calculate the median and rotate it by 90 degrees;
2. load a set of images, listed in the file "list.txt";
3. the images have been collected in groups of 3. For each group, I want to produce an image whose intensity values are 3 times the image median above a certain threshold and otherwise equal to the sum of the intensity values;
4. for each group of three images, subtract the open beam median (calculated in 1.) from the combined image (calculated in 3.)

Considering one of the tifs produced using the process above, I have that the maximum value is 65211, which is not 3* the median for the three images of the corresponding group (I checked considering the pixel position). Do you have any suggestion on why this happens, and how I could fix it?

The code is reported below. Thanks!

%Here we calculate the average for the open beam
clear;
j = 0;
for i=1:5
   s = sprintf('/Users/Alberto/Desktop/Midi/17_OB_2.75/midi_%04i.fits',i);
   j = j+1;
   A(j,:,:) = uint16(fitsread(s));
end

OB_median = median(A,1);
OB_median = squeeze(OB_median);
OB_median_rot=rot90(OB_median);

%Here we calculate, for each projection, the average value from the three datasets

%Read list of images from text file
fid = fopen('/Users/Alberto/Desktop/Midi/list.txt', 'r');
a = textscan(fid, '%s');
fclose(fid);

%load images
j = 0;
for i = 1:1:42 %556 entries; 543 valid values
      s = sprintf('/Users/Alberto/Desktop/Midi/%s',a{1,1}{i,1});
      j = j+1;
      A(j,:,:) = uint16(fitsread(s));
end

threshold = 80 %This is a discretional number. I put it after noticing 
%that we get the same number of pixels with a value >100 if we use 80 or 50.

k = 0;
for ii = 1:3:42
    N(1,:,:) = A(ii,:,:);
    N(2,:,:) = A(ii+1,:,:);
    N(3,:,:) = A(ii+2,:,:);
    median_N = median(N,1);
    median_N = squeeze(median_N); 
    B(:,:) = zeros(2160,2592);
    for i = 1:1:2160
        for j = 1:1:2592
            RMS(i,j) = sqrt((double(N(1,i,j).^2) + double(N(2,i,j).^2) + double(N(3,i,j).^2))/3);
            if RMS(i,j) > threshold
               %B(i,j) = 30;
               B(i,j) = 3*median_N(i,j);
            else
               B(i,j) = A(ii,i,j) + A(ii+1,i,j) + A(ii+2,i,j);
               %B(i,j) = A(ii,i,j);
            end
        end
    end
    k = k+1;
    filename = sprintf('/Users/Alberto/Desktop/Midi/Edited_images/Despeckled_images/despeckled_image_%03i.tif',k);
    %Now we rotate the matrix
    B_rot=rot90(B);
    imwrite(B_rot, filename);
    %imwrite(uint16(B_rot), filename);
    %Now we subtract the OB median
    B_final_rot = double(B_rot) - 3*double(OB_median_rot);
    filename = sprintf('/Users/Alberto/Desktop/Midi/Edited_images/Final_image/final_image_%03i.tif',k);
    imwrite(uint16(B_final_rot), filename);
end

Answer 1

The maximum integer that can be represented by the uint16 data type is

>> a=100000; uint16(a)

ans =

  65535

To circumvent this limitation you need to rescale your data as type double and adjust the range (the image contrast) to agree with the limits imposed by the uint16 data type, before saving as uint16.