Pydicomを使ってRDSRを読んで必要な情報を取り出してみる
RDSRはRadiation Dose Structured Reportの略で、CTや血管造影などX線を用いた検査でどのくらX線量（Radiation dose）を照射したのか、構造化された書式（Structured report）にて出力したものです。
RDSRはDICOM規格で規定されています。詳細は規格書を読みましょう。
今回は時間がなく規格書を理解せず、帰納的にSRの解析に取り組みました。
解析結果が正しいかどうか、医療現場で使用しているソフトウェアの結果で検証しています。
検証した範囲では問題なさそうですが、利用は自己責任でお願いします。
プログラムはこちらから
RDSRファイルはこちらから
（この記事は2022年4月24日に日本放射線技術学会 北海道支部 第78回春季学術大会 医療情報セミナーで発表したものです）
PydicomはDICOM画像だけでなくSRも読むことができます。RDSRファイルを読み込んで表示するだけであれば数行で書くことができます。
import pydicom
fnm = "./CT_rdsr.SR"
ds = pydicom.dcmread(fnm)
print(ds)
fnm = "./CT_rdsr.SR"
ds = pydicom.dcmread(fnm)
print(ds)
このたった4行でRDSRのDICOMファイル"rdsr.sr"の内容をpydicomを用いて変数dsに読み込み、それをプリントすることでRDSRの内容を確認することができる。
DICOM RDSRの内容を見てみる
DICOMのファイルはデータセット呼ばれるデータエレメントの集合でできています。
データエレメントは、タグと呼ばれる16進数のグループ番号とエレメント番号から始まり、キーワード、値表現（VR）、値の長さ（VL）と値（VF, Value field）で構成されています。
データエレメントは、タグと呼ばれる16進数のグループ番号とエレメント番号から始まり、キーワード、値表現（VR）、値の長さ（VL）と値（VF, Value field）で構成されています。
print(ds)で出力された内容を見てみしょう。
表示された各行はデータエレメントを表していて次のように表示されています。
(グループ番号, エレメント番号) キーワード VR:値
(0002,0002)と(0008,0016)のSOP Class UIDが”X-Ray Radiation Dose SR Storage"と記述されています。
また、(0008,0060)のModalityは”SR”とあり、このファイルがStructured reportであることがわかります。
Dataset.file_meta -------------------------------
(0002, 0000) File Meta Information Group Length UL: 196
(0002, 0001) File Meta Information Version OB: b'\x00\x01'
(0002, 0002) Media Storage SOP Class UID UI: X-Ray Radiation Dose SR Storage
(0002, 0003) Media Storage SOP Instance UID UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000249
(0002, 0010) Transfer Syntax UID UI: Explicit VR Little Endian
(0002, 0012) Implementation Class UID UI: 1.3.12.2.1107.5.4.7
(0002, 0013) Implementation Version Name SH: 'SIEMENSXLEOVB21C'
-------------------------------------------------
(0008, 0005) Specific Character Set CS: 'ISO_IR 100'
(0008, 0016) SOP Class UID UI: X-Ray Radiation Dose SR Storage
(0008, 0018) SOP Instance UID UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000249
(0008, 0020) Study Date DA: '20220224'
(0008, 0021) Series Date DA: '20220224'
(0008, 0023) Content Date DA: '20220224'
(0008, 0030) Study Time TM: '074907.240000'
(0008, 0031) Series Time TM: '075458.275000'
(0008, 0033) Content Time TM: '075458.315000'
(0008, 0050) Accession Number SH: '6010016637053'
(0008, 0060) Modality CS: 'SR'
上記以外の情報は通常のDICOM画像とそう変わりありません。
しかし、もう少し下のほうを見てみると、(0040,0275)のRequest Attributes Sequence以下で字下げ（インデント）が見られ、
さらに(0040,0008) Scheduled Protocol Code Sequence以下でインデントが発生しています。
これはVRが"SQ"に設定されているためです。
”SQ”は「ある意味のある集合のシーケンス（連続した項目）を含み、このシーケンスは入れ子状態をとることができる」ものです。
詳しくはDICOM Part５（http://www.jira-net.or.jp/dicom/dicom_data_02_01.html）を参照してください。
(0029, 1008) [syngo Report Type] CS: 'CTDOSEREPORT'
(0029, 1009) [syngo Report] LO: '1.0'
(0029, 1015) [SR Variant] US: 6
(0029, 1017) [SC SOP Instance UID] UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000247
(0040, 0275) Request Attributes Sequence 1 item(s) ----
(0040, 0008) Scheduled Protocol Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '3003100550020000'
(0008, 0102) Coding Scheme Designator SH: '99'
(0008, 0103) Coding Scheme Version SH: '6010016637053'
(0008, 0104) Code Meaning LO: 'AG?F?~???ETAE?F?S?g'
---------
(0040, 0009) Scheduled Procedure Step ID SH: '6010016637053'
(0040, 1001) Requested Procedure ID SH: '6010016637053'
---------
この"SQ"をさらに追いかけてみます。
すると、(0040, a730) Content Sequence が20項目あると表示されています。
つまり、これ以下に何らかの意味ある項目が20項目あることを示しています。
(0040, a050) Continuity Of Content CS: 'SEPARATE'
(0040, a372) Performed Procedure Code Sequence 0 item(s) ----
(0040, a491) Completion Flag CS: 'COMPLETE'
(0040, a493) Verification Flag CS: 'UNVERIFIED'
(0040, a504) Content Template Sequence 1 item(s) ----
(0008, 0105) Mapping Resource CS: 'DCMR'
(0040, db00) Template Identifier CS: '10011'
---------
(0040, a730) Content Sequence 20 item(s) ----
(0040, a010) Relationship Type CS: 'HAS CONCEPT MOD'
(0040, a040) Value Type CS: 'CODE'
(0040, a043) Concept Name Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '121058'
(0008, 0102) Coding Scheme Designator SH: 'DCM'
(0008, 0104) Code Meaning LO: 'Procedure reported'
---------
(0040, a168) Concept Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: 'P5-08000'
(0008, 0102) Coding Scheme Designator SH: 'SRT'
(0008, 0104) Code Meaning LO: 'Computed Tomography X-Ray'
---------
(0040, a730) Content Sequence 1 item(s) ----
(0040, a010) Relationship Type CS: 'HAS CONCEPT MOD'
(0040, a040) Value Type CS: 'CODE'
(0040, a043) Concept Name Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: 'G-C0E8'
(0008, 0102) Coding Scheme Designator SH: 'SRT'
(0008, 0104) Code Meaning LO: 'Has Intent'
---------
さらに、下を見てみると次のように表示されています。
(0008, 0104) Code Meaningには'CT Accumulated Dose Data'とあり、
その下に
(0008, 0104) Code Meaningには'Total Number of Irradiation Events'とあり、
またその下に
(0008, 0104) Code Meaningには'events'とあり、
次に
(0040, a30a) Numeric Valueには"7.0"とあります。
さらに、その下には
(0008, 0104) Code Meaning には'CT Dose Length Product Total'とあり、
その下に
(0008, 0104) Code Meaningには'mGy.cm'と単位があり、
次に
(0040, a30a) Numeric Valueには"377.94"と数値があります。
(0040, a010) Relationship Type CS: 'CONTAINS'
(0040, a040) Value Type CS: 'CONTAINER'
(0040, a043) Concept Name Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '113811'
(0008, 0102) Coding Scheme Designator SH: 'DCM'
(0008, 0104) Code Meaning LO: 'CT Accumulated Dose Data'
---------
(0040, a050) Continuity Of Content CS: 'SEPARATE'
(0040, a730) Content Sequence 2 item(s) ----
(0040, a010) Relationship Type CS: 'CONTAINS'
(0040, a040) Value Type CS: 'NUM'
(0040, a043) Concept Name Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '113812'
(0008, 0102) Coding Scheme Designator SH: 'DCM'
(0008, 0104) Code Meaning LO: 'Total Number of Irradiation Events'
---------
(0040, a300) Measured Value Sequence 1 item(s) ----
(0040, 08ea) Measurement Units Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '{events}'
(0008, 0102) Coding Scheme Designator SH: 'UCUM'
(0008, 0103) Coding Scheme Version SH: '1.4'
(0008, 0104) Code Meaning LO: 'events'
---------
(0040, a30a) Numeric Value DS: "7.0"
---------
---------
(0040, a010) Relationship Type CS: 'CONTAINS'
(0040, a040) Value Type CS: 'NUM'
(0040, a043) Concept Name Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '113813'
(0008, 0102) Coding Scheme Designator SH: 'DCM'
(0008, 0104) Code Meaning LO: 'CT Dose Length Product Total'
---------
(0040, a300) Measured Value Sequence 1 item(s) ----
(0040, 08ea) Measurement Units Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: 'mGy.cm'
(0008, 0102) Coding Scheme Designator SH: 'UCUM'
(0008, 0103) Coding Scheme Version SH: '1.4'
(0008, 0104) Code Meaning LO: 'mGy.cm'
---------
(0040, a30a) Numeric Value DS: "377.94"
---------
---------
---------
このように、Code Meaningの値やNumeric Valueの値を抽出することで必要な情報を得ることができるのではないかと考えました。
少し解析のフローをまとめると
- データエレメントのKeywordが“Content Sequence”が見つかったら解析を始める
- データエレメントのVRが“SQ“のとき、KeywordとValueを取得する
- データエレメントのkeywordが“Code Meaning”のとき、Valueの値を出力する
- 取得したValue（データエレメント）について２.の解析を繰り返す。
つまり、VRに"SQ"が見つかったらその下の階の項目を取り出し、その中のVRに"SQ"が見つかったら、さらにその下の。。。と繰り返すことで各階のCode Meaningの情報を取り出すこと考えました。
しかし、階層が3階層より下が上手く取り出すことができませんでした。
これはPydicomの仕様（バグ？）かわかりません。
そこで、タグ、キーワード、VRと値を取り出すメソッド（getInfoDCM）を作成して対応しました。
RDSRを読むPythonプログラム
RDSRを読み込むための準備を最初に書きます。
インポートするのPydicomだけです。
RDSRのファイル名を変数fnmに代入します。
rstxt, dstxtとsmryはSRの処理結果を格納する変数です。
nLvは"SQ"による階層の深さ、インデントの深さを代入する変数です。
prevLは前の"SQ"の行（データエレメント）の深さを格納する変数です。
inkeywordsには、抽出すべき値のキーワードを格納しています。
import pydicom
fnm = "./CT_rdsr.SR" # RDSR DICOM file name
rstxt = "" # For result text using getInfoDCM method
dstxt = "" # For DICOM dataset
smry = "" # For target text
nLv = 0 # Depth Level
prevL = 0 # For depth level of previous SQ line
# For extraction keyword
inkeywords = "CodeMeaning,TextValue,NumericValue,UID,DateTime"
データエレメントの深さ求めるメソッドとデータエレメントの各要素求めるメソッドについて
１．データエレメントの深さを求めるgetDepthメソッド
インデントは各行（データエレメント）の前の空白の数を数えて返します。
引数のdataはひとつのデータエレメントです。
def getDepth(data):
cnt = 0
for i in range(len(data)):
if data[i] == ' ':
cnt += 1
else:
break
return cnt
２．データエレメントの各要素求めるgetInfoDCMメソッド
このメソッドはデータエレメント中のタグ、キーワード、VRと値を抽出して返します。
引数のdataには (grp,elm) Key word VR: 'value' のようにひとつのデータエレメントが格納されています。
Pydicomでは"SQ"のシーケンスの終わりに'---------'が表示されます。
メソッドの最初にこの処理を記述します。
与えれたdataの内容が'---------'であった場合、
tagとVRには何も値を入れず、
キーワードには、”NoKeyword”と
値には、'---------'をそのまま入れて返します。
def getInfoDCM(data):
# ex) data is "(grp,elm) Key word VR: 'value' "
data = data.strip(" ") # Remove blanks on both sides of string
if data == '---------':
tag = ["",""]
kwd = "NoKeyword"
vr = ""
val = '---------'
return tag, kwd, vr, val
# Tag
tag = data.split(")")[0] # "(grp, elm"
tag = tag.replace("(", "").replace(" ","") # "grp, elm"
tag = tag.split(",") # [group, elemnt]
# vr
vr = ""
if "item(s)" in data: # ex) (0040, a043) Concept Name Code Sequence 1 item(s) ----
vr = "SQ"
else:
vr = data.split(":")[0][-2:] # ex) "(grp,elm) Key word VR"->"VR"
# value
val = ""
if ":" in data:
val = data.split(":")[-1] # [・・・ VR:value] -> value
elif "item(s)" in data:
val = ' '.join(data.split(" ")[-3:-1]) # -> n item(s)
val = val.replace("\'","").replace("\"","") # remove " and '
# keyword
kwd = ""
stat = data.find(")")+1 # Start position in keyword
end = data.find(vr) # End position in keyword
kwd = data[stat:end] # Extract keyword (remove (g,e))
if "item(s)" in kwd: # Case of VR="SQ"
kwd = kwd.split("item(s)")[0] # [kwd n,----] -> (g,e) kwd n
kwd = kwd.strip(" ").split(" ")[0:-1] # Extract the string before "n items"
kwd = ''.join(kwd).strip(" ") # Joint keyword and remove blanks
else:
kwd = kwd.replace(" ","") # Remove blank in keyword
return tag, kwd, vr, val
はじめにTagの処理を見ます。
(grp,elm) Key word VR: 'value'
グループ番号とエレメント番号は()の中にありますのでこれを頼りに２つの番号を抽出し、リストtagに代入します。
次に、VRです。
(grp,elm) Key word VR: 'value'
(grp,elm) Key word 99 item(s) ----
変数dataの中に"item(s)"があればVRは”SQ”です。
"SQ"以外の場合（else)、コロン（：）で変数dataの値を２つ分割し、前半[0]の最後の2文字[-2:]を抽出し変数vrに代入します。
次に、値です。
(grp,elm) Key word VR: 'value'
(grp,elm) Key word 99 item(s) ----
変数data中にコロン（：）があれば2つに分割し後半[1]を値として変数valに代入します。
コロンがなく、”item(s)”があれば、変数dataを空白で分割し”99 items(s)"のみ抽出しvalに代入します。
最後にvalの値の前後のシングル/ダブルクォーテーション（ ’ / " ）を削除します。
最後にキーワードです。
(grp,elm) Key word VR: 'value'
(grp,elm) Key word 99 item(s) ----
データエレメントの”）”の後ろの文字（stat番目）からVRの前（end)の文字列を抽出data[stat:end]して変数kwdに代入します。
しかし、この処理では、VRが"SQ"のときキーワードに"n item(s)"が存在します。
これを除く処理を行いました。
また、抽出したキーワード中の空白は除くようにしました。つまり”Study Date”は”StudyDate”になります。
これはPydicomのキーワードの表現に揃えるためです。
メインのプログラムを読む
はじめにpydicomを使ってRDSRのファイルを読み込みデータセットを変数dsに代入します。
SOPClassUIDなど必要な基本情報を表示します。
smryは必要な抽出項目の要約を代入する変数です。
変数dsは'pydicom.dataset.FileDataset'の型を取りますが、これを文字列型に変換して変数txdsに代入します。
txdsを改行（\n)で分割して行単位に分割してリストlndsに代入します。
読み込んだ情報などを表示していったん処理が停止（pause）します。
Enterキーで処理が進みます。
for文でlndsから1行（1データエレメント）ずつ取り出して処理します。
valueに'ContentSequence'が見つかったら、smryflagをTrueにして必要情報を抽出するように処理します。
変数rstxtには行番号、深さのレベル数とともに抽出した結果を追加代入(+=)しています。
一方、dstxtにはデータセットに行番号とデータエレメントを追加代入しています。
この２つの出力ファイルを比較することにより、間違いなく抽出できているか確認できます。
変数smryには注目するキーワードのデータエレメントの値を追加代入しています。
# ===== MAIN =====
ds = pydicom.dcmread(fnm) # Read DICOM SR
#print("ds=\n",ds)
print("SOP Class UID:", ds.SOPClassUID)
print("Modality :", ds.Modality)
print("Manufacturer :", ds.Manufacturer)
print("Study Date :", ds.StudyDate)
print("Pt Name :", ds.PatientName)
print("Pt ID :", ds.PatientID)
print("Modality :", ds.Modality)
print("Content Sequence:", len(ds.ContentSequence))
input("pause. Type ENTER key")
smry = "Study Date," + str(ds.StudyDate) + "\n"
smry += "Study Time," + str(ds.StudyTime) + "\n"
txds = str(ds) # Convert pydicom.dataset to String
print(txds)
print("\nData type :",type(ds),"\n\t\t-->",type(txds))
lnds = txds.split("\n") # Split text to lines
print("len(ds)=",len(lnds))
input("pause. Type ENTER key")
for i, ln in enumerate(lnds):
tag, keyword, vr, value = getInfoDCM(ln)
nLv = getDepth(ln)//3
if value in 'ContentSequence': # Start anlysis from this line
smryflag = True
sIndent = " "*(nLv*3) # Indent
rstxt += str(i)+",L:"+str(nLv)+sIndent+"("+str(tag)+"),"+keyword+","+vr+","+value+"\n"
dstxt += str(i)+","+ln + "\n"
print("value = >",value,"<")
# CR flag
if prevL > nLv or nLv <= 1:
crflag = True
else:
crflag = False
if smryflag and keyword in inkeywords:
#smry += "["+str(i)+"L"+str(+nLv)+"]"+str(value)+", "
smry += str(value)+", "
if vr == "SQ" and crflag: # VR is SQ, crflag is True
smry += "\n" # Then add CR in smry
if vr == "SQ": # VR is SQ,
prevL = nLv # Then assign now level to prevL
smry += rmCRline(smry)
with open("rdsrResult.txt","w") as f:
f.write(rstxt)
with open("rdsrDataset.txt","w") as f:
f.write(dstxt)
with open("rdsrSummary.csv","w") as f:
f.write(smry)
さて、ここで問題です。変数smryに値を追加していくだけですと、抽出したい情報がたった1行で出力されてしまいます。
そこで、意味ある情報の塊の境界で改行を入れる必要があります。
改行するかどうかを変数crflagを用いて制御しています。
if prevL > nLv or nLv <= 1:
crflag = True
else:
crflag = False
これは、1つ前の"SQ"のデータエレメントの階層レベルprevLが現在の階層レベルnLvより大きい（つまり現在の階層は前の階層より上がるとき）、または現在の階層が1以下のとき（つまり階層が0または1の深さにある）ときはcrflagをTrueにしています。
それ以外はcrflagをFalseに設定します。
そして、VRの値とcrflagの状態を見て変数smryに改行を入れています。
if vr == "SQ" and crflag: # VR is SQ, crflag is True
smry += "\n" # Then add CR in smry
このままでは改行が至る所に入ってしまうため、無駄な改行を削除するためrmCRlineメソッドを作成しました。
smry += rmCRline(smry)
実行してみる
プログラム実行してみましょう
初めに基本的な情報を表示して停止します。
SOP Class UID: 1.2.840.10008.5.1.4.1.1.88.67
Modality : SR
Manufacturer : SIEMENS
Study Date : 20220224
Pt Name : NO^Name
Pt ID : 12345678
Modality : SR
Content Sequence: 20
pause. Type ENTER key
ここで、「ENTER」キーを打つと、処理が次に進みます。
・・・・・・・
(0008, 0104) Code Meaning LO: 'Source of Dose Information'
---------
(0040, a168) Concept Code Sequence 1 item(s) ----
(0008, 0100) Code Value SH: '113856'
(0008, 0102) Coding Scheme Designator SH: 'DCM'
(0008, 0104) Code Meaning LO: 'Automated Data Collection'
---------
---------
Data type : <class 'pydicom.dataset.FileDataset'>
--> <class 'str'>
len(ds)= 3159
pause. Type ENTER key
読み込んだデータセットのデータの型が'pydicom.dataset.FileDataset'から'str'に変換されました。
「ENTER」キーを打つと、処理が次に進みます。
・・・・・・・
value = > CONTAINS <
value = > CODE <
value = > 1 item(s) <
value = > 113854 <
value = > DCM <
value = > Source of Dose Information <
value = > --------- <
value = > 1 item(s) <
value = > 113856 <
value = > DCM <
value = > Automated Data Collection <
value = > --------- <
value = > --------- <
最後に３つのファイルが出力されます。
rdsrDataset.txt データエレメントに行番後を付加して出力します。
rdsrResult.txt getInfoDCMメソッドで抽出した結果を行番号と一緒に出力します。
rdsrSummary.txt 必要な情報を行単位で出力します。
summayの結果を見てみる
rdsrSummary.csvファイルの出力結果は次のようになります。
最初に基本的な情報が表示されます。
以下ではCTのX線量に関する情報が単位とともに表示されています。
おわりに
今回、2022年4月24日に日本放射線技術学会 北海道支部 第78回春季学術大会 医療情報セミナーで発表したものを掲載しました。
DICOM SRの読み込みはPydicomのサンプルrtplan.srで試していましたが、今回RDSRは初めての経験でした。
読み込むまでは良かったのですが、階層が6階層まで下りて行く必要があり、Pydicomの標準機能ではうまく情報を取り込むことができませんでした。
CTのほかにAG（血管造影）のRDSRも解析しましたが、同じ結果でした。
本来はDICOM規格をしっかり理解したうえでプログラミングすべきでしたが、短期間に結果を出す必要があり、力ずくの対応でした。
したがってまだ気づかぬバグなどが潜んでいるかもしれません。プログラムの使用には注意をお願いします。
0 件のコメント:
コメントを投稿