Pylinac's main purpose is to make doing routine quality assurance easier and as automatic as possible. The main modules can be utilized in only a few lines of code and does a complete analysis of your QA images. However, when researching new tools, comparing specific algorithms, or doing custom testing, the default tools aren't appropriate. Pylinac's modules are built on simpler tools that can be used directly or easily extended. In this tutorial, we'll introduce some of pylinac's potential as a toolkit for customizing or creating new tests.
This is also viewable as a Jupyter notebook <http://nbviewer.ipython.org/github/jrkerns/pylinac/blob/master/docs/source/pylinac_core_hacking.ipynb>_.
In many applications of QA, the physicist has images acquired via the EPID, kV imager, or scanned films. The image module provides easy tools for loading, manipulating, and converting these images. For this tutorial we'll use a few demo images and showcase the various methods available.
There are a few simple tips for using the image module:
load, load_url, or load_multiples.Dicom, File, and ArrayThe docs for the image API are here.
First, our imports for the whole tutorial:
In [1]:
%matplotlib inline
from urllib.request import urlretrieve
import matplotlib.pyplot as plt
import numpy as np
from pylinac.core import image
# pylinac demo images' URL
PF_URL = 'https://s3.amazonaws.com/pylinac/EPID-PF-LR.dcm'
STAR_URL = 'https://s3.amazonaws.com/pylinac/starshot.tif'
# local downloaded images
PF_FILE, _ = urlretrieve(PF_URL)
STAR_FILE, _ = urlretrieve(STAR_URL)
# sample numpy array
ARR = np.arange(36).reshape((6,6))
Let's load the picket fence demo image:
In [2]:
pfimg = image.load_url(PF_URL)
type(pfimg)
Out[2]:
We can load the local file just as easily:
In [4]:
pfimg = image.load(PF_FILE)
type(pfimg)
Out[4]:
The load() function can take an string pointing to a file, a data stream, or a numpy array:
In [6]:
arr = np.arange(36).reshape((6,6))
img2 = image.load(arr) # .pixel_array is a numpy array of DICOM file
type(img2)
Out[6]:
Additionally, multiple images can be loaded and superimposed. For example multiple collimator star shots. Just pass a list of the images. In this case I'm passing the same image twice, but the point is this can be done with any set of images:
In [9]:
superimposed_img = image.load_multiples([STAR_FILE, STAR_FILE])
While the load() function will always do smart image inference for you, if you already know the file type you can instantiate directly. Furthermore, keyword arguments can be passed to FileImage and ArrayImage if they are known.
In [11]:
dcm_img = image.DicomImage(PF_FILE)
arr_img = image.ArrayImage(ARR, dpi=30, sid=500)
file_img = image.FileImage(STAR_FILE, dpi=50, sid=1000)
file_img2 = image.load(STAR_FILE, dpi=30, sid=500)
type(file_img) == type(file_img2)
Out[11]:
Images can be easily plotted using the plot() method:
In [12]:
dcm_img.plot()
Out[12]:
The plot can also be passed to an existing matplotlib axes:
In [13]:
fig, ax = plt.subplots()
dcm_img.plot(ax=ax)
Out[13]:
The image classes contain a number of useful attributes for analyzing and describing the data:
In [14]:
dcm_img.dpi
Out[14]:
In [15]:
dcm_img.dpmm
Out[15]:
In [16]:
dcm_img.sid
Out[16]:
In [17]:
dcm_img.shape
Out[17]:
There are also attributes that are useful in radiation therapy:
In [18]:
dcm_img.cax # the beam CAX
Out[18]:
In [19]:
dcm_img.center # the center location of the image
Out[19]:
Above, the values are the same because the EPID was not translated, which would move the CAX but not the image center.
The image values can also be sampled by slicing and indexing:
In [20]:
dcm_img[12, 60]
Out[20]:
In [21]:
dcm_img[:100, 82]
Out[21]:
Now the really fun stuff! There are many methods available to manipulate the data.
First, let's smooth the data:
In [24]:
file_img.filter(kind='median', size=3)
file_img.plot()
Out[24]:
Sometimes starshots from scanned film have edges that are very high or low value (corners of the film can be bent or rounded). We can easily trim the edges:
In [25]:
file_img.remove_edges(pixels=100)
file_img.plot()
Out[25]:
The data can also be explicitly inverted (EPID images oftentimes need this), or rolled on an axis:
In [26]:
file_img.invert()
file_img.plot()
Out[26]:
In [28]:
file_img.roll(direction='y', amount=300)
file_img.plot()
Out[28]:
We can also rotate and resize the image:
In [29]:
file_img.rot90(n=1)
file_img.plot()
Out[29]:
In [30]:
file_img.resize(size=(1000, 1100))
file_img.plot()
Out[30]:
Scanned film values can be very high, even in low dose areas. We can thus "ground" the image so that the lowest value is zero; this will help us later on when detecting profiles.
In [31]:
np.min(file_img) # before grounding
Out[31]:
In [33]:
file_img.ground()
np.min(file_img) # after grounding
Out[33]:
In [34]:
file_img.plot()
Out[34]:
We can also apply a high-pass filter to the image:
In [35]:
thresh_val = np.percentile(file_img, 95)
file_img.threshold(thresh_val)
file_img.plot()
Out[35]:
The image can also be converted to binary, which can be used later for ROI detection. Note that unlike any other method, this returns a new ArrayImage (hinted by the as_)...
In [36]:
new_img = file_img.as_binary(thresh_val)
new_img.plot()
Out[36]:
...and leaves the original unchanged.
In [37]:
file_img.plot()
Out[37]:
Physicists often need to evalute a profile, perhaps from a linac beam EPID image, or some fluence profile. The profile module allows the physicist to find peaks in a 1D array and determine beam profile information (FWHM, penumbra, etc). There are two main profile classes:
SingleProfile - This class is for profiles with a single peak; e.g. an open beam delivered to a film or EPID. The main goal of this class is to describe the profile (FWHM, penumbra, etc). MultiProfile - This class is for profiles with multiple peaks. The main goal of this class is to find the peak or valley locations. A MultiProfile can be broken down into SingleProfiles.CircleProfile - A MultiProfile, but in the shape of a circle. CollapsedCircleProfile - A CircleProfile that "collapses" a thick ring of pixel data to create a 1D profile.The profile API docs are here.
For this demonstration we'll find some peaks and then determine profile information about one of those peaks. Let's use the starshot demo image since it contains all the types of profiles:
In [38]:
from pylinac.core.profile import SingleProfile, MultiProfile, CircleProfile, CollapsedCircleProfile
star_img = image.load_url(STAR_URL)
Let's start by sampling one row from the starshot image:
In [39]:
row = star_img[800, :]
plt.plot(row)
Out[39]:
So, judging by the profile, it needs to be filtered for the spurious signals, it has multiple peaks, and it's upside down.
Let's make a MultiProfile and clean up the data.
In [40]:
mprof = MultiProfile(row)
mprof.plot()
First, let's invert it so that pixel value increases with dose.
In [41]:
mprof.invert()
mprof.plot()
We've loaded the profile and inverted it; let's run a filter over it.
In [42]:
mprof.filter(size=6)
mprof.plot()
The profile could probably be filtered more since there's still a few spurious signals, but this will work nicely for our demonstration.
First, we want to find the peak locations:
In [43]:
mprof.find_peaks()
Out[43]:
The method has found the 3 major peaks of the profile. Note that there are actually 5 peaks if we count the spurious signals near indices 800 and 1200.
For fun, let's see if we can detect these peaks. We can change the parameters to find_peaks() to optimize our search.
In [44]:
mprof.find_peaks(threshold=0.1)
Out[44]:
By lowering the peak height threshold we've found another peak; but the peak near 1200 wasn't found. What gives?
The find_peaks() method also eliminates peaks that are too close to one another. We can change that:
In [45]:
mprof.find_peaks(threshold=0.1, min_distance=0.02)
Out[45]:
By changing the minimum distance peaks must be from each other, we've found the other peak.
But, let's say we need to use these settings for whatever reason. We can additionally limit the number of peaks using max_number.
In [46]:
mprof.find_peaks(threshold=0.1, min_distance=0.02, max_number=3)
Out[46]:
Now, we can visualize where these peaks are by using the plot() method, which shows the peaks if we've searched for them; note the green dots at the detected peak locations.
In [47]:
mprof.plot()
We can also search a given portion of the region; for example if we only wanted to detect peaks in the first half of the profile we can easily add a search_region. Note that the last peak was not detected.
In [48]:
mprof.find_peaks(search_region=(0, 0.6)) # search the left 60% of the profile
Out[48]:
We can search not simply for the max value peaks, but for the FWHM peaks. Note that these values are slightly different than the max value peaks we found earlier.
In [49]:
mprof.find_fwxm_peaks(x=50) # 50 is 50% height
Out[49]:
Finally, we can subdivide the profile into SingleProfile's to further describe single peaks:
In [50]:
single_profiles = mprof.subdivide() # returns a list of SingleProfile's
SingleProfiles are useful to describe profiles with a single peak. It can describe the FWXM (X=any height), penumbra on each side, field width, and calculations of the field. Continuing from above:
In [51]:
sprof = single_profiles[0]
sprof.plot()
The multiprofile has been cut into multiple single profiles, of which this is the first.
Let's first find the FWHM, and the center of the FWHM:
In [52]:
sprof.fwxm(x=50)
Out[52]:
In [53]:
sprof.fwxm_center(x=50)
Out[53]:
Note that this is the same value as the first FWHM peak value we found in the MultiProfile.
We can now find the penumbra values:
In [54]:
sprof.penumbra_width(side='left', upper=80, lower=20)
Out[54]:
In [55]:
sprof.penumbra_width(upper=90, lower=10) # default is average of both sides
Out[55]:
The careful reader will notice that the profiles, since we created them, has not had a minimum value of 0. Normally, this would cause problems and sometimes it does, but pylinac normalizes the FWXM and penumbra search. However, just to make sure all is well we can easily shift the values so that the lower bound is 0:
In [56]:
sprof.ground()
sprof.plot()
In [57]:
sprof.penumbra_width(upper=90, lower=10)
Out[57]:
The average penumbra is the same as we found earlier.
We can also normalize and stretch the profile values. Let's first get the original maximum value so we know what we need to restore the profile:
In [58]:
np.max(sprof)
Out[58]:
In [59]:
sprof.normalize()
sprof.plot()
We can also stretch the values to be bounded by any values we please:
In [60]:
sprof.stretch(min=30, max=284)
sprof.plot()
Let's restore our profile based on the earlier max and min values:
In [61]:
sprof.stretch(max=47, min=0)
sprof.plot()
Circular profiles are useful for concentric profiles; starshots are great examples. Let's explore the two circular profiles as they relate to a starshot image.
Let's once again load the starshot demo image:
In [62]:
star_img.plot()
Out[62]:
We saw above from the profile that the image is actually inverted (pixel value decreases with dose), so we need to invert the image:
In [63]:
star_img.invert()
To use a CircleProfile we need to specify the center and radius of the circle, as well as the array to operate over. The approximate center of the starshot is (1450, 1250). Let's also use a radius of 300 pixels.
In [64]:
cprof = CircleProfile(center=(1250, 1450), radius=300, image_array=star_img) # center is given as (x, y)
cprof.plot()
We have a nice profile showing the starshot peaks. It appears we've cut one of the peaks in half though; this is because the profile starts at 0 radians on the unit circle (directly to the right) and there is a starshot line right at 0. We can also change the direction of the profile from the default of counter-clockwise to clockwise:
In [65]:
cprof = CircleProfile(center=(1250, 1450), radius=300, image_array=star_img, start_angle=0.2, ccw=False)
cprof.plot()
Alternatively, we can roll the profile directly:
In [66]:
cprof.roll(amount=50)
cprof.plot()
Now, because CircleProfile is a subclass of MultiProfile we can search for the peaks:
In [67]:
cprof.find_peaks()
cprof.plot()
The profile is 1D, but was derived from a circular sampling. How do we know what the locations of the sampling is? We have x and y attributes:
In [68]:
cprof.x_locations
Out[68]:
In [69]:
cprof.y_locations
Out[69]:
We can also add this profile to a plot to show where it is:
In [70]:
ax = star_img.plot(show=False)
cprof.plot2axes(ax)
Looks good! Now let's take a CollapsedCircleProfile. The advantage of this class is that a thick ring is sampled, which averages the pixel values. Thus, noise and spurious signals are reduced. Beyond CircleProfile there are 2 more keyword arguments:
In [71]:
ccprof = CollapsedCircleProfile(center=(1250, 1450), radius=300, image_array=star_img, width_ratio=0.2, num_profiles=10)
ccprof.plot()
Note that this profile looks smoothed; this comes from averaging over the num_profiles within the ring. The width_ratio is a function of the radius, so in this case the actual ring width is 300*0.2 = 60 pixels, and 10 equally distributed profiles are taken within that ring.
Let's find the peaks and then plot the ring to the starshot image:
In [72]:
ccprof.find_peaks()
ccprof.plot()
In [73]:
ax = star_img.plot(show=False)
ccprof.plot2axes(ax)
We now have a good start on a starshot algorithm!
In [ ]: