Exercício #3

Autor: Guilherme O. Franchi Junior
Data: 11/03/2011

iamosaic Review

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.

Synopsis

It creates a mosaic of images from a volumetric image (3D).

  • g = iamosaic(img, N_img, I_img, r_img, plane='XY')
    • g: Output image
    • img : Image. input image;
    • N_img: number of images;
    • I_img: Initial/Starting image;
    • r_img: Images per row;
    • plane: String/chosen plane. Possible values: 'XY', 'XZ', 'YZ'.
 1 def iamosaic(f, N, X, n, plane='XY'):
 2 
 3     import numpy as np
 4 
 5     slice,row,column = f.shape
 6 
 7     # XY plane
 8     if plane == 'XY':
 9         Xini = min(X,slice)-1
10         Xend = min(X+N,slice+1)-1
11         froi = f[Xini:Xend,:,:]
12     # XZ plane
13     elif plane == 'XZ':
14         Xini = min(X,row)-1
15         Xend = min(X+N,row+1)-1
16         froi = f[:,Xini:Xend,:]
17         # Transformations to correctly display the XZ plane
18         #temp = froi[::-1,:,:]
19         temp = froi[:,:,:]
20         froi = temp.transpose(1,2,0)
21     # YZ plane
22     elif plane == 'YZ':
23         Xini = min(X,column)-1
24         Xend = min(X+N,column+1)-1
25         froi = f[:,:,Xini:Xend]
26         # Transformations to correctly display the YZ plane
27         #temp = froi[::-1,:,:]
28         temp = froi[:,:,:]
29         froi = temp.transpose(2,1,0)
30 
31     sliceS,rowS,columnS = froi.shape
32 
33     Nlinhas = np.ceil(float(sliceS)/n)
34     Ifalt = Nlinhas*n-sliceS
35 
36     temp0 = froi.reshape((sliceS*rowS,columnS))
37     temp1 = np.hstack((np.vsplit(temp0,sliceS)))
38 
39     if Ifalt:
40         background = np.zeros((rowS,columnS*Ifalt))
41         temp2 = np.hstack((temp1,background))
42     else:
43         temp2 = temp1
44 
45     g = np.vstack((np.hsplit(temp2,Nlinhas)))
46     g = g.astype(np.int32)
47 
48     return g

Description

This function returns an image containing juxtaposed slices of the input volume chosen by the user.

This function receives the following parameters:

  • An input image (3D Image);
  • Total of images which will be displayed on the mosaic image - Integer;
  • An initial image, from the list of images, which will be set to be first displayed - Index/Integer;
  • Total of images which will be displayed by each row - Integer;
  • From which plane the new mosaic image will be created - String;

The firts process of this method is to get how many slices, rows and collumns the input image does have. With that in place, the next step gets which plane will be visualized, and gets only the "region of interest" of the input image based on starting image plus total of images which will be displayed. An inversion of images position from axis X and Y is also made, since their positions are inverted. Then, a transpose of desired axis is done to bring them to the front.

New variables are created containing the new shape of that new 3D image with desired region of interest (froi variable). Some math calculations are made to determine how many lines the mosaic image should have, based on number of images to be displayed on each line plus total of images being displayed. An Ifalt variable is created to determine if that new image will have "blank" fields in its last line.

The temp0 variable will reshape the froi image to bring all images from each slice, concatenated vertically (2D image). Another variable (temp1) will contain a new image, based on temp0. First of all, each image on that "reshaped" image will be separated on individual arrays (using vsplit) and then, each of those arrays are going to be concatenated in a single line of images using hstack.

The next step is to verify if there is any blank field which will need to be filled with "black" values to complete our image, based on Ifalt variable. Otherwise, temp2 will be referencing temp1 image.

Then, the final image is generated based on temp2 image (with ROI). That single line of images is splitted (cut) based on the total lines which will be displayed on the mosaic image. Each single line is divided using hsplit which will return one array for each line. Then, those arrrays are concatenated vertically, one bellow the other, generating our final image, using vstack.

As a final step, that final image is converted as int32 type and returned from the iamosaic function.

Numerical Example

