Hvordan vi ser det skjulte

Her finde du information fra øvelsesdagen.

Filer

SCRIPT: getImages.m - show
function im = getImages(list)
% im = getImages;           - get all
% im = getImages([10:100]); - get projection0010-0100

if nargin == 0    % If no input arguments are given:
  list = 1:200;   % Select all images
end               %

im = zeros(500,500,length(list)); % Preallocate matrix
for i = list
    % Read image files and store in matrix
    file = fullfile('img_org',sprintf('projection%04d.tif',i-1));
    im(:,:,i) = imread(file);
end

SCRIPT: sinogram.m - show
function sg = sinogram(im,line)
% sg = sinogram(im)     - get all (!!!very memory intensive!!!)
% sg = sinogram(im, 10) - get line 10

if nargin < 2             % If line argument is not set:
    line = 1:size(im,1);  % Get all sinograms
end                       %

sg = zeros(size(im,2),size(im,3),length(line)); % Preallocate matrix
for k = 1:length(line)
    for i = 1:size(im,3)              % For each images:
        sg(:,i,k) = im(line(k),:,i);  % Put the selected line in the matrix
    end
end

SCRIPT: reconstruction_roll.m - show
function reconstruction_roll
% Animated roll through the sinograms

pixelshift = -9;

figure                        % Now figure
im = getImages;               % Load in all imagesc

subplot(1,2,1)                % Split figure into 1x2 grid, select element 1
imagesc(im(:,:,1));           % Show image 
colormap('gray');             % Set colorscale to gray
axis image                    % Set axis
hold on                       % Make sure image is not replaced by plot
p = plot([1,size(im,2)],[1,1],'r-'); % Plot red line

theta = linspace(0,360,size(im,3));

subplot(1,2,2)                % 1x2 grid, select element 2
sg = sinogram(im,1);          % Compute new sinogram
if pixelshift <= 0
  sg = [zeros(-2*pixelshift,size(sg,2)); sg]; % Compensate for shift in the physical setup
else
  sg = [sg; zeros(2*pixelshift,size(sg,2))]; % Compensate for shift in the physical setup
end
I = iradon(sg, theta);        % Compute reconstruction
h = imagesc(I);               % Show reconstruction image
colormap('gray');             % Set colorscaling to gray
axis image                    % Set axis
ht = title('line #1');        % Set title

for i = 1:size(im,1)
    set(p,'YData',i*[1,1]);   % Update red line Y-position
    sg = sinogram(im,i);      % Compute new sinogram
    if pixelshift <= 0
      sg = [zeros(-2*pixelshift,size(sg,2)); sg]; % Compensate for shift in the physical setup
    else
      sg = [sg; zeros(2*pixelshift,size(sg,2))]; % Compensate for shift in the physical setup
    end
    I = iradon(sg, theta);    % Compute reconstruction
    set(h,'Cdata',I);         % Update image
    set(ht,'String',sprintf('line #%d',i)) % Update title
    drawnow;                  % Make sure image is displayed before proceeding
end

SCRIPT: reconstruction_roll_interactive.m - show
function reconstruction_roll_interactive
% interactive roll through the reconstructions

global ax1 p im theta h ht line_selected pixelshift
pixelshift = -9;

fig = figure;                 % Now figure
im = getImages;               % Load in all imagesc

ax1 = subplot(1,2,1);         % Split figure into 1x2 grid, select element 1
imagesc(im(:,:,1));           % Show image 
colormap(ax1,'gray');         % Set colorscale to gray
axis(ax1,'image');            % Set axis
hold(ax1,'on');               % Make sure image is not replaced by plot
p = plot(ax1,[1,size(im,2)],[1,1],'r-'); % Plot red line
line_selected = 1;

ax2 = subplot(1,2,2);         % 1x2 grid, select element 2
sg = sinogram(im,1);          % Compute new sinogram
if pixelshift <= 0
  sg = [zeros(-2*pixelshift,size(sg,2)); sg]; % Compensate for shift in the physical setup
else
  sg = [sg; zeros(2*pixelshift,size(sg,2))]; % Compensate for shift in the physical setup
end
theta = linspace(0,360,size(im,3));

I = iradon(sg, theta);
h = imagesc(I);               % Show sinogram image
colormap(ax2,'gray');         % Set colorscaling to gray
axis(ax2,'image');            % Set axis
ht = title(ax2,'1');          % Set title

% Dynamic callbacks (Advanced)
set(fig,'windowbuttondownfcn',@(s,e) mouse_track(s,e,true));
set(fig,'windowbuttonmotionfcn',@(s,e) mouse_track(s,e,false));

end

function plot_sinogram(i)
  global p im theta h ht pixelshift
  set(p,'YData',i*[1,1]);  
  sg = sinogram(im,i); 
  if pixelshift <= 0
    sg = [zeros(-2*pixelshift,size(sg,2)); sg]; % Compensate for shift in the physical setup
  else
    sg = [sg; zeros(2*pixelshift,size(sg,2))]; % Compensate for shift in the physical setup
  end
  I = iradon(sg, theta);
  set(h,'CData',I);
  set(ht,'String',sprintf('line #%d',i))
  drawnow;
end

