Generate Spoke Pattern image

Usage

Edit spoke_pattern.py using any text editor

then

giplt spoke_pattern.py <output image>

Description

This script can be used to create a “spoke pattern” image in the style of the ones featured in the following paper:

Downing K.H., Glaeser R.M., Restoration of weak phase-contrast images recorded with a high degree of defocus: The “twin image” problem associated with CTF correction, Ultramicroscopy 108 (2008), 921-928, PMID: 18508199, PMCID: PMC2694513

All tweakable parameters are defined in the first lines of the script, which can be edited before the script is executed. The script generates the image, opens a viewer to show it and saves it in the current folder with the name specified by the user. See figure 1 for a sample image.

Figure 1: Spoke pattern

Script

import math
import sys
from ost.img import *

# Tweakable parameters
wedge_start = 80 # Width of wedge at the narrower point
wedge_end = 2000 # Width of wedge at the broader point
number_of_bands = 11 # Number of alternating bands (must be odd)
size_of_image_x = 2000 # Size of the image (X)
size_of_image_y = 2000 # Size of the image (Y)
pixel_sampling = 1.0 * Units.A # Pixel width
threshold = 2.0 * Units.A # Threshold for low pass filtering
amplitude = 255 # amplitude of the obscillations

# Image is created
image=CreateImage(Size(size_of_image_x,size_of_image_y))
image.CenterSpatialOrigin()
image.SetSpatialSampling(pixel_sampling)

# For a cleaner code, 4 shorthands are defined
image_extent=image.GetExtent()
start_y=image_extent.GetStart()[1]
end_y=image_extent.GetEnd()[1]
start_x=image_extent.GetStart()[0]
end_x=image_extent.GetEnd()[0]

# Wedge is written in the image
for y in range (start_y,end_y+1):
  wedge_width=wedge_start+(y-start_y)*(wedge_end-wedge_start)/(end_y-start_y)
  for x in range (start_x,end_x+1):
    half_wedge=wedge_width/2.0
    outer_bands=(number_of_bands-1)/2
    factor=(0.5+outer_bands)*math.pi/half_wedge
    if float(x)>-half_wedge and float(x)<half_wedge:
      value=255*math.cos(float(x)*factor)
      image.SetReal(Point(x,y),value)

# Image is low-pass filtered
filter=alg.GaussianLowPassFilter(threshold)
image.ApplyIP(filter)

# Image is saved
ost.io.SaveImage(image,sys.argv[1])

# Viewer is launched to show the result
Viewer(image)