; $Id: read_dicom.pro,v 1.12 2004/01/21 15:55:00 scottm Exp $
;
; Copyright (c) 1993-2004, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.

;
;       Return the first datum reference in iref before the reference imgid.
;
FUNCTION ReadDicomGetFirstBefore, oDicom, iref, imgid

	COMPILE_OPT hidden, strictarr

	iIndex = iref[0]

	parent = oDicom->GetParent(imgid)
	par_iref = oDicom->GetParent(iref)
	v = WHERE(par_iref EQ parent)
	IF (v[0] NE -1) THEN BEGIN
		v_iref = iref[v]
		w = WHERE(v_iref LT imgid)
		IF (w[0] NE -1) THEN BEGIN
			iIndex = v_iref[w[N_ELEMENTS(w)-1]]
		END
	END

	p = oDicom->GetValue(REFERENCE=iIndex,/NO_COPY)

	RETURN, p[0]
END

FUNCTION READ_DICOM, file, red, green, blue, IMAGE_INDEX = iIndex
;+
; NAME:
;	READ_DICOM
;
; PURPOSE:
; 	This function reads an image from a DICOM format file using
;	the IDLffDICOM object interface.
;
; CATEGORY:
;   	Input/Output
;
; CALLING SEQUENCE:
;   	Result = READ_DICOM(File)
;
; INPUTS:
; 	File: The full path name of the file to read.
;       [Red, Green, Blue] = Vectors which return the color palette
;                            (if any)
;
; OPTIONAL KEYWORDS:
;       IMAGE_INDEX - Set this keyword to the index of the image to
;                     read from the file.
;
; OUTPUTS:
;	This function returns a 2D array containing the image data
;	from the file.
;
; KEYWORDS:
;	None
;
; SIDE EFFECTS:
;   	IO is performed.
;
; RESTRICTIONS:
;       Only uncompressed data format is supported (as per current DICOM obj).
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;   RJF, RSI.   Sep, 1998. Original version.
;   RJF, RSI.   Jan, 1999. Filter searches by sequence value.
;-
COMPILE_OPT strictarr

; Set up error handling
CATCH, errorStatus
IF (errorStatus NE 0) THEN BEGIN
	CATCH,/CANCEL
	PRINT, !ERROR_STATE.msg
        IF (OBJ_VALID(oDicom)) THEN OBJ_DESTROY, oDicom
        RETURN, 0L
ENDIF

IF (N_ELEMENTS(iIndex) EQ 0) THEN iIndex=0
IF (iIndex LT 0) THEN iIndex = 0

; Create the DICOM object
oDicom = OBJ_NEW('IDLffDICOM')
IF (NOT OBJ_VALID(oDicom)) THEN BEGIN
	MESSAGE,'IDLffDICOM object not supported on this platform.'
END

; Open the file
IF (oDicom->Read(file) NE 1) THEN BEGIN
	MESSAGE,'The file '+file+' is not in a supported DICOM format.'
END

; Get a list of the images
ref = oDicom->GetReference('7FE0'x,'0010'x)
IF (SIZE(ref,/N_DIMENSIONS) EQ 0) THEN BEGIN
	MESSAGE,'No images could be found in the DICOM file.'
END
IF (iIndex GE N_ELEMENTS(ref)) THEN BEGIN
	MESSAGE,'There are only '+STRING(N_ELEMENTS(ref))+ $
		' images in the DICOM file.'
END

; Get the image in question
img = oDicom->GetValue(REFERENCE=ref[iIndex],/NO_COPY)
ret = *(img[0])
;
; Get additional image info (samples/pixel and palette)
;
iSamp = 1
iref = oDicom->GetReference('0028'x,'0002'x)
IF (SIZE(iref,/N_DIMENSIONS) NE 0) THEN BEGIN
        p = ReadDicomGetFirstBefore(oDicom,iref,ref[iIndex])
        IF (PTR_VALID(p)) THEN iSamp = *p
END

iPlanar = 0
IF (iSamp GT 1) THEN BEGIN
        iref = oDicom->GetReference('0028'x,'0006'x)
        IF (SIZE(iref,/N_DIMENSIONS) NE 0) THEN BEGIN
                p = ReadDicomGetFirstBefore(oDicom,iref,ref[iIndex])
                IF (PTR_VALID(p)) THEN iPlanar = *p
        END
END

;
;	This interface always returns pixel interleaved data...
;
IF (iPlanar NE 0) THEN BEGIN
	ret = TRANSPOSE(ret,[2,0,1])  ; planar to pixel interleave
END

iPalette = 0
iref = oDicom->GetReference('0028'x,'0004'x)
IF (SIZE(iref,/N_DIMENSIONS) NE 0) THEN BEGIN
        p = ReadDicomGetFirstBefore(oDicom,iref,ref[iIndex])
        IF ((PTR_VALID(p)) AND (STRPOS(*p,"PALETTE COLOR") NE -1)) THEN BEGIN
;
;	Read the palette vectors
;
	    FOR i=0,2 DO BEGIN
;
;	Get the offset value
;
		iOffset = 0
		iref = oDicom->GetReference('0028'x,'1101'x+i)
		IF (SIZE(iref,/N_DIMENSIONS) NE 0) THEN BEGIN
			p = ReadDicomGetFirstBefore(oDicom,iref,ref[iIndex])
			IF (PTR_VALID(p)) THEN iOffset = (*p)[1]

		END
;
;	Get the actual palette
;
		iref = oDicom->GetReference('0028'x,'1201'x+i)
		IF (SIZE(iref,/N_DIMENSIONS) NE 0) THEN BEGIN
			p = ReadDicomGetFirstBefore(oDicom,iref,ref[iIndex])
			IF (PTR_VALID(p)) THEN BEGIN
				v = REPLICATE((*p)[0],iOffset+N_ELEMENTS(*p))
				v[iOffset:*] = *p
				CASE i OF
					0 : red = v
					1 : green = v
					2 : blue = v
				ENDCASE
			END
		END

	    END
        END
END

; Done!!
OBJ_DESTROY, oDicom

RETURN, ret

END