mar_3

Autor:marianapbento
Data:10/03/2011

1. Criar uma lição sobre a função iamosaic, explicando passo a passo sua implementação. Use como exemplo as lições da toolbox IA636 (Lessons). Fique à vontade para alterar a implementação, caso ache necessário. Crie um exemplo numérico para testar a função.

Function Implementation - step by step

  1. FIRST STEP: ROIs choice

    • Compute Xini, that is equal to X-1(inicial image input parameter less 1). In case of X value be bigger than slice-1(maximum number os slices less 1), Xini receives slices-1 value.

    • Compute Xend, that is equal to X+N(inicial image input parameter plus number of images input parameter). In case of X value be bigger than slice(error on input parameters), Xend receives slice value.

    • In order to visualize planes XZ and YZ, some addition transformations are required:
      • First Transformation: rows inversion to correct acquisition procedure.
      • Second Transformation: planes transposition, that change the visualize plane.
  2. SECOND STEP: Slices exibition

    • Compute the necessary number of rows Nlinhas to exhibit the selected slices from step 1. Nlinhas = sliceS/n
    • Create a 1-column vector with all slices.
    • Transforme the 1-column vector on 1-row vector.
    • Split the 1-row vector on Nlinhas rows. Those rows only contain n slices, where n is equal to the number of slices per row input parameter.
    • Stack up all rows vertically, creating the result image.
  1 #import dicom
  2 import numpy as np
  3 from ia636 import ianormalize
  4 
  5 # reading dicom file
  6 filename = 'PHILIPS/DICOM/IM_0007'
  7 dataset = dicom.read_file(find_attachment_file(filename))
  8 data = dataset.pixel_array.astype(float64)
  9 
 10 def iamosaic(f, N, X, n, plane='XY'):
 11 
 12     import numpy as np # import numpy library
 13 
 14     slice,row,column = f.shape  # read input image shape
 15 
 16     # 1. FIRST STEP: ROIs choice
 17 
 18     # 1.1 Xini = X-1 or slice-1 (case X>slice)
 19     # 1.2 Xend = X+N ou slice (case input parameter are incorrect)
 20     # 1.3 froi = Interest slices, extracted by slicing the input image.
 21     # 1.4 Planes XZ e YZ transformation´s to visualization correction.
 22     # 1.4.1 First Transformation: rows inversion.
 23     # 1.4.2 Transformação 2: planes transpode: change the visualization plane.
 24 
 25     # XY plane
 26     if plane == 'XY':
 27         Xini = min(X,slice)-1
 28         Xend = min(X+N,slice+1)-1
 29         froi = f[Xini:Xend,:,:] # extraction of all rows and columns from interest slices.
 30         adshow(ianormalize(froi[(Xini+Xend)/2,:,:],"roi XY")) # selected slice sample exibition
 31     # XZ plane
 32     elif plane == 'XZ':
 33         Xini = min(X,row)-1
 34         Xend = min(X+N,row+1)-1
 35         froi = f[:,Xini:Xend,:] #  extraction of all slices and column from interest rows.
 36         adshow(ianormalize(froi[:,(Xini+Xend)/2,:],"roi XZ")) # selected slice sample exibition
 37         # Transformations to correctly display the XZ plane
 38         temp = froi[::-1,:,:]
 39         adshow(temp[:,(Xini+Xend)/2,:],"XZ inversion")
 40         froi = temp.transpose(1,2,0)
 41         adshow(ianormalize(froi[:,(Xini+Xend)/2,:],"XZ transpose"))
 42     # YZ plane
 43     elif plane == 'YZ':
 44         Xini = min(X,column)-1
 45         Xend = min(X+N,column+1)-1
 46         froi = f[:,:,Xini:Xend] # extraction of all slices and rows from interest columns.
 47         adshow(ianormalize(froi[:,:,(Xini+Xend)/2],"roi YZ")) # selected slice sample exibition
 48         # Transformations to correctly display the YZ plane
 49         temp = froi[::-1,:,:]
 50         adshow(ianormalize(temp[:,:,(Xini+Xend)/2],"YZ inversion"))
 51         froi = temp.transpose(2,1,0)
 52         adshow(ianormalize(froi[:,:,(Xini+Xend)/2],"YZ transpose"))
 53 
 54 
 55     # 2. SECOND STEP: Slices exibition
 56 
 57     # 2.1 Compute Nlinhas = sliceS/n
 58     # 2.2 Check if Nlinhas*n is equal to number of selected slices
 59     #     2.2.1 If different, complete the difference with zeros (black value).
 60     # 2.2 Create a 1-column vector with all slices.
 61     # 2.3 Transforme the 1-column vector on 1-row vector.
 62     # 2.4 Split the 1-row vector into Nlinhas(2.1) rows
 63     # 2.5 Stack up vertically all rows(calculate on 2.4)
 64 
 65     sliceS,rowS,columnS = froi.shape # read input image shape
 66 
 67     Nlinhas = np.ceil(float(sliceS)/n) # Compute the number of necessary lines to correct exibition
 68     Ifalt = Nlinhas*n-sliceS # Check if Nlinhas*n is equal to number of selected slices
 69     print "Ifalt =",Ifalt # show Ifalt value
 70     print "Num de Linhas = ",Nlinhas # show Nlinhas value
 71 
 72     temp0 = froi.reshape((sliceS*rowS,columnS)) # temp0 = slices 1-column vetor
 73     temp1 = np.hstack((np.vsplit(temp0,sliceS))) # transform temp0 in 1-row vector
 74     # hstack: stack input elements horizontally
 75     # vsplit: split temp0 into sliceS equal parts
 76     # temp1: create a 1-row vector with all slices
 77     adshow(temp1,"temp1 - 1-row vector")
 78 
 79     if Ifalt: # check if Ifalt ==1, it true complete the rest with zeros
 80         background = np.zeros((rowS,columnS*Ifalt))
 81         temp2 = np.hstack((temp1,background))
 82     else:
 83         temp2 = temp1
 84 
 85 
 86     # create rows vector with N elements, and stack them vertically
 87     g = np.vstack((np.hsplit(temp2,Nlinhas)))
 88     adshow(g,"Final Result")
 89     g = g.astype(np.int32) # conversion to int32 type
 90 
 91     return g
 92 
 93 # check data dimension
 94 linhas = dataset.Rows
 95 colunas = dataset.Columns
 96 fatias = dataset.NumberofFrames
 97 
 98 # iamosaic example
 99 #g = iamosaic(data,fatias,1,15,'XY')
