Matlab - DICOM

From XennisWiki
Jump to: navigation, search

Digital Imaging and Communications in Medicine (DICOM) ist ein offener Standard zur Speicherung und zum Austausch von Informationen im medizinischen Bilddatenmanagement. (Wikipedia)

Grundlagen

grayValueTransformations.m

%-----Grundlagen--------

% Einlesen und Darstellen des Bildes "CTScull_01.dcm"
ctImage = dicomread('CTScull_01.dcm');
figure('Name', 'Grundlagen');
subplot(2,3,1);
imshow(ctImage);
title('Originalbild ctImage');
minGrayCT = min(ctImage(:));
maxGrayCT = max(ctImage(:));
classCT = class(ctImage);
fprintf('CTScull.dcm\n    Typ: %s\n    Grauwertintervall: [%.0f, %.0f]\n\n' , classCT, minGrayCT, maxGrayCT)

% Ausgabe
%   Typ: int16
%    imMin: -1024
%    imMax:  1962
%
% int16 (range: -32,768 to 32,767)
%
% Die Darstellung ist nicht zufriedenstellend,
% denn die Grauwertinformationen, nutzen nur
% einen sehr kleinen Teil ihres Wertebereichs.
% Dynamik: 2986 Max_Dynamik: 65535

% Verbesserung des Kontrastes
%
% imshow(I,[low high])
% 
% "imshow(I,[low high]) displays the grayscale image I, specifying the
% display range for I in [low high]. [...] If you use an empty matrix ([])
% for [low high], imshow uses [min(I(:)) max(I(:))]; that is, the minimum
% value in I is displayed as black, and the maximum value is displayed as
% white."

subplot(2,3,2);
imshow(ctImage,[]);
title('automatische Anpassung ctImage');

%
% Konvertierung des Datentyps (Casting)
%

%Casting zu "uint8"
ctImage_uint = cast(ctImage,'uint8');
% Plot
subplot(2,3,3);
imshow(ctImage_uint);
title('ctImage gecastet zu uint8');


% Einlesen und Darstellen des Bildes "MRIScull_01.dcm"

mriImage = dicomread('MRIScull_01.dcm');

subplot(2,3,4);
imshow(mriImage);
title('Originalbild mriImage');
minGrayMR = min(mriImage(:));           % 0
maxGrayMR = max(mriImage(:));           % 178
classMR = class(mriImage);              % uint16 (range: -32,768 to 32,767)
fprintf('MRIScull.dcm\n    Typ: %s\n    Grauwertintervall: [%.0f, %.0f]\n\n' , classMR, minGrayMR, maxGrayMR)

% Auf dem Bild ist nichts zu erkennen, denn es wird
% nur ein Bruchteil des eigendlichen Wertebereichs 
% genutzt. Dynamik: 179 Max_Dynamik: 65535

subplot(2,3,5);
imshow(mriImage,[]);
title('automatische Anpassung mriImage');

mriImage_uint = cast(mriImage,'uint8');
subplot(2,3,6);
imshow(mriImage_uint);
title('mriImage gecastet zu uint8');

% Grauchwertintervall nach dem Casting
ctImage_uintMin = min( ctImage_uint(:) );
ctImage_uintMax = max( ctImage_uint(:) );
fprintf('CTScull.dcm - Casting\n    Typ: uint8\n    Grauwertintervall: [%.0f, %.0f]\n\n' , ctImage_uintMin, ctImage_uintMax)

mriImage_uintMin = min( mriImage_uint(:) );
mriImage_uintMax = max( mriImage_uint(:) );
fprintf('MRIScull.dcm - Casting\n    Typ: uint8\n    Grauwertintervall: [%.0f, %.0f]\n\n' , mriImage_uintMin, mriImage_uintMax)

% Das MR-Bild ist nach dem Casten wesendlich besser,
% als das gecastete CT-Bild, denn es hat ursprünglich
% einen Wertebereich, der sich gut mit 8 Bit darstellen
% lässt. Wogegen das CT-Bild einen viel zu großen Bereich
% abdeckt. Dadurch gehen Informationen verloren.
%
% Grauwertintervall CT-Originalbild: (-1024 bis 1962)
% Grauwertintervall CT-cast:         (0 bis 255)
% Grauwertintervall MR-Originalbild: (0 bis 179)
% Grauwertintervall MR-cast:         (0 bis 179)

Window-Leveling

Fortsetzung grayValueTransformations.m

%------Window-Leveling---------

