CAT scanner
From Noah.org
Contents |
CAT Scanning
I'd like to build a CAT Scanner. It turns out that the software isn't that hard. The algorithm is quite simple. Here are some sample images that illustrate the process. Since I don't yet have the hardware to generate x-ray slices of an object I have synthesized what an x-ray would yield if it scanned an object. These scans are then reassembled to yield the original target image.
In a real CAT Scan system the 1 dimensional slices would be taken from the horizontal row of a series of x-rays. In this demo I don't yet have the x-rays to work with so I synthesize the 1D bands from the target image that I want to regenerate. So given a target image I generate a series of 1D radial slices by rotating the target image and then averaging all values in the rows of the image. Then I rotate the slice back to the original angle.
Slices are synthesized from a 180 degree rotation of the target image.
The target and four 1-dimensional sections. The position of the sections corresponds to the angle of projection. | File:target.png | File:section 4 0.png |
| File:section 4 3.png | File:section 4 2.png | File:section 4 1.png |
Add (overlay) the four 1D radial sections together to get the CAT scan image at the right.
File:section 4 0.png+File:section 4 3.png+File:section 4 2.png+File:section 4 1.png = File:target 4 out.png
The resulting CAT scan is reminiscent of the target, but it is ambiguous. I found that you need at least 8 sections to get a recognizable image. The more sections you use the better. Here is a composite of 8 radial sections:
File:targetd animated.gif File:target 8 out.png
And finally, here is 32 radial sections. This seems to be a good number for an image this resolution:
This simple algorithm will work with very complex images. Given a photograph I synthesized 32 1-dimensional scans and then regenerated the photograph using the CAT algorithm.
A more complex example
File:target n.png File:target n out.png File:target n out enhanced.png
The real world of x-rays is not so simple. In the experiments above, I actually synthesize the radial sections by averaging all the pixels in the rows. This roughly approximates how an x-ray reveals average density of a line through the target. But in the real world there can be materials inside the target that are so dense that they totally block x-ray energy. This reconstruction technique assumes that every part of a target is at least somewhat transparent. If there are parts of a target that are totally opaque even to x-rays then this will result in ambiguous, hidden sections. These hidden sections not only hide what is inside of them, but they also cause shadows that distort areas outside.
The following target is similar to the one used before except that it now has a screen added. The algorithm that synthesizes the radial sections was modified to treat any red are as totally opaque. The result is that anything inside the cup shaped screen is totally hidden. The dense area also throws off the contrast so that it is difficult to see the notch at the top of the target, but you can more or less make it out.
File:targeto.png File:targeto out.png
Image processing
The image that results from the composite of the 1D sections has very low contrast. It is simple to expand the dynamic range of the image, but also note that the contrast is weighted towards the center. This is because the radial sections favor the center of the image. The center of the target has the most overlapping sections so the pixels near the center contribute more signal to the average. It's difficult to apply a uniform contrast enhancement over the entire image because the result will leave the edges too dark or the center too light. What is needed is contrast enhancement that will be weighted based on the distance from the center of the image.
Source code
This is a rough draft of the source code written in Python. It requires the PIL module for image manipulation.
Click here to download: cat_scanner.py
#!/usr/bin/env python
"""This demomstrates the CAT Scan algorithm. This takes an image and
synthesizes 1-Dimensional radial slices. These slices approximate what an x-ray
machine would see if it scanned the image. The 1D slices are then recombined to
show how the CAT Scan algorithm reconstructs the image from 1D slices.
Noah Spurrier
"""
import Image
import math
from pprint import pprint
# This is the name of the target PNG format image file that will be used
# to synthesize the radial sections. Give this name without the .PNG extension.
TARGET_NAME = 'targetc'
# This sets the number of radial sections to synthesize.
SECTIONS = 32
def enhance (value, x,y, cent_x, cent_y):
"""This enhanced the contrast of the image with a radial dampening factor towards the center.
In other words, the center of the image is less enhanced than the edge.
"""
xd = x - cent_x
yd = y - cent_y
distance = math.sqrt (xd * xd + yd * yd)
max_distance_center = max(cent_x, cent_y)
factor = 1.0 + (0.282828 * distance / max_distance_center)
return value * factor
def band_horizontal (im):
"""This takes an image and turns each horizontal row into a homogenous band
by averaging all pixels in the row. This creates a 1D section of the image.
The final image is square, but all pixels in a column are redundant. A
single 1 pixel high image could be used, but it is easier to work with
square images. """
im_new = Image.new ("L", im.size, 255)
width = im.size[0]
height = im.size[1]
s = list(im.getdata())
k = 0
for y in range (0, height):
total = 0
for x in range (0, width):
total = total + s[y*height+x]
gray = total / width
for x in range (0, width):
im_new.putpixel((x,y), gray)
return im_new
def average_pixels (bands, (x, y)):
"""This takes a set of images and returns the average pixel value averaged
over the same (x,y) cordinate in each image. This is used to overlay a set
of images into one image. This assume the set of images are the same size.
"""
a = 0
for b in bands:
p = b.getpixel((x,y))
a = a + p
l = len(bands)
value = a / l
value = enhance (value, x,y, bands[0].size[0]/2, bands[0].size[1]/2)
return int(value)
print 'step 1 - Synthesizing 1D radial bands'
im = Image.open(TARGET_NAME+'.png')
# If not gray scale, then convert to gray scale.
if im.mode != 'L':
im = im.convert('L')
for a in range (SECTIONS):
print '\tBand #%0d' % a
print '\t\tangle:', a*(180.0/SECTIONS)
im2 = im.rotate(a*(180.0/SECTIONS), Image.BICUBIC)
im3 = band_horizontal(im2)
im4 = im3.rotate(a*(-180.0/SECTIONS), Image.BICUBIC)
im4.save (TARGET_NAME+'_%02d.png'%a)
print 'step 2 - Collecting 1D radial bands'
imc = Image.new ("L", im.size, (255))
imband=[]
for a in range (SECTIONS):
imband.append(Image.open(TARGET_NAME+'_%02d.png' % a))
print 'step 3 - Merging 1D radial bands into image.'
for y in range (im.size[0]):
print '\tRow #%0d' % y
for x in range (im.size[1]):
av = average_pixels (imband, (x,y))
imc.putpixel ((x,y), av)
imc.save (TARGET_NAME+'_out.png')