This single example shows iamosaic usage, based in a numerical matrix with sequential numbers.

 1 # Creates an array of numbers, starting on 0 till 63
 2 tmpArray = arange(64, dtype='uint64')
 3 print "Array =\n", tmpArray
 4 # Transforms it in a three dimensional array
 5 imgMatrix = tmpArray.reshape(4,4,4)
 6 print "3D Array=\n", imgMatrix
 7 # Using iamosaic to display that 3D array in mosaic mode
 8 finalArray = iamosaic(imgMatrix, size(imgMatrix), 1, 2, 'XY')
 9 print "Final Array (XY Plane)=\n", finalArray
10 finalArray = iamosaic(imgMatrix, size(imgMatrix), 1, 2, 'YZ')
11 print "Final Array (YZ Plane)=\n", finalArray
12 finalArray = iamosaic(imgMatrix, size(imgMatrix), 1, 2, 'XZ')
13 print "Final Array (XZ Plane)=\n", finalArray
Array =
[ 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 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
 50 51 52 53 54 55 56 57 58 59 60 61 62 63]
3D Array=
[[[ 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 30 31]]

 [[32 33 34 35]
  [36 37 38 39]
  [40 41 42 43]
  [44 45 46 47]]

 [[48 49 50 51]
  [52 53 54 55]
  [56 57 58 59]
  [60 61 62 63]]]
Final Array (XY Plane)=
[[ 0  1  2  3 16 17 18 19]
 [ 4  5  6  7 20 21 22 23]
 [ 8  9 10 11 24 25 26 27]
 [12 13 14 15 28 29 30 31]
 [32 33 34 35 48 49 50 51]
 [36 37 38 39 52 53 54 55]
 [40 41 42 43 56 57 58 59]
 [44 45 46 47 60 61 62 63]]
Final Array (YZ Plane)=
[[ 0 16 32 48  1 17 33 49]
 [ 4 20 36 52  5 21 37 53]
 [ 8 24 40 56  9 25 41 57]
 [12 28 44 60 13 29 45 61]
 [ 2 18 34 50  3 19 35 51]
 [ 6 22 38 54  7 23 39 55]
 [10 26 42 58 11 27 43 59]
 [14 30 46 62 15 31 47 63]]
Final Array (XZ Plane)=
[[ 0 16 32 48  4 20 36 52]
 [ 1 17 33 49  5 21 37 53]
 [ 2 18 34 50  6 22 38 54]
 [ 3 19 35 51  7 23 39 55]
 [ 8 24 40 56 12 28 44 60]
 [ 9 25 41 57 13 29 45 61]
 [10 26 42 58 14 30 46 62]
 [11 27 43 59 15 31 47 63]]

Image Example

It presents a mosaic for a given volume, the first eigenvalue of a Diffusion Tensor Image of the brain.

 1 # Read DICOM file
 2 filename = 'PHILIPS/DICOM/IM_0007'
 3 dataset = dicom.read_file(find_attachment_file(filename))
 4 
 5 # verify dimension of input image (3D)
 6 imglines = dataset.Rows
 7 imgcolumns = dataset.Columns
 8 imgslices = dataset.NumberofFrames
 9 
10 print "Total Number of Lines = ", imglines
11 print "Total Number of Columns = ", imgcolumns
12 print "Total Number of Slices = ", imgslices
13 
14 data = dataset.pixel_array.astype(float64)
15 
16 mAxi = iamosaic(data, imgslices, 1, 20, 'XY')
17 mCor = iamosaic(data, imgslices, 1, 20, 'XZ')
18 adshow(mAxi, title = "Axial Mosaic")
19 adshow(mCor, title = "Coronal Mosaic")
Total Number of Lines =  224
Total Number of Columns =  224
Total Number of Slices =  180

Axial Mosaic

Coronal Mosaic

Proposed Modifications

  1. Remove the line where the inversion happens (temp = froi[::-1,:,:]). This is not necessary as the image is already positioned and ordered correctly on each plane.
  2. Create a new way of copying that image inside of a black image with width/height resolution based on number of lines + total of images being displayed (capturing the size of each individual image to generate that new image) and, then, using slices.
  3. Change plane naming convention. Instead of using XY, should be easier to use Z, instead of XZ, Y, and instead of YZ, X.

