In [ ]:
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rc('image', aspect='equal')
%matplotlib inline
import SimpleITK as sitk
# Download data to work on
from downloaddata import fetch_data as fdata
Let us begin by developing a convenient method for displaying images in our notebooks.
In [ ]:
img = sitk.GaussianSource(size=[64]*2)
plt.imshow(sitk.GetArrayFromImage(img))
In [ ]:
img = sitk.GaborSource(size=[64]*2, frequency=.03)
plt.imshow(sitk.GetArrayFromImage(img))
In [ ]:
def myshow(img):
nda = sitk.GetArrayFromImage(img)
plt.imshow(nda)
myshow(img)
If you are familiar with numpy, sliced index then this should be cake for the SimpleITK image. The Python standard slice interface for 1-D object:
Operation | Result |
d[i] | ith item of d, starting index 0 |
d[i:j] | slice of d from i to j |
d[i:j:k] | slice of d from i to j with step k |
With this convient syntax many basic tasks can be easily done.
In [ ]:
img[24,24]
In [ ]:
myshow(img[16:48,:])
In [ ]:
myshow(img[:,16:-16])
In [ ]:
myshow(img[:32,:32])
In [ ]:
img_corner = img[:32,:32]
myshow(img_corner)
In [ ]:
myshow(img_corner[::-1,:])
In [ ]:
myshow(sitk.Tile(img_corner, img_corner[::-1,::],img_corner[::,::-1],img_corner[::-1,::-1], [2,2]))
In [ ]:
img = sitk.GaborSource(size=[64]*3, frequency=0.05)
# Why does this produce an error?
myshow(img)
In [ ]:
myshow(img[:,:,32])
In [ ]:
myshow(img[16,:,:])
In [ ]:
myshow(img[:,::3,32])
Most python mathematical operators are overloaded to call the SimpleITK filter which does that same operation on a per-pixel basis. They can operate on a two images or an image and a scalar.
If two images are used then both must have the same pixel type. The output image type is ussually the same.
As these operators basically call ITK filter, which just use raw C++ operators, care must be taked to prevent overflow, and divide by zero etc.
Operators |
+ |
- |
\* |
/ |
// |
** |
In [ ]:
img = sitk.ReadImage(fdata("cthead1.png"))
img = sitk.Cast(img,sitk.sitkFloat32)
myshow(img)
img[150,150]
In [ ]:
timg = img**2
myshow(timg)
timg[150,150]
All three Python division operators are implemented __floordiv__
, __truediv__
, and __div__
.
The true division's output is a double pixle type.
See PEP 238 to see why Python changed the division operator in Python 3.
In [ ]:
img = sitk.ReadImage(fdata("cthead1.png"))
myshow(img)
In [ ]:
img = sitk.ReadImage(fdata("cthead1.png"))
myshow(img)
In [ ]:
myshow(img>90)
In [ ]:
myshow(img>150)
In [ ]:
myshow((img>90)+(img>150))