100 g = iamosaic(data,fatias,1,15,'XZ')
101 #g = iamosaic(data,fatias,1,15,'YZ')
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 100, in <module>
  File "<string>", line 36, in iamosaic
  File "/awmedia/www/packages/ia636/ianormalize.py", line 15, in ianormalize
    lower = range[0]
IndexError: 0-d arrays can't be indexed

------------------------------------------------------------

Possible Modifications

  • Split iamosaic function into three different functions. First one, responsible to select interesting slices, second function transpose axis, and the last would be responsible for display slices.
  • Take off row inversion procedure(froi(::-1,:,:))
  • Correct transpose procedure and it´s correspond plane
  • Split functions: http://parati.dca.fee.unicamp.br/adesso/wiki/courseIA369O1S2011/iamosaic/view/

Numeric Validation

 1 #import dicom
 2 import numpy as np
 3 #from ia636 import *
 4 
 5 a = arange(30).reshape(3,2,5)
 6 print "a = ",a
 7 print
 8 
 9 # check data dimension
10 fatias,linhas,colunas = a.shape
11 
12 # PLANO XY
13 
14 b = iamosaic(a,fatias,1,3,'XY') # 3 slices per row
15 print "b = ",b
16 print
17 
18 c = iamosaic(a,fatias,1,1,'XY') # 1 slice per row
19 print "c =",c
20 print
21 
22 d = iamosaic(a,2,1,3,'XY') # shows only two slices: first and second one
23 print "d =",d
24 print
25 
26 e = iamosaic(a,2,2,2,'XY') # shows only two slices: the second and third one
27 print "e =",e
28 print
29 
30 # PLANO YZ
31 
32 f = iamosaic(a,fatias,1,3,'YZ') # 3 slices per row - plane XY
33 print "f = ",f
34 print
35 
36 # PLANO XZ
37 
38 g = iamosaic(a,fatias,1,3,'XZ') # 3 slices per row - plane XZ
39 print "g = ",g
40 print
a =  [[[ 0  1  2  3  4]
  [ 5  6  7  8  9]]

 [[10 11 12 13 14]
  [15 16 17 18 19]]

 [[20 21 22 23 24]
  [25 26 27 28 29]]]

b =  [[ 0  1  2  3  4 10 11 12 13 14 20 21 22 23 24]
 [ 5  6  7  8  9 15 16 17 18 19 25 26 27 28 29]]

c = [[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]
 [20 21 22 23 24]
 [25 26 27 28 29]]

d = [[ 0  1  2  3  4 10 11 12 13 14  0  0  0  0  0]
 [ 5  6  7  8  9 15 16 17 18 19  0  0  0  0  0]]

e = [[10 11 12 13 14 20 21 22 23 24]
 [15 16 17 18 19 25 26 27 28 29]]

f =  [[20 10  0 21 11  1 22 12  2]
 [25 15  5 26 16  6 27 17  7]]

