Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
1.0k views
in Technique[技术] by (71.8m points)

r - What could be the reason for bad dicom image plot

What could be the reason that the dicom file of this usual x-ray is getting plotted in a messed up manner:

enter image description here

The algorithm used is as follows:

The original image matrix is 3d:

int [1:2014, 1:2014, 1:3] 110 51 99 113 52 101 111 53 102 110 ...

This rgb is converted to gray scale by formula:

gray = 0.3*mat[,,1] + 0.59*mat[,,2] + 0.11*mat[,,3] ; 

And then it is plotted after specifying colors as:

grey(0:64/64)

Where could be the error?

I am using oro.dicom package in R with function:

jj = readDICOMFile(fname, endian = "little", flipud = TRUE, DICM = TRUE, skipSequence = FALSE, pixelData = TRUE, warn = -1, debug = FALSE) 

and it returns a the matrix jj$img whose structure is:

int [1:2014, 1:2014, 1:3] 110 51.... 

I then convert it to gray and plot it. If it was rgba, the matrix would have been 2014*2014*4 rather than *3. The header of dicom image mentions "PhotometricInterpretation" as "RGB". The header also mentions rows and columns as 2014 each. Could it be related to bit problem: leadtools.com/sdk/medical/dicom-spec17.htm

Edit: Bits allocated is 8, bits stored is 8 and highBit is 7.

Following is the link of sample dicom image which has similar image matrix and give similar error: http://www.barre.nom.fr/medical/samples/files/US-RGB-8-esopecho.gz

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

The readDICOMFile might have a bug. You can fix by rearrange the image array:

jj = readDICOMFile(fname, flipud = FALSE, DICM = TRUE, skipSequence = FALSE, pixelData = TRUE, warn = -1, debug = FALSE)
img <- jj$img # extract image part
img <- aperm(array(c(aperm(img, c(2, 1, 3))), c(3, 256, 120)), c(3, 2, 1)) # rearrange dimension
img <- img[120:1,,] # flip ud
grid::grid.raster(scales::rescale(img))

enter image description here


UPDATE

readDICOMFile has another bug. This is what you want. You may be better to report this but to the authors of oro.dicom.

img <- jj$img # extract image part
img <- aperm(array(c(aperm(img, c(2, 1, 3))), c(3, 256, 120)), c(3, 2, 1)) # rearrange dimension

# conversion b/w unsigned and signed
img <- ifelse(img > 0, img, 256+img)

# window-ing
wc <- 127
ww <- 255

ymin <- 0
ymax <- 1

img2 <- ifelse(img <= wc - 0.5 - (ww-1)/2, ymin, 
               ifelse(img > wc - 0.5 + (ww-1)/2, ymax,
                      ((img - (wc - 0.5)) / (ww - 1) + 0.5) * (ymax - ymin) + ymin
                      ))

grid::grid.raster(img2)

enter image description here


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...