function mouse_track(s,e,update)
  global ax1 p im theta h ht line_selected
  pos = get(ax1,'currentpoint');
  pos = pos(1,1:2);
  xlims = xlim(ax1);
  ylims = ylim(ax1);
  if pos(1) < xlims(1) || pos(1) > xlims(2) || pos(2) < ylims(1) || pos(2) > ylims(2)
    set(p,'YData',line_selected*[1,1]);
    return
  end
  line = floor(pos(2));
  set(p,'YData',line*[1,1]);
  if update
    line_selected = line;
    plot_sinogram(line);
  end
end

SCRIPT: sinogram_roll.m - show
function sinogram_roll
% Animated roll through the sinograms

figure                        % Now figure
im = getImages;               % Load in all imagesc

subplot(1,2,1)                % Split figure into 1x2 grid, select element 1
imagesc(im(:,:,1));           % Show image 
colormap('gray');             % Set colorscale to gray
axis image                    % Set axis
hold on                       % Make sure image is not replaced by plot
p = plot([1,size(im,2)],[1,1],'r-'); % Plot red line

subplot(1,2,2)                % 1x2 grid, select element 2
for i = 1:size(im,1)
    set(p,'YData',i*[1,1]);   % Update red line Y-position
    sg = sinogram(im,i);      % Compute new sinogram
    imagesc(sg);              % Show sinogram image
    colormap('gray');         % Set colorscaling to gray
    axis image                % Set axis
    title(sprintf('line #%d',i)) % Set title
    drawnow;                  % Make sure image is displayed before proceeding
end

SCRIPT: sinogram_roll_interactive.m - show
function sinogram_roll_interactive
% interactive roll through the sinograms

global im ax1 p h ht line_selected

fig = figure;                 % Now figure
im = getImages;               % Load in all imagesc

ax1 = subplot(1,2,1);         % Split figure into 1x2 grid, select element 1
imagesc(im(:,:,1));           % Show image 
colormap(ax1,'gray');         % Set colorscale to gray
axis(ax1,'image');            % Set axis
hold(ax1,'on');               % Make sure image is not replaced by plot
p = plot(ax1,[1,size(im,2)],[1,1],'r-'); % Plot red line
line_selected = 1;

ax2 = subplot(1,2,2);         % 1x2 grid, select element 2
sg = sinogram(im,1);          % Compute new sinogram
h = imagesc(sg);              % Show sinogram image
colormap(ax2,'gray');         % Set colorscaling to gray
axis(ax2,'image');            % Set axis
ht = title(ax2,'1');          % Set title
   
% Dynamic callbacks (Advanced)
set(fig,'windowbuttondownfcn',@(s,e) mouse_track(s,e,true));
set(fig,'windowbuttonmotionfcn',@(s,e) mouse_track(s,e,false));

end

function plot_sinogram(i)
  global im ax1 p h ht
  set(p,'YData',i*[1,1]);  
  sg = sinogram(im,i);     
  set(h,'CData',sg);
  set(ht,'String',sprintf('line #%d',i))
  drawnow;
end

function mouse_track(s,e, update)
  global im ax1 p h ht line_selected
  pos = get(ax1,'currentpoint');
  pos = pos(1,1:2);
  xlims = xlim;
  ylims = ylim;
  if pos(1) < xlims(1) || pos(1) > xlims(2) || pos(2) < ylims(1) || pos(2) > ylims(2)
    set(p,'YData',line_selected*[1,1]);
    return
  end
  line = floor(pos(2));
  set(p,'YData',line*[1,1]);
  if update
    line_selected = line;
    plot_sinogram(line);
  end
end

Video

Her vil ligge skanningsdata fra Nanoteket. Det ligger i fem formater: (1) rådata fra programmet, (2) png-version af projektionerner, (3) sinogrammer/radontransformationen for alle vandrette snit i billedet, (4) rekonstruktion med back projection for hvert sinogram, (5) rekonstruktion med filtered back projection for hvert sinogram.

Dertil er tre informationspunkter om dataene: (1) Akseforskydning i pixels, (2) Udfyldningen i pixels og X-ray kilden i maskinen.

Akseforskydning: Når man bruger iradon funktionen i MATLAB (eller Octave) til at rekonstruere et billede fra et sinogram er det vigtigt at rotationsaksen er i midten af sinogrammet. Hvis det ikke er tilfældet må vi forskyde det ved at indsætte nogle rækker af pixels enten i toppen eller i bunden af sinogrammet. Negativ pixelværdi betyder de er indsat i toppen, positiv i bunden. Forskydningen er her bestemt ved visuel inspektion af rekonstruktionen, men vi kan også bestemme den ud fra 360-graders sinogrammet ved at finde ud af hvilken linie der ligger mellem top og bund af bølgen. Eksempel: Billedet er 500px højt, men rotationsaksen ligger omkring række 241, dvs. forskydning på -9px fra midten af billedet, række 250. For at korrigere dette tilføjer vi 2×9px = 18px i toppen af billedet.

Udfyldning: På grund af den måde iradon funktionen fungerer på er det rekonstruerede område bestemt ud af højden af sinogrammet. Så for at få mest muligt billede med blev det tilføjet nogle rækker af fyldpixels i toppen og bunden af sinogrammet.

Fjerbold data

Information

Akseforskydning: -9px
Udfyldning: 40px
X-ray kilde: ???

Øvelsesdag

Filer

Octave

Skanningsdata

Fjerbold data