RAW形式ボリュームデータをVTI形式に変換する

2020年5月18日資料：医用画像,資料：生体シミュレーションソフトウェア開発

signed 16bit で格納されているRAW形式のボリュームデータを，VTK ImageDataファイル (vti) 形式に変換する Python プログラム．

VTI形式には Legacy と XML とがありますが，これはXML．
XMLにもテキストとバイナリとがありますが，これはテキスト．

色々と試行錯誤のあとがありますが，ご容赦ください．

使い方は

$ python3 raw2vtitxt.py in.raw out.vti width height slice

です．

コード

import sys
import numpy as np
import xml.etree.ElementTree as ET

def generateVtiTree():
    # XML ドキュメントの作成
    vtket = ET.Element('VTKFile',{'type':'ImageData'})
#    imagedata = ET.SubElement(vtkfile, 'ImageData', {'WholeExtent':'0 2 0 2 0 2'})
    imagedata = ET.SubElement(vtket,'ImageData')
    imagedata.set('WholeExtent','0 2 0 2 0 2')
    imagedata.set('Origin','0 0 0')
    imagedata.set('Spacing','1.0 1.0 1.0')
#    imagedata.text = '<Hello> world!'
    piece = ET.SubElement(imagedata, 'Piece', {'Extent':'0 2 0 2 0 2'})

    celldata = ET.SubElement(piece,'CellData',{'Scalars':'cell_values'})

    dataarray = ET.SubElement(celldata,'DataArray',{'Name':'cell_values'})
    dataarray.set('type','Int16')
    dataarray.set('format','ascii')

    dataarray.text = '1 2 3 4 5 6 7 8'

    return vtket



def setDataarray(vtket,vtext):
    for dataarray in vtket.iter('DataArray'):
        print;

#    dataarray.text = '8 7 6 5 4 3 2 1'
    dataarray.text = vtext

    return

def setSize(vtket,width,height,slice):

    strSize = '0 {0} 0 {1} 0 {2}'.format(width,height,slice)
#    print(strSize)

    for imagedata in vtket.iter('ImageData'):
        print;

    imagedata.set('WholeExtent',strSize)

    for piece in vtket.iter('Piece'):
        print;

    piece.set('Extent',strSize)


    return


def loadRawVolume( fp ):
    inVolume = np.fromfile(fp,np.int16,-1)
    fp.close()

    return inVolume


def convertNumpyToText(inVolume):

    mapped_list = map(str,inVolume)
    vtext = ' '.join(mapped_list)

    return vtext


def run(pathInfile,pathOutfile,width,height,slice):

    # load RAW to numpy
    fpInfile = open(pathInfile, 'rb')
    inVolume = loadRawVolume(fpInfile)
    print(inVolume.size)

    # convert numpy to text
    vtext = convertNumpyToText(inVolume)
    print(len(vtext))

    vtket = generateVtiTree()
#    ET.dump(vtket)
    setDataarray(vtket,vtext)
    setSize(vtket,width,height,slice)
#    ET.dump(vtket)
    etree = ET.ElementTree(vtket)
    etree.write(pathOutfile)


if __name__ == '__main__':
    args = sys.argv

    if len(args) == 6:
        pathInfile = args[1]
        pathOutfile = args[2]
        argsWidth = args[3]
        argsHeight = args[4]
        argsSlice = args[5]

        run(pathInfile,pathOutfile,int(argsWidth),int(argsHeight),int(argsSlice))

    else:
        print('Usage: python3 raw2vtitxt.py <infile> <outfile> <width> <height> <slice>')
        quit()

参考

XMLベースのvtk形式に関する覚書
https://qiita.com/hsimyu/items/a8fba8866c5f5c2d8b78