{
Dicom.p
  by: Jim Nash, Synergistic Research Systems (jim.nash@his.com)
  Reads and decodes the DICOM header so that NIH Image can
  import DICOM images. DICOM (Digital Imaging and Communications
  in Medicine) is a format popular in the medical imaging
  community. This code is in the public domain.
}


unit DICOM;

interface

	uses
		Types, Memory, QuickDraw, Packages, Menus, Events, Fonts, Scrap, ToolUtils, Resources, 
		Errors, Palettes, Printing, StandardFile, Folders, TextUtils, Files,
		globals, Utilities, Text, Graphics, Utilities, file2;


	procedure ImportDICOMImages (fname: Str255; RefNum: integer; ImportAll: boolean; {}
									function ImportFile (fname: str255; vnum: integer): boolean);

implementation

	const
		dDicomName = 'DICOM dictionary';
		maxElements = 1000;
		elemNameLength = 50;
	type
		DataKind = (kUnknown, kString, kInteger, kLongint, kReal, kUInteger, kULongint);

		DataElement = record
				group, element: integer;
				code: packed array[1..2] of char;
				list: boolean;
				name: string[elemNameLength];
			end;
		ElemArray = array[1..maxElements] of DataElement;
		ElemArrayPtr = ^ElemArray;

		DataDictionary = record
				number: integer;
				elem: ElemArrayPtr;
			end;
		DataDictionaryPtr = ^DataDictionary;
	var
		dictionary: DataDictionaryPtr;
		loaded: boolean;
		mySliceSpacing: real;


{ **************  Utility routines ***************** }

	procedure StringToBase (s: Str255; base: integer; var value: longint);
{converts a string in some base to longint.  Typically}
{base = 2,8,10,16 to represent binary, octal, decimal and hexadecimal}
		var
			ch: char;
			good: boolean;
			len, digit: integer;
			i: longint;
	begin
		i := 1;
		value := 0;
		len := length(s);
		while (i <= len) do begin
			good := true;
			ch := s[i];
			if ch in ['A'..'Z'] then
				digit := ord(ch) - ord('A') + 10
			else if ch in ['a'..'z'] then
				digit := ord(ch) - ord('a') + 10
			else if ch in ['0'..'9'] then
				digit := ord(ch) - ord('0')
			else
				good := false;
			if good then
				value := value * base + digit;
			i := i + 1;
		end;
	end;


	procedure BaseToString (value: longint; base: integer; var s: Str255);
{converts a long integer to a string in any base.  Typically}
{base = 2,8,10,16 to represent binary, octal, decimal and hexadecimal.}
{Ignores the sign bit unless base=10.}
		var
			sign, decimal: boolean;
			digit: integer;
			ch: char;
	begin
		decimal := (base = 10);
		s := '';
		sign := (value < 0);
		if decimal then
			value := abs(value)
		else
			value := BAND(value, $7FFFFFFF);
		if value = 0 then
			s := '0'
		else
			while (value <> 0) do begin
				digit := value mod base;
				value := value div base;
				if (digit >= 0) and (digit <= 9) then
					ch := chr(digit + ord('0'))
				else
					ch := chr(digit - 10 + ord('A'));
				s := concat(ch, s);
			end;
		if sign then
			if decimal then
				s := concat('-', s)
			else if s[1] < '2' then
				s[1] := chr(ord(s[1]) + 8)
			else
				s[1] := chr(ord(s[1]) - ord('2') + ord('A'));
	end;


	function htos (i: longint): Str255;
{A convenience function to replace BaseToString (hexadecimal) }
		var
			s: Str255;
	begin
		BaseToString(i, 16, s);
		htos := s;
	end;


	function itos (i: longint): Str255;
{A convenience function to replace NumToString}
		var
			s: Str255;
	begin
		NumToString(i, s);
		itos := s;
	end;