g =  [[20 10  0 25 15  5  0  0  0]
 [21 11  1 26 16  6  0  0  0]
 [22 12  2 27 17  7  0  0  0]
 [23 13  3 28 18  8  0  0  0]
 [24 14  4 29 19  9  0  0  0]]

2. Estude as funções da toolbox IA6363: iaramp, iapplylut, ianormalize. Isso irá ajudá-los a entender a função adshow. Crie uma documentação para a função adshow, explorando o que acontece quando chamamos esta função com imagens de tipos diferentes (uint8, uint16, int8, int16, booleano, ponto flutuante, etc). Tente sugerir ao usuário qual a forma de usar a função adshow com os diversos tipos de variáveis da maneira mais adequada.

Functions iaramp, iapplylut, ianormalize

  • iaramp: create an image with vertical bands with crescent values of input level.
1 from ia636 import iaramp
2 
3 F = iaramp([5,3], 3,[0,255]) # matriz size: [5,3]; number of levels: 3; level range: [0,255]
4 G = iaramp([3,2], 2, [0,1])  # matriz size: [3,2]; number of levels: 2; level range: [0,1]
5 print "Gray Level",F
6 print "Boolean",G
Gray Level [[  0 127 255]
 [  0 127 255]
 [  0 127 255]
 [  0 127 255]
 [  0 127 255]]
Boolean [[0 1]
 [0 1]
 [0 1]]
  • iaapplylut: Apllies the transformation into input image
 1 from ia636 import iaapplylut
 2 from numpy import *
 3 
 4 f = arange(10).reshape(2,5)
 5 print 'f=\n',f
 6 it = array(range(10)) # transformation = identity
 7 print 'it =',it  # visualize transformation
 8 print
 9 g = iaapplylut(f, it) # transformation application
10 print 'g=\n',g # visualize transformation result
11 print
12 
13 itn = 5 - it  # negative of a image
14 print 'itn =',itn  # visualize transformation
15 print
16 gn = iaapplylut(f, itn) # transformation application
17 print 'gn=\n',gn  # visualize transformation result
18 print
19 
20 itv = it
21 itv[0] = 123 # change of specific value
22 print 'itv =',itv  # visualize transformation
23 gv = iaapplylut(f, itv) # transformation application
24 print 'gv=\n',gv # visualize transformation result
f=
[[0 1 2 3 4]
 [5 6 7 8 9]]
it = [0 1 2 3 4 5 6 7 8 9]

g=
[[0 1 2 3 4]
 [5 6 7 8 9]]

itn = [ 5  4  3  2  1  0 -1 -2 -3 -4]

gn=
[[ 5  4  3  2  1]
 [ 0 -1 -2 -3 -4]]

itv = [123   1   2   3   4   5   6   7   8   9]
gv=
[[123   1   2   3   4]
 [  5   6   7   8   9]]
  • ianormalize: normilize pixels from input image into a fixed range of values
 1 from ia636 import ianormalize
 2 
 3 f = array([1., 50., 100.]) # positive values normalization example
 4 g1 = ianormalize(f, [0,255]) # example1: normalization into gray level value
 5 print 'g1 =',g1
 6 g2 = ianormalize(f, [-1,1]) # example2: normalization into a small range of negative and positive values
 7 print 'g2 =',g2
 8 g3 = ianormalize(f, [0,1]) # example 3: normalization into a [0,1] range
 9 print 'g3 =',g3
10 
11 f = array([-100., 0., 100.]) # positive and negative values normalization example
12 g4 = ianormalize(f, [0,255])
13 print 'g4 =',g4
14 g5 = ianormalize(f, [-1,1])
15 print 'g5 =',g5
g1 = [   0.          126.21212121  255.        ]
g2 = [-1.         -0.01010101  1.        ]
g3 = [ 0.          0.49494949  1.        ]
g4 = [   0.   127.5  255. ]
g5 = [-1.  0.  1.]

Função adshow

Sinopse

Show gray level images.

  • adshow(f,title = '')
    • f: Image. input image.
    • title: image title.

Description

The function adshow is a image visualization tool. The basic exibition scale is gray level [0,255].

Examples

Auxiliary function: create edges on input images

1 def edge(x):
2     a,b = x.shape
3     x[0::a-1,:] = 0
4     x[:,0::b-1] = 0
5     return x

Example 1: using function ramp

 1 from ia636 import iaramp
 2 
 3 F = iaramp([100,100], 5, [0,255])
 4 G = iaramp([100,100], 2, [0,1])
 5 H = iaramp([100,100], 5, [0.,255.])
 6 I = iaramp([100,100], 5, [-100,125])
 7 Y = iaramp([100,100], 5, [-100,100])
 8 adshow(edge(F),"UINT8 ramp")
 9 adshow(edge(G),"BOOLEAN ramp")