ia636 Functions

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.

Function iaramp

This function creates a ramp of values, based on input provided. It accepts three parameters:

  • A height/width for the result image;
  • Number of bands (vertical) that the image will have;
  • Minimum and maximum values that will be using between bands;

Numeric Example

1 from ia636 import *
2 from numpy import *
3 
4 ramp = iaramp([10,10], 4, [2,8])
5 print "Ramp Array=\n", ramp
Ramp Array=
[[2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]
 [2 2 2 4 4 6 6 6 8 8]]

Image Example

1 ramp = iaramp([200,180], 10, [30,180])
2 adshow(ramp, title='Ramped Image')

Ramped Image

Function iaapplylut

This function will apply an intensity image transform to the input image, used as parameter. This function accepts both colored images or images in gray scale. The intensity transform is represented by a table where the input (gray scale) color address the table line and its column contents indicates the output (gray scale) image color. In resume, the image received as parameter will be used as indices of the new image, and the table will contains the values which will be used by those indices from image.

Numeric Example

This example will demonstrate how iaapplylut is used and what happens when a new table is used to multiply each value of the original bi-dimensional array by 2. What happens is that each value in the table created is mapped to the value from the original array.

 1 arr = arange(256)
 2 arr[80] = 255
 3 arr[200] = 30
 4 mat = arr.reshape(16,16)
 5 print "Original 2D Array=\n", mat
 6 
 7 
 8 # create a new array that will be used as index
 9 # This new table will multiply each value of original array by two
10 tmp = arange(256) * 2
11 # Make some changes on this table to see the difference
12 tmp[30] = 0
13 tmp[255] = 1
14 
15 final_arr = iaapplylut(mat, tmp)
16 print "Final 2D Array (Multiplied by 2)=\n", final_arr
Original 2D Array=
[[  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  30  31]
 [ 32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47]
 [ 48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63]
 [ 64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79]
 [255  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95]
 [ 96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111]
 [112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127]
 [128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143]
 [144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159]
 [160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175]
 [176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191]
 [192 193 194 195 196 197 198 199  30 201 202 203 204 205 206 207]
 [208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223]
 [224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239]
 [240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255]]
Final 2D Array (Multiplied by 2)=
[[  0   2   4   6   8  10  12  14  16  18  20  22  24  26  28  30]
 [ 32  34  36  38  40  42  44  46  48  50  52  54  56  58   0  62]
 [ 64  66  68  70  72  74  76  78  80  82  84  86  88  90  92  94]
 [ 96  98 100 102 104 106 108 110 112 114 116 118 120 122 124 126]
 [128 130 132 134 136 138 140 142 144 146 148 150 152 154 156 158]
 [  1 162 164 166 168 170 172 174 176 178 180 182 184 186 188 190]
 [192 194 196 198 200 202 204 206 208 210 212 214 216 218 220 222]
 [224 226 228 230 232 234 236 238 240 242 244 246 248 250 252 254]
 [256 258 260 262 264 266 268 270 272 274 276 278 280 282 284 286]
 [288 290 292 294 296 298 300 302 304 306 308 310 312 314 316 318]
 [320 322 324 326 328 330 332 334 336 338 340 342 344 346 348 350]
 [352 354 356 358 360 362 364 366 368 370 372 374 376 378 380 382]
 [384 386 388 390 392 394 396 398   0 402 404 406 408 410 412 414]
 [416 418 420 422 424 426 428 430 432 434 436 438 440 442 444 446]
 [448 450 452 454 456 458 460 462 464 466 468 470 472 474 476 478]
 [480 482 484 486 488 490 492 494 496 498 500 502 504 506 508   1]]

Image Example

This example shows what happens when a table with negated values is used to generate a new image. This demonstrate how to use it on image files.

1 imagem = adreadgray('cowboy_color.jpg')
2 adshow(imagem, title="Original Image")
3 
4 # Creates a new array of values but now with inverted values
5 n_vector = uint8(255 - arange(256))
6 
7 # Uses the created vector as parameter for index generation, and thos new values
8 n_imagem = iaapplylut(imagem, n_vector)
9 adshow(n_imagem, title="Transformed Image")