{ **************  DICOM routines ***************** }


	procedure InitDICOM;
		var
			err: integer;
	begin
		dictionary := nil;
		loaded := false;
		DicomInitialized:=true;
	end;



	function FindDicomElement (group, element: integer): integer;
		var
			i, index: integer;
	begin
		index := 0;
		if loaded and (dictionary <> nil) then
			with dictionary^ do
				if (elem <> nil) then begin
					i := 1;
					while (group > elem^[i].group) and (i < number) do
						i := i + 1;
					if (i <= number) then
						while (element > elem^[i].element) and (group = elem^[i].group) and (i < number) do
							i := i + 1;
					if (i <= number) then
						if (element = elem^[i].element) and (group = elem^[i].group) then
							index := i;
				end;
		FindDicomElement := index;
	end;


	procedure InitUserChoice;
{selected elements to list in text window, short form.}
{This a minimum set.  If a user wants more information, they do a full dump.}

		procedure Select (group, element: integer);
			var
				index: integer;
		begin
			with dictionary^ do begin
				index := FindDicomElement(group, element);
				if index > 0 then
					elem^[index].list := true;
			end;
		end;

	begin
		with dictionary^ do begin
			Select($8, $20);
			Select($8, $30);
			Select($8, $60);
			Select($8, $1030);
			Select($8, $103E);
			Select($8, $1070);

			Select($10, $10);
			Select($10, $20);
			Select($10, $21B0);

			Select($18, $10);
			Select($18, $50);
			Select($18, $88);

			Select($20, $10);
			Select($20, $11);
			Select($20, $12);
			Select($20, $13);

			Select($28, $10);
			Select($28, $11);
			Select($28, $30);
			Select($28, $100);
		end;
	end;


	procedure LoadDataDictionary;
		type
			CharBuf = packed array[0..100000] of char;
			CharBufPtr = ^CharBuf;
		var
			err, refnum, len, i1, i2, n: integer;
			index1, index2, logEOF, count, num, theSize: longint;
			f: text;
			sp: StringPtr;
			str: Str255;
			s1: Str255;
			buf: CharBufPtr;

	begin
		if dictionary = nil then begin
			dictionary := DataDictionaryPtr(NewPtr(sizeof(DataDictionary)));
			if dictionary <> nil then
				dictionary^.elem := ElemArrayPtr(NewPtr(sizeof(ElemArray)));
			loaded := false;
		end;
		if (not loaded) and (dictionary <> nil) then
			with dictionary^ do begin
				err := HSetVol(nil, StartupSpec.vRefNum, StartupSpec.parID);
				err := FSOpen(dDicomName, 0, refnum); 	{check that file is present}
				if (err = 0) and (elem <> nil) then begin
					err := GetEOF(refnum, logEOF);
					buf := CharBufPtr(NewPtr(logEOF + 10));
					if (buf <> nil) then begin
						loaded := true;
						number := 0;
						count := logEOF;
						err := FSRead(refnum, count, ptr(buf));
						err := FSClose(refnum);
						index1 := 0;
						repeat
							index2 := index1;
							str := '';
							while (buf^[index2] <> cr) and (index2 < logEOF) and (length(str) < 255) do begin
								str := concat(str, buf^[index2]);
								index2 := index2 + 1;
							end;
							index1 := index2 + 1;
							len := length(str);
							if len > 0 then
								if str[1] = '{' then begin
									number := number + 1;
									if (number mod 10) = 0 then
										ShowAnimatedWatch;
									with elem^[number] do begin
										list := false;

										i1 := pos('x', str);
										s1 := copy(str, i1 + 1, 4);
										StringToBase(s1, 16, num);
										group := num;
										str := copy(str, i1 + 6, length(str)-(i1 + 6));

										i1 := pos('x', str);
										s1 := copy(str, i1 + 1, 4);
										StringToBase(s1, 16, num);
										element := num;
										str := copy(str, i1 + 6, length(str)-(i1 + 6));

										i1 := pos('''', str);
										if length(str) >= (i1 + 2) then begin
											code[1] := str[i1 + 1];
											code[2] := str[i1 + 2];
											str := copy(str, i1 + 5, length(str)-(i1 + 5));
										end
										else
											str := '';

										i1 := pos('"', str);
										if i1 > 0 then
											str[i1] := ' ';
										i2 := pos('"', str);
										if i2 = 0 then
											number := number - 1
										else begin
											n := i2 - i1 - 1;
											if n > elemNameLength then
												n := elemNameLength;
											name := copy(str, i1 + 1, n);
										end;
									end;
								end;
						until (index1 >= logEOF);
					end;
					DisposePtr(ptr(buf));
				end;
				InitUserChoice;
			end;
	end;
{$R+}


	function GetDataKind (index: integer): DataKind;
		var
			kind: DataKind;
	begin
		kind := kUnknown;
		if (dictionary <> nil) and (index > 0) then
			with dictionary^.elem^[index] do begin
				if (code = 'AE') or (code = 'AS') or (code = 'CS') or (code = 'DA') or (code = 'DS') then
					kind := kString
				else if (code = 'DT') or (code = 'IS') or (code = 'LO') or (code = 'LT') or (code = 'PN') then
					kind := kString
				else if (code = 'SH') or (code = 'ST') or (code = 'TM') or (code = 'UI') then
					kind := kString
				else if (code = 'SS') then
					kind := kInteger
				else if (code = 'SL') then
					kind := kLongint
				else if (code = 'US') then
					kind := kUInteger
				else if (code = 'UL') then
					kind := kULongint;
			end;
		GetDataKind := kind;
	end;

	procedure ImportDICOMImages (fname: Str255; RefNum: integer; ImportAll: boolean; {}
									function ImportFile (fname: str255; vnum: integer): boolean);
		var
			enable_text, enable_open_text, first_image, sw, listAll, UseFixedScale: boolean;
			ImageNumber:integer;
			myIntercept, myScale: extended;

		function GetDICOMParams (fname: Str255; vNum: integer): integer;
			const
				id_offset = 128;			{location of "DICM"}
				firstDicomElement = 132;	{first element}
				maxbuf = 20000;
			type
				name4 = packed array[1..4] of char;
				name4ptr = ^name4;
			var
				open_sw, done, window_sw: boolean;
				f, err, index, len: integer;
				groupWord, elementWord, lastGroup, FirstElement: integer;
				height, width, bits_alloc, bits_stored, high_bit, representation, offset, bitsAllocated: integer;
				seriesMin, seriesMax: integer;           {intensity range}
				scale, aspect, units: Str255;       {spatial}
				s, imgNumString, sliceSpacingStg, rescaleInterceptStg, rescaleSlopeStg: Str255;
				buflen, elementLength, groupLength: longint;
				buf: packed array[0..maxbuf] of byte;
				kind: DataKind;
				dictionaryIndex: integer;
				xStr,yStr, vr:str255;

			procedure MyWriteElement (str: Str255);
				const
					spaces = '                                                                                               ';
					padWidth = 4;
					nameWidth = 32;
				var
					s1, s2: str255;

				function pad (s: Str255): Str255;
					const
						width = 4;
				begin
					while length(s) < width do
						s := concat(' ', s);
					pad := s;
				end;

			begin
				with dictionary^.elem^[dictionaryIndex] do begin
					s2 := name;
					if listAll then begin
						s1 := concat('(', pad(htos(groupWord)), ',', pad(htos(elementWord)), ')  (', pad(itos(elementLength)), ')');
						s2 := copy(concat(name, spaces), 1, nameWidth);
						s2 := concat(s1, '  ', code[1], code[2], '  ', s2);
					end;
					str := concat(s2, ':  ', str);
					if enable_text then begin
						if groupWord <> lastGroup then
							InsertText('', true);
						lastGroup := groupWord;
						InsertText(str, true);
					end;
				end;
			end;

			function GetInteger (index: integer): integer;
				var
					i: integer;
			begin
				i := buf[index] + $100 * buf[index + 1];
				GetInteger := i;
			end;

			function GetLongint (index: integer): longint;
				var
					i: longint;
			begin
				i := Ord4(buf[index]) + $100 * (buf[index + 1] + $100 * (buf[index + 2] + $100 * buf[index + 3]));
				GetLongint := i;
			end;

			function GetLength (i: integer): longint;
			begin
				vr[1] := chr(buf[i]);
				vr[2] := chr(buf[i + 1]);
				if vr[2] < 'A' then {implecit vr with 32-bit length}
					GetLength := Ord4(buf[i]) + $100 * (buf[i+1] + $100 * (buf[i+2] + $100 * buf[i+3]))
				else if (vr = 'OB') or (vr = 'OW') or (vr = 'SQ') then begin  {explicit VR with 32-bit length}
				    i := i + 4;  {skip 2 byte string and 2 reserved bytes}
				    index := index + 4;
					GetLength := Ord4(buf[i]) + $100 * (buf[i+1] + $100 * (buf[i+2] + $100 * buf[i+3]))
				end else  {explicit VR with 16-bit length}
					GetLength := Ord4(buf[i+2]) + $100 * (buf[i+3]);
			end;

			function GetUInteger (index: integer): longint;
				var
					i: integer;
			begin
				i := buf[index] + $100 * buf[index + 1];
				GetUInteger := BAND(i, $FFFF);
			end;

			function GetULongint (index: integer): longint;
			{does not correctly report numbers > $7FFFFFFF}
			begin
				GetULongint := GetLongint(index);
			end;

			function GetString (index: integer): Str255;
				var
					i: integer;
					s: Str255;
			begin
				s := '';
				for i := 0 to elementLength - 1 do
					s := concat(s, chr(buf[index + i]));
				GetString := s;
			end;

		
            function htos2(i: LongInt): str255;
            {Converts an integer to hex using a fixed field width of 6}
            var
            	s: str255;
            begin
	            s := htos(i);
	            while length(s) < 6 do
	            	s := concat(' ', s);
	            htos2 := s;
            end;


			procedure DoByteSwap (var i: LongInt);
				var
					a: ostype;
					c: char;
			begin
				a := ostype(i);
				c := a[1];
				a[1] := a[2];
				a[2] := c;
				c := a[3];
				a[3] := a[4];
				a[4] := c;
				i := LongInt(a)
			end;
			
			
	procedure Swap(var i: LongInt);
		var
			a: ostype;
			c: char;
	begin
		a := ostype(i);
		a[1] := a[3];
		a[2] := a[4];
		a[3] := chr(0);
		a[4] := chr(0);
		i := LongInt(a)
	end;


			procedure GetNextElement;
				var
					i: longint;
					str: Str255;
			begin
				if index = 0 then
					index := firstElement
				else
					index := index + elementLength;
				if (index < 0) or (index >= buflen) then
					exit(GetNextElement);
				groupWord := GetInteger(index);
				elementWord := GetInteger(index + 2);
				elementLength := GetLength(index + 4);
				if elementLength = 13 then {hack needed to read some GE files}
					elementLength := 10;
				if ControlKeyDown then
                	InsertText(stringOf(index:6, htos2(groupWord), htos2(elementWord), htos2(elementLength)), true);
                index := index + 8;
				dictionaryIndex := FindDicomElement(groupWord, elementWord);
				if dictionaryIndex > 0 then
					with dictionary^ do
						if elem^[dictionaryIndex].list or listAll then
							with elem^[dictionaryIndex] do begin
								kind := GetDataKind(dictionaryIndex);
								case kind of
									kString: 
										str := GetString(index);
									kInteger: 
										str := itos(GetInteger(index));
									kLongint: 
										str := itos(GetLongint(index));
									kUInteger: 
										str := itos(GetInteger(index));
									kULongint: 
										str := itos(GetULongint(index));
									otherwise
										str := 'unknown format';
								end;
								MyWriteElement(str);
							end;
			end;

			function IsElement (group, element: integer): boolean;
			begin
				IsElement := (group = groupWord) and (element = elementWord);
			end;

			procedure TestError (err1: integer; str: Str255);
				var
					str1: Str255;
			begin
				err := err1;
				if err1 <> 0 then begin
					if err <> 1 then
						str := concat(str, ' - error ', itos(err));
					PutMessage(str);
					if open_sw then
						err1 := fsclose(f);
					GetDICOMParams := err;
					exit(GetDICOMParams);
				end;
			end;

			procedure OpenDicomTextWindow;
				var
					width, height: integer;
			begin
				if listAll then begin
					width := 500;
					height := 400;
				end
				else begin
					width := 350;
					height := 300;
				end;
				if enable_open_text then
					window_sw := MakeNewTextWindow(concat(fname, ' header'), width, height);
				CurrentFontID := monaco;
				CurrentSize := 9;
				ChangeFontOrSize;
				enable_open_text := false;
				if enable_text then begin
					if loaded then
						InsertText('Selected fields from the DICOM file header', true)
					else begin
						InsertText(concat('Can''t find file: ', dDicomName, '.'), true);
						InsertText('', true);
						InsertText('This file is required to decode the DICOM header. It is available from: ftp://zippy.nimh.nih.gov/pub/nih-image/documents/dicom-dict.hqx. It must be located in the same folder as NIH Image or in the System folder.', true)
					end;
					InsertText('', true);
				end;
			end;

		begin
		    vr :='--';
			err := 0;
			buflen := maxbuf + 1;
			open_sw := false;
			TestError(fsopen(fname, vNum, f), 'Open');
			open_sw := true;
			TestError(FSRead(f, buflen, @buf), 'Read');
			TestError(fsclose(f), 'Close');
			if name4ptr(longint(@buf) + id_offset)^ = 'DICM' then
				FirstElement:=FirstDicomElement
			else if name4ptr(longint(@buf))^ = 'DICM' then
				FirstElement:=4
			else
				FirstElement:=0; {TestError(1, 'This is not a DICOM file.');}
			OpenDicomTextWindow;
			sliceSpacingStg := '1.0';
			rescaleInterceptStg := '0.0';
			rescaleSlopeStg := '1.0';
			imgNumString := '';
			height := -1;
			width := -1;
			offset := -1;
			index := 0;
			lastGroup := $8;
			done := false;
			scale := '0.0';
			aspect := '1.0';
			representation := 0;
			bitsAllocated := 16;
			seriesMin := 0;
			seriesMax := 0;
			units := '';
			repeat
				GetNextElement;
				if (index < 0) or (index >= buflen) then
					leave;
				if (elementWord = 0) and (elementLength = 0) then
					leave;

				if IsElement($18, $88) then
					sliceSpacingStg := GetString(index)
				else if IsElement($20, $13) then
					imgNumString := GetString(index)
				else if IsElement($28, $10) then
					height := GetInteger(index)
				else if IsElement($28, $11) then
					width := GetInteger(index)
				else if IsElement($28, $30) then
					scale := GetString(index)
				else if IsElement($28, $34) then
					aspect := GetString(index)
				else if IsElement($28, $100) then
					bitsAllocated := GetInteger(index)
				else if IsElement($28, $103) then
					representation := GetInteger(index)
				else if IsElement($28, $108) then
					seriesMin := GetInteger(index)
				else if IsElement($28, $109) then
					seriesMax := GetInteger(index)
				else if IsElement($28, $1052) then
					rescaleInterceptStg := GetString(index)
				else if IsElement($28, $1053) then
					rescaleSlopeStg := GetString(index)
				else if IsElement($7FE0, $10) then begin
					offset := index;
					done := true;
				end else if IsElement($7F88, $10) then begin
					offset := index;
					done := true;
				end;
				if CommandPeriod then
					listAll := false;
			until done;
			if (width = -1) or (height = -1) or (offset = -1) then begin
				InsertText('', true);
				InsertText('Unable to decode DICOM header.', true);
				GetDICOMParams := -1;
				exit(GetDICOMParams)
			end;

			{Image dimension information}
			ImportCustomWidth := width;
			ImportCustomHeight := height;
			ImportCustomOffset := offset;
			ImportSwapBytes := true; {(representation = 1);}

			{Intensity information}
			if bitsAllocated = 8 then begin
				ImportCustomDepth := EightBits;
				if ImageNumber=1 then
					InsertText('', true);
				GetDICOMParams := err;
				exit(GetDICOMParams);
			end else
				ImportCustomDepth := SixteenBitsSigned;
			if not ((seriesMin = 0) and (seriesMax = 0)) then begin
				ImportAutoScale:=false;
				ImportMin:=seriesMin;
				ImportMax:=seriesMax;
			end else if (ImageNumber>1) and UseFixedScale then begin
				ImportAutoScale:=false;
				ImportMin:=info^.CurrentMin;
				ImportMax:=info^.CurrentMax;
			end else
				ImportAutoScale:=true;

			{convert from scaled units to independent units}
			myIntercept := StringToReal(rescaleInterceptStg);
			myScale := StringToReal(rescaleSlopeStg);

			{Spatial scale information}
			with Info^ do begin
				PixelAspectRatio := StringToReal(aspect);
				xScale := 1;
				yScale := 1;
				zScale := 1.0 / StringToReal(sliceSpacingStg);

				xUnit := '';
				SpatiallyCalibrated := false;
				if scale <> '' then begin
					xStr:=copy(scale, pos('\', scale) + 1, length(scale) - pos('\', scale));
					xScale := StringToReal(xStr);
					yStr:=copy(scale, 1, pos('\', scale) - 1);
					yScale := StringToReal(yStr);
					xUnit := 'mm';
					SpatiallyCalibrated := (xScale <> 0.0) and (yScale <> 0.0);
					if SpatiallyCalibrated then begin
						xScale := 1.0 / xScale;
						yScale := 1.0 / yScale;
					end;
				end;
			end; {with}
			if ImageNumber=1 then
				InsertText('', true);
			GetDICOMParams := err;
		end;


		procedure UpdateCoefficients;
		{Scale coefficients given Dicom Rescale Intercept and Rescale Slope}
		begin
		{
			with info^ do begin
				info^.Coefficient[1] := myIntercept + myScale * info^.Coefficient[1];
				info^.Coefficient[2] := myScale * info^.Coefficient[2];
				fit := StraightLine;
				GenerateValues;
			end;
		}
		end;
		

		procedure ImportAllDicomFiles (RefNum: integer);
			var
				OpenedOK: boolean;
				index: integer;
				name: Str255;
				ftype: OSType;
				err: OSErr;
				PB: HParamBlockRec;
		begin
			index := 0;
			while true do begin
				index := index + 1;
				with PB do begin
					ioCompletion := nil;
					ioNamePtr := @name;
					ioVRefNum := RefNum;
					ioVersNum := 0;
					ioFDirIndex := index;
					err := PBGetFInfoSync(@PB);
					if err = fnfErr then
						exit(ImportAllDicomFiles);
					ftype := ioFlFndrInfo.fdType;
				end;

				if GetDICOMParams(name, RefNum) <> 0 then
					exit(ImportAllDicomFiles);
				WhatToImport := ImportCustom;
				if not ImportFile(name, RefNum) then
					exit(ImportAllDicomFiles);
				WhatToImport := ImportDicom;
				if (myIntercept <> 0.0) or (myScale <> 1.0) then
					UpdateCoefficients;
				with info^ do InsertText(StringOf(ImageNumber:3, ': "', title, '", min=', CurrentMin:1, ', max=', CurrentMax:1), true);
				ImageNumber:=ImageNumber+1;
				enable_text := false;		{text saved for first image only}
				first_image := false;
				if CommandPeriod then begin
					beep;
					exit(ImportAllDicomFiles);
				end;
			end;
		end;

	begin		{ImportDICOMImages}
		if not DicomInitialized then
			InitDICOM;
		listAll := OptionKeyDown or OptionKeyWasDown;
		enable_open_text := true;
		enable_text := true;
		first_image := true;
		ImageNumber:=1;
		UseFixedScale:=ShiftKeyDown;
		LoadDataDictionary;
		ImportingDicom := true;
		if ImportAll then
			ImportAllDicomFiles(RefNum)
		else begin
			if GetDICOMParams(fname, RefNum) <> 0 then
				exit(ImportDICOMImages);
			WhatToImport := ImportCustom;
			if ImportFile(fname, RefNum) then
				if (myIntercept <> 0.0) or (myScale <> 1.0) then
					UpdateCoefficients;
			WhatToImport := ImportDicom;
			with info^ do InsertText(StringOf('file "', title, '", min=', CurrentMin:1, ', max=', CurrentMax:1), true);
		end;
		ImportingDicom := false;
	end;

end.