10 adshow(edge(H),"FLOAT ramp")
11 adshow(edge(I),"SIGNED INT ramp")
12 adshow(edge(Y),"SIGNED INT ramp 2")
ERROR execute

------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 11, in <module>
  File "xs_runner_python.py", line 194, in adshow
    raise IOError('Image is out of range: min=%d max=%d, use ianormalize' % (m,M))
IOError: Image is out of range: min=-100 max=125, use ianormalize

------------------------------------------------------------

UINT8 ramp

BOOLEAN ramp

FLOAT ramp

Example 2: using function normalize

1 from ia636 import iaramp
2 from ia636 import ianormalize
3 
4 # utilização da função normalize
5 adshow(edge(ianormalize(G, [0,255])),"BOOLEAN normalize")
6 adshow(edge(ianormalize(I, [0,255])),"SIGNED INT normalize")
7 adshow(edge(ianormalize(H, [0,255])),"SIGNED INT 2 normalize")

BOOLEAN normalize

SIGNED INT normalize

SIGNED INT 2 normalize

Example 3: using function iaapplylut

 1 from ia636 import iaramp
 2 from ia636 import iaapplylut
 3 from ia636 import ianormalize
 4 
 5 F = iaramp([100,100], 10, [0,255])
 6 
 7 it = arange(256)
 8 G = iaapplylut(F, it)
 9 adshow(edge(G), title='Identity Transformation')
10 
11 itn = 255 - it
12 H = iaapplylut(F, itn)
13 adshow(edge(H), title='Inverse Transformation')
14 
15 itv = it
16 itv[255] = 128
17 I = iaapplylut(F, itv)
18 adshow(edge(I), title='Specific color transformation')
19 
20 itt = uint8(zeros(256))
21 itt[0:100] = 1
22 itt[101:180] = 2
23 itt[181:256] = 3 # Generalized Threshold
24 
25 V = iaapplylut(F, itt) # result has 3 possible levels [1,2,3]
26 adshow(edge(ianormalize(V)), title='Generalized Threshold')# receives a normalization process to visualization of 3 different ranges.

Identity Transformation

Inverse Transformation

Specific color transformation

Generalized Threshold

Example 4: using numerical samples

 1 from ia636 import ianormalize
 2 from ia636 import iaapplylut
 3 from numpy import *
 4 from ia636 import iaramp
 5 
 6 a = arange(10000).reshape(100,100)
 7 
 8 thr = arange(10000) # transformada identidade
 9 k = uint8(iaapplylut(a, thr)) # tipo UINT8
10 print "tipo =",type(k[(0,0)])
11 
12 x = int8(k*-1) # tipo INT8
13 print "tipo =",type(x[(0,0)])
14 
15 g = uint16(a) # tipo UINT16
16 print "tipo =",type(g[(0,0)])
17 
18 y = int16(g*-1) # tipo SIGNED UINT16
19 print "tipo =",type(y[(0,0)])
20 
21 h = zeros( (100,100) )
22 for i in range (0, 100):
23     for j in range (0, 100):
24         h[(i,j)] = float(a[(i,j)])  # type= FLOAT
25 
26 print "tipo =",type(h[(0,0)])
27 
28 f = a>4999
29 print "tipo =",type(f[(0,0)]) # type= BOOLEAN
30 
31 # adshow application
32 
33 adshow(edge(k),"matrix uint8")
34 adshow(edge(x),"matrix signed uint8")
35 adshow(edge(g),"matrix uin16")
36 adshow(edge(y),"matrix signed uint16")
37 adshow(edge(h),"matrix float")
38 adshow(edge(f),"matrix boolean")
39 
40 #m = iaramp([100,100],5,[100,255])
41 #adshow(edge(m),"teste")
ERROR execute
tipo = <type 'numpy.uint8'>
tipo = <type 'numpy.int8'>
tipo = <type 'numpy.uint16'>
tipo = <type 'numpy.int16'>
tipo = <type 'numpy.float64'>
tipo = <type 'numpy.bool_'>
------------------------------------------------------------
*** Exception while evaluating code:
  File "<string>", line 34, in <module>
  File "xs_runner_python.py", line 194, in adshow
    raise IOError('Image is out of range: min=%d max=%d, use ianormalize' % (m,M))
IOError: Image is out of range: min=-128 max=127, use ianormalize

------------------------------------------------------------

matrix uint8

Final Considerations

  1. Iamosaic function can be diveded into two main blocks: the first one is to choose the imagens that will be visualized. The second main point is actually show the imagens side by side. The procedure validation is provided by the numerical example.
  2. Adshow function is used in many applications to visualize imagens into a gray level range. In case of small or limited range input images, it is required to use a previous normalization process, other way, it will be impossible to recognize different regions. Another verification is that the same data stored into different type variables, generates similar images.