% "Knochenfenster": width = 1500, level = 300
% "Lungenfenster": width = 2000, level = -200
% "Weichteilfenster": width = 350, level = 50
% "Gehirnfenster": width = 100, level = 35

% Vektoren für die Parameter der verschiedenen Fenster
width = [1500,2000,350,100];
level = [300 -200 50 35];

figure('Name', 'Window-Leveling');
for i=1 : 4
   gMin = level(i)-width(i)/2; % Minimaler Grauwert
   gMax = level(i)+width(i)/2; % Maximaler Grauwert
   subplot(2,2,i);
   imshow(ctImage,[gMin gMax]);
   title(sprintf('Width: %d Level: %d',width(i),level(i)));
end

Histogrammdarstellung und Histogrammausgleich

Fortsetzung grayValueTransformations.m

%---Histogrammdarstellung-und-Histogrammausgleich---

%
% Histogramm berechnen und darstellen
%
figure('Name', 'Histogrammdarstellung und -ausgleich');
subplot(2,2,1);
hist_uint = imhist(mriImage_uint);
bar(hist_uint);
title('Histogramm MR Image casted');

%
% Kontrastverbesserung mittels Histogramm?äquilisation
%
subplot(2,2,2);
mriImage_Eq = histeq(mriImage_uint); % Histogrammausgleich
hist_uint = imhist(mriImage_Eq);
bar(hist_uint);
title('Histogramm MR Image casted + Äquilisiert');
subplot(2,2,3);
imshow(mriImage_uint);
title('mriImage gecastet zu uint8');
subplot(2,2,4);
imshow(mriImage_Eq);
title('mriImage casted + Äquilisiert');

Histogrammhyperbolisation

Fortsetzung grayValueTransformations.m

%-------Histogrammhyperbolisation---------

alpha = [0, -1/3, -0.5, -2/3, -0.8];
figure('Name', 'Histogrammhyperbolisation');
for i=1 : length(alpha)
   mriImage_hyp = histoHyperbuint8(mriImage_uint, alpha(i));
   subplot(2,5,i);
   imshow(mriImage_hyp);
   title(sprintf('Alpha = %.2f',(alpha(i))));
   mriImage_hypName = sprintf('MRIScull_hyperb-%03d.png', abs( round( alpha(i) * 100 ) ) );
   imwrite(mriImage_hyp, mriImage_hypName, 'PNG');
   subplot(2,5,i+length(alpha));
   hist = imhist(mriImage_hyp);
   bar(hist);
   title( sprintf('imhist (alpha %.2f)', alpha(i)) )
end

% Aufg. 5
% Unserer Meinung fuehrt  -2/3 <= alpha <= -0.5 zu dem besten (subjektiven)
% Ergebnis.

histoHyperbuint8.m

function [ outImage ] = histoHyperbuint8( inImage, alpha )
% inImage: zu transformierende Quellbild des Datentyps uint8
% alpha: Parameter im Intervall (-1, 0]
% outImage: transformiertes Bild

	% Variablen initialisieren + Hilfsvariablen
	outImage = [];
	allowedType = 'uint8';
	typeMin = 0.0;
	typeMax = 255.0;

	% Falsche und fehlende Eingaben abfangen
	if( nargin < 1 )
		fprintf('Not enough input arguments!\n');
		return;
	end
	if( ~isa(inImage, allowedType) )
		fprintf('Input has to be of type %s!\n', allowedType);
		return;
	end
	if( alpha > 0 || alpha <= -1 )
		fprintf('-1 < alpha <= 0');
		return;
	end

	% Konvertierung zu 1D Vektor
	internalImage = inImage(:);
	imageSize = length(internalImage);

	% Lookup-Tabelle (LUT) initialisieren und füllen
	lut = zeros(2, typeMax + 1);
	lut(1, :) = typeMin : typeMax;

	% Berechnung der absoluten Haeufigkeit
	frequency = zeros( size(lut, 2), 1 );
	for idx = 1 : imageSize
		value = internalImage(idx);
		frequency(value + 1) = frequency(value + 1) + 1;
	end

	% Berechnung der relativen Haeufigkeit
	frequency = cumsum(frequency) / imageSize;

	% Berechnung der transformierten Graufwerte
	lut(2, :) = round( typeMax * frequency.^(1 / (alpha + 1) ) );

	% Speicherreservierung
	outImage = zeros( size(inImage), allowedType );

	% Aufg. 2 (g)
	for idx = 1 : imageSize
		value = internalImage(idx);
		internalImage(idx) = lut(2, value + 1);
	end
	outImage(1 : imageSize) = internalImage;

end