Mediasys Software
Relazione Script MatLab su Sottocampionamento e Interpolazione 3D

Contenuti

Parte 1: sottocampionamento 8-1. Carico la matrice originale.

clear all;
clc;

D = zeros (256,256,20);

for i = 1:20
filename = ['brainMR/brain_' num2str(i,'%03d') '.dcm'];
  D(:,:,i) = dicomread(filename);
end

Parte 1: Ricavo le dimensioni della matrice e calcolo la Matrice sottocampionata SENZA cicli FOR.

  x=numel(D(1,:,1));  % ricavo le dimensioni della matrice
  y=numel(D(:,1,1));
  z=numel(D(1,1,:));
  x_new=x/2;          % calcolo le dimensioni della matrice sottocampionata
  y_new=y/2;
  z_new=z/2;

  P = zeros(y_new,x_new,z_new);  %  alloco la matrice sottocampionata

  % Versione di calcolo senza cicli FOR.

P =  (D(:,:,1:2:end-1) + D(:,:,2:2:end)) ;
P =  (P(1:2:end-1,:,:) + P(2:2:end,:,:)) ;
P =  (P(:,1:2:end-1,:) + P(:,2:2:end,:));

% Versione OLD con 3 Cicli FOR.
 % for zi=1:z_new
   %   for yi=1:y_new
       %   for xi=1:x_new
         %        P(yi,xi,zi)= (D(2*yi-1,2*xi-1,2*zi-1)+ D(2*yi,2*xi-1,2*zi-1)+ D(2*yi-1,2*xi,2*zi-1)+ D(2*yi,2*xi,2*zi-1)+  D(2*yi-1,2*xi-1,2*zi)+  D(2*yi,2*xi-1,2*zi)+  D(2*yi-1,2*xi,2*zi)+   D(2*yi,2*xi,2*zi));
       %   end
   %   end
 % end

   P=P./8;

Parte 1: Calcolo della griglia su cui interpolare

   [x_int,y_int,z_int] = meshgrid(0.75:0.5:y_new+0.25, 0.75:0.5:x_new+0.25, 0.75:0.5:z_new+0.25);

% Calcolo delle interpolazioni con i vari metodi

   D_int_nearest = uint16(interp3(P,x_int,y_int,z_int,'nearest'));
   D_int_linear = uint16(interp3(P,x_int,y_int,z_int,'linear'));
   D_int_cubic = uint16(interp3(P,x_int,y_int,z_int,'cubic'));

% Spiegare bene il motivo del passaggio a double, in pratica prima si e'
% portatto tutto a uint16 per liberarsi dei NaN , poi si e' tornati a
% double per fare i conti meglio.

   D_int_nearest = double(D_int_nearest);
   D_int_linear = double(D_int_linear);
   D_int_cubic = double(D_int_cubic);

Parte 1: Calcolo del PSNR nel caso di interpolazione con metodo nearest

   MSE_nearest_temp = (D(:) - D_int_nearest(:));
   MSE_nearest_temp = (MSE_nearest_temp).^2;
   MSE_nearest_temp = (MSE_nearest_temp)./(x*y*z);
   MSE_nearest = sum(MSE_nearest_temp);

   PSNR_nearest = 10 * log10((max(D(:))^2)/MSE_nearest);

   clear MSE_nearest_temp;

Parte 1: Calcolo del PSNR nel caso di interpolazione con metodo linear

  D_temp = D(:,:,2:1:end-1);
  D_int_linear = D_int_linear (:,:,2:1:end-1);
  MSE_linear_temp = (D_temp(:) - D_int_linear(:));
  MSE_linear_temp = (MSE_linear_temp).^2;
  MSE_linear_temp = (MSE_linear_temp)./(x*y*(z-2));
  MSE_linear = sum(MSE_linear_temp);

   PSNR_linear = 10 * log10((max(D(:))^2)/MSE_linear);

   clear MSE_linear_temp;

Parte 1: Calcolo del PSNR nel caso di interpolazione con metodo cubic

  D_int_cubic = D_int_cubic(:,:,2:1:end-1);
  MSE_cubic_temp = (D_temp(:) - D_int_cubic(:));
  MSE_cubic_temp = (MSE_cubic_temp).^2;
  MSE_cubic_temp = (MSE_cubic_temp)./(x*y*(z-2));
  MSE_cubic = sum(MSE_cubic_temp);

   PSNR_cubic = 10 * log10((max(D(:))^2)/MSE_cubic);

   clear MSE_cubis_temp;

Parte 1: Un pò di immagini

% montage immagine  originale
% Mi occorre una immagine copia in uint16 per plottare D, perchè D mi
% occorre in originale in Double per i calcoli della seconda parte

D_montage = uint16(D);

for i = 1:size(D_montage,3)
    test(:,:,1,i) = gray2ind(imadjust(D_montage(:,:,i)),256);
end
figure(1);
montage(test(:,:,1,1:2:20),gray(256));
set(gcf,'position',[10 100 500 500]);


% montage immagine  interpolata con metodo nearest
D_int_nearest = uint16(D_int_nearest);

for i = 1:size(D_int_nearest,3)
    test(:,:,1,i) = gray2ind(imadjust(D_int_nearest(:,:,i)),256);