Original Image

Transformed Image

Function ianormalize

This function is used to normalize an imagem between a frame of values. It will readjust your image with the values between a minimun and maximum values passed as parameters.

Numeric Example

 1 arr01 = array([4., 8., 78., 10., 15])
 2 print "Original Array 1=\n", arr01
 3 n_arr01 = ianormalize(arr01, [30, 40])
 4 print "Normalized Array 1=\n", n_arr01
 5 
 6 arr02 = array([10, 30, 60, 90])
 7 print "Original Array 2=\n", arr02
 8 n_arr02 = ianormalize(arr02, [-1, 1])
 9 print "Normalized Array 2=\n", n_arr02
10 
11 arr03 = array([1, 0, 0, 1])
12 print "Original Array 3=\n", arr03
13 n_arr03 = ianormalize(arr03, [0, 255])
14 print "Normalized Array 3=\n", n_arr03
Original Array 1=
[  4.   8.  78.  10.  15.]
Normalized Array 1=
[30 30 40 30 31]
Original Array 2=
[10 30 60 90]
Normalized Array 2=
[-1  0  0  1]
Original Array 3=
[1 0 0 1]
Normalized Array 3=
[255   0   0 255]

Image Example

1 image = adreadgray('cowboy_color.jpg')
2 adshow(image, title="Original Image")
3 
4 # Normalize our original image with values between 50 ~ 180
5 n_image = ianormalize(image, [50,180])
6 adshow(n_image, title="Normalized Image")

Original Image

Normalized Image

Function adshow

The main purpose of this function is to display an image to the user's page. This is a function from adpil module, from Adesso Wiki. The adpil module provides a link between numpy arrays and Python Imaging Library (PIL) images.

Synopsis

Display an image from a numpy array in a Tk toplevel with title given by argument 'title' and window id given by argument 'id'.

  • g = adshow(arr, title='adshow', id=0)
    • g: Output image
    • arr : numpy array to be displayed.
    • title : Title of the view window, which will be added as legend at the bottom of the image. Default = 'adshow'
    • id : Id of the Window being displayed. Default = 0

Examples

 1 # Displaying a colored image
 2 image = adread('cowboy_color.jpg')
 3 adshow(image, title='Colored Image')
 4 
 5 # Displaying a gray image
 6 image = adreadgray('cowboy_color.jpg')
 7 adshow(image, title='Gray Image')
 8 
 9 # Displaying array images created from iaaplylut function
10 adshow(mat, title='Normalized Array 02')
11 adshow(final_arr, title='Normalized Array 03')
12 
13 # Using a new image with negative numbers
14 tmp_img = image * -1.0
15 adshow(tmp_img, title='Gray Image with negative numbers')
16 
17 # Using ianormalize to fix -255 numbers
18 tmp_img = ianormalize(tmp_img, [0,255])
19 adshow(tmp_img, title='Gray Image with negative numbers')
20 
21 # Use a generated image (from ianormalize bellow) and check max and min values
22 norm_img = adreadgray('cowboy_mod.png')
23 print "Min Value from Image: ", min(ravel(norm_img))
24 print "Max Value from Image: ", max(ravel(norm_img))
Min Value from Image:  50
Max Value from Image:  180

Colored Image

Gray Image

Normalized Array 02

Normalized Array 03

Gray Image with negative numbers

Gray Image with negative numbers

Final Conclusions

  1. Copying iamosaic function to this page and having it detailed, step by step, was easier to understand. Some numeric examples was used here to make it "easier" to see what happens/changes with that function call. I have also proposed some modifications to make it easier to understand (like plane name convention), and some code changes to probably make it faster and with less code.
  2. Used iamosaic, iaapplylut and ianormalize, in conjunction with adshow function. I have also made some calls to adshow using negative numbers. I was able to see that the adshow function uses only numbers ranged between 0 and 255 so, every number that is before or after that range, is truncated to a new number, for example, 256 becames 0. Some internal treatment (probably normalization or mod division?) is done inside of that function, but since no documentation is available, I cannot affirm that this is what really happens.