|
FUNCTION Z, T, NT1, TR, SE
B1 = fltARR(10)
C1 = fltARR(10)
Z = 0.
For M = 0,NT1 - 1 do
begin
EE = EXP(-TR(M)/T)
B1(M) = TR(M)*EE
C1(M) = 1. - EE
endfor
For M = 0, NT1 - 1 do begin
For N = 0, NT1 -1 do begin
Z = Z+B1(N)*C1(M)*(C1(N)*SE(M)-C1(M)*SE(N))
endfor
endfor
Return, Z
END
;***************************************************
FUNCTION Z1, T, NT2, TE, SE
C1 = fltARR(10)
B1 = fltARR(10)
D1 = fltARR(10)
Z1 = 0.
For M = 0,NT2 - 1 do begin
B1(M) = EXP(-TE(M)/T)
C1(M) = TE(M) * B1(M)
D1(M) = B1(M) * B1(M)
endfor
For M = 0, NT2 - 1 do begin
For N = 0, NT2 -1 do begin
Z1 = Z1+SE(M)*B1(M)*TE(N)*D1(N)-$
D1(M)*SE(N)*C1(N)
endfor
endfor
Return, Z1
END
;****************************************************
pro t1t2rhov5
RHS = 1.
T1MIN = 0.11
T1MAX = 20.0
T2MIN = 0.011
T2MAX = 100.
STT = fltarr(65536)
ST = fltarr(65536)
IRHO = fltarr(256,256)
SE = fltarr(10)
IJK = fltarr(256)
;************* Open DAT file and read in Data ******
ff=dialog_pickfile(/read, filter="*.DAT")
mytemp=CREATE_STRUCT('Version',float(1),'datastart',$
long(0),'delimiter',byte(44),'missingvalue',$
!Values.F_NaN,'commentsymbol','','fieldcount',$
long(2),'fieldtypes',[7,7],'fieldnames',$
['field1','field2'],'fieldlocations',$
long([0,2]),'fieldgroups',[0,1])
data = READ_ASCII(ff, TEMPLATE=mytemp)
print, data.field1
print, data.field2
constantTR=data.field2[0]
ITRC = constantTR
NT2=data.field1[0] ;number
of t2 images
t2images = data.field2[1:NT2]
T2outputfile = data.field1[NT2+1]
NT1 = data.field1[NT2+2] ;number
of t1 images
ConstantTE = data.field2[NT2+2]
ITEC = ConstantTE
t1images = data.field2[NT2+3:*]
xy = fix(NT2)
xx = fix(NT1)
T1outputfile = data.field1[xx+xy+3]
RHOoutputfile = data.field1[xx+xy+4]
SET2 = fltarr(xy,256,256)
SET1 = fltarr(xx,256,256)
cd, dialog_pickfile(/directory)
;change to directory with the MRI data images
;***************opening t2images*****************
swap = DIALOG_MESSAGE('Do you want to swap bytes?'$
,question=1)
for i=0,xy-1 do begin
openr, unit, t2images[i] ,/get_lun
image2 = intarr(256,256)
readu, unit, image2
close, unit
free_lun, unit
SET2[i,*,*] = image2[*,*]
if swap eq 'Yes' then $
SET2[i,*,*] = swap_endian(image2[*,*])
endfor
window, 2, xsize = 256, ysize = 256
tvscl, SET2[XY-1,*,*]
;**************opening t1images***********************
for i=0,xx-1 do begin
openr, unit, t1images[i] ,/get_lun
image1 = intarr(256,256)
readu, unit, image1
close, unit
free_lun, unit
SET1[i,*,*] = image1[*,*]
if swap eq 'Yes' then SET1[i,*,*]
= $
swap_endian(image1[*,*])
endfor
window, 1, xsize = 256, ysize = 256
tvscl, SET1[XX-1,*,*]
;****Calculate MIN Signal Cut Off for T2 & T1 ********
;Using BOX_CURSOR
;Once the box cursor has been realized, hold down the
;left mouse button to move the box by dragging.
;Hold down the middle mouse button to resize the box
;by dragging. (The corner nearest the initial mouse
;position is moved.)
;Press the right mouse button to exit the procedure
;and return the current box parameters.
box_cursor, llx, lly, xsize, ysize
llx=fix(llx)
lly=fix(lly)
xs=fix(xsize)
ys=fix(ysize)
window, title = 'histogram of background'
plot,histogram(SET1[*,llx:llx+xs,lly:lly+ys]+$
SET2[*,llx:llx+xs,lly:lly+ys],
min = 0, $
max =
max(SET1[*,llx:llx+xs,lly:lly+ys]+$
SET2[*,llx:llx+xs,lly:lly+ys]) ),
$
ytitle = 'count', xtitle = 'pixel
value'
scot1 = max(SET1[*,llx:llx+xs,lly:lly+ys]) *NT1
;scot1 t1 signal cut off * number of t1 images
scot2 = max(SET2[*,llx:llx+xs,lly:lly+ys]) *NT2
;scot2 t2 signal cut off * number of t2 images
print, 'scot1 = ' ,scot1
print, 'scot2 = ' ,scot2
TE=data.field1[1:xy]/1000.
TR=data.field1[xy+3:xy+2+xx]/1000.
set2 = reform(set2,nt2,n_elements(set2)/nt2)
set1 = reform(set1,nt1,n_elements(set1)/nt1)
;*******Calculate T2 DATA Matrix *********************
IF(nt2 EQ 0) then goto, jumpt1
IF(nt2 EQ 1) then begin
stt=set2[0,*]
goto,jumpt1
endif
TimeT2 = SYSTIME(1)
print, 'calculating T2'
For I=0L,65535 do begin
TMIN = T2MIN
TMAX = T2MAX
Sum = 0.
se = set2[*,i]
SUM = total(se)
A = Z1(TMIN, NT2, TE, SE)
B = Z1(TMAX, NT2, TE, SE)
IF((A*B GE 0.) or (SUM LE SCOT2))
then begin
T2 = 0.
endif else begin
jumphere:
T = (TMIN*B-TMAX*A)/(B-A)
C = Z1(T, NT2,TE, SE)
IF((A*C LT 0.)AND(ABS(C) GT
0.0001))$
then begin
B = C
TMAX = T
goto, jumphere
endif else begin
IF((B*C LT 0.) AND $
(ABS(C) GT
0.0001)) $
then
begin
A = C
TMIN = T
goto, jumphere
endif else begin
T2 = T
endelse
endelse
endelse
STT(I) = T2
print, I
endfor
print, 'Finished calculating T2.'
PRINT, SYSTIME(1) - TimeT2, ' Seconds'
TimeT2 = SYSTIME(1) - TimeT2
;*****storing T2 data ********************************
stt2=stt*10000
stt3=long(stt2)
stt3 = reform(stt3,256,256)
WRITE_TIFF,T2outputfile,stt3,/LONG
;***************Calculate T1*************************
jumpt1:
IF (NT1 EQ 0) then goto, jumpend
print, 'calculating T1'
TimeT1 = SYSTIME(1)
For I=0L,65535 do begin
TMIN = T1MIN
TMAX = T1MAX
Sum = 0.
se = set1[*,i]
SUM = total(se)
A = Z(TMIN, NT1, TR, SE)
B = Z(TMAX, NT1, TR, SE)
IF((A*B GE 0.) or (SUM LE SCOT1))
then begin
T1 = 0.
endif else begin
jumphere1:
T = (TMIN*B-TMAX*A)/(B-A)
C = Z(T, NT1, TR, SE)
IF((A*C LT 0.)AND(ABS(C) GT
0.0004))$
then begin
B = C
TMAX
= T
goto, jumphere1
endif else begin
IF((B*C LT 0.) AND $
(ABS(C) GT 0.0004)) $
then begin
A = C
TMIN = T
goto, jumphere1
endif else begin
T1 = T
endelse
endelse
endelse
ST(I) = T1
print, I
endfor
print, 'Finished calculating T1.'
PRINT, SYSTIME(1) - TimeT1, ' Seconds'
TimeT1 = SYSTIME(1) - TimeT1
;*****storing T1 data ********************************
st2=st*1000
st3=fix(st2)
openw, lun, T1outputfile, /GET_LUN
writeu, lun, st3
close, lun
free_lun, lun
IF(NT2 EQ 0) then goto, jumpend
;**************Calculate RHO *************************
jumprho:
print, 'calculating rho'
TEC = Float(ITEC)/1000. ;itec
is 15
I = 0L
For KI = 0, 255 do begin
For KII = 0, 255 do begin
I = I+1
IF (I EQ 65536) then goto, endrho
IF((ST(I) LE 0.001) OR (STT(I) LE
0.001)) $
then
begin
IRHO(KI,KII) = 0
goto, endrho
endif
SY=0.
SYY=0.
For KKK = 0, NT1-1 do begin
PAR = 1. -
EXP(-TR(KKK)/ST(I))
SY = Temporary(SY) +
set1[KKK,I]*PAR
SYY
= Temporary(SYY)+PAR*PAR
endfor
IRHO(KI,KII) = fix(((SY/SYY)/$
EXP(-TEC/STT(I)))*RHS)
endrho:
endfor
endfor
IRHO = fix(temporary(IRHO))
print, 'Finished calculating rho.'
irho = rotate(irho,4)
;***********storing rho ******************************
openw, lun, RHOoutputfile, /GET_LUN
writeu, lun, irho
close, lun
free_lun, lun
;*****************************************************
jumpend:
end
|