end
figure(2);
montage(test(:,:,1,1:2:20),gray(256));
set(gcf,'position',[210 100 500 500]);

% montage immagine  interpolata con metodo linear
D_int_linear = uint16(D_int_linear);

for i = 1:size(D_int_linear,3)
   test(:,:,1,i) = gray2ind(imadjust(D_int_linear(:,:,i)),256);
 end
 figure(3);
 montage(test(:,:,1,1:2:18),gray(256));
 set(gcf,'position',[410 100 500 500]);

%  montage immagine  interpolata con metodo cubic
D_int_cubic = uint16(D_int_cubic);

 for i = 1:size(D_int_cubic,3)
     test(:,:,1,i) = gray2ind(imadjust(D_int_cubic(:,:,i)),256);
 end
  figure(4);
montage(test(:,:,1,1:2:18),gray(256));
set(gcf,'position',[610 100 500 500]);

montage immagine originale, con passo 2 (vengono visualizzate solo slices dispari)

montage immagine interpolata con metodo nearest

montage immagine interpolata con metodo linear

montage immagine interpolata con metodo cubic

PARTE 2. Sottocampionamento con z/2

  % remember D is the original matrix containing the Brain MRI

    P_parte_2 = zeros(256,256,z_new);

    P_parte_2  = (D(:,:,1:2:end-1) + D(:,:,2:2:end))./2;

    [x_int_2,y_int_2,z_int_2] = meshgrid(1:1:256, 1:1:256, 0.75:0.5:z_new+0.25);

  % Calcolo delle interpolazioni con i vari metodi

   D_int_nearest_2 = uint16(interp3(P_parte_2,x_int_2,y_int_2,z_int_2,'nearest'));
   D_int_linear_2 = uint16(interp3(P_parte_2,x_int_2,y_int_2,z_int_2,'linear'));
   D_int_cubic_2 = uint16(interp3(P_parte_2,x_int_2,y_int_2,z_int_2,'cubic'));

   D_int_nearest_2 = double(D_int_nearest_2);
   D_int_linear_2 = double(D_int_linear_2);
   D_int_cubic_2 = double(D_int_cubic_2);

PARTE 2: Calcolo del PSNR nel caso di interpolazione con metodo nearest

   MSE_nearest_temp = (D(:) - D_int_nearest_2(:));
   MSE_nearest_temp = (MSE_nearest_temp).^2;
   MSE_nearest_temp = (MSE_nearest_temp)./(x*y*z);
   MSE_nearest_2 = sum(MSE_nearest_temp);

   PSNR_nearest_2 = 10 * log10((max(D(:))^2)/MSE_nearest_2);

   clear MSE_nearest_temp;

PARTE 2: Calcolo del PSNR nel caso di interpolazione con metodo linear

  D_int_linear_2 = D_int_linear_2 (:,:,2:1:end-1);
  MSE_linear_temp = (D_temp(:) - D_int_linear_2(:));
  MSE_linear_temp = (MSE_linear_temp).^2;
  MSE_linear_temp = (MSE_linear_temp)./(x*y*(z-2));
  MSE_linear_2 = sum(MSE_linear_temp);

   PSNR_linear_2 = 10 * log10((max(D(:))^2)/MSE_linear_2);

   clear MSE_linear_temp;

PARTE 2: Calcolo del PSNR nel caso di interpolazione con metodo cubic

  D_int_cubic_2 = D_int_cubic_2(:,:,2:1:end-1);
  MSE_cubic_temp = (D_temp(:) - D_int_cubic_2(:));
  MSE_cubic_temp = (MSE_cubic_temp).^2;
  MSE_cubic_temp = (MSE_cubic_temp)./(x*y*(z-2));
  MSE_cubic_2 = sum(MSE_cubic_temp);

   PSNR_cubic_2 = 10 * log10((max(D(:))^2)/MSE_cubic_2);

   clear MSE_cubis_temp;

PARTE 2: montage immagine interpolata con metodo nearest

D_int_nearest_2 = uint16(D_int_nearest_2);

for i = 1:size(D_int_nearest_2,3)
    test(:,:,1,i) = gray2ind(imadjust(D_int_nearest_2(:,:,i)),256);
end
figure(5);
montage(test(:,:,1,1:2:20),gray(256));
set(gcf,'position',[10 200 500 500]);

% montage immagine  interpolata con metodo linear

D_int_linear_2 = uint16(D_int_linear_2);

for i = 1:size(D_int_linear_2,3)
    test(:,:,1,i) = gray2ind(imadjust(D_int_linear_2(:,:,i)),256);
end
figure(6);
montage(test(:,:,1,1:2:18),gray(256));
set(gcf,'position',[210 300 400 400]);

% montage immagine  interpolata con metodo cubic

D_int_cubic_2 = uint16(D_int_cubic_2);

for i = 1:size(D_int_cubic_2,3)
    test(:,:,1,i) = gray2ind(imadjust(D_int_cubic_2(:,:,i)),256);
end
figure(7);
montage(test(:,:,1,1:2:18),gray(256));
set(gcf,'position',[410 300 400 400]);

montage immagine interpolata con metodo nearest:

montage immagine interpolata con metodo linear:

montage immagine interpolata con metodo cubic :

Parte 3:Risultati ottenuti:

Published with MATLAB® 7.8