; b00.radon
; Copyright 2002 Hanspeter Niederstrasser
; version 1.0 2002-07-29

; This script is released under a modified version of the GNU General
; Public License (GPL). The only difference from the GPL is that if you
; use the output from this script in any forum (such as a, but not limited
; to, presentation, book, poster, web site, etc) you must credit the 
; Electron Microscopy for Dummies web site as the source.

;    This program is free software; you can redistribute it and/or modify
;    it under the terms of the GNU General Public License as published by
;    the Free Software Foundation; either version 2 of the License, or
;    any later version.
;   
;    This program is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;    GNU General Public License for more details.
;   
;    You should have received a copy of the GNU General Public License
;    along with this program; if not, write to the Free Software
;    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


; Makes 1D projections from a 2D image (using the radon program)
; extracts each individual line from the sinogram
; makes a 2D 'stretched' image out of each 1D image, 
; rotates it according to its originial projection angle,
; and averages the rotated projections to reconstruct the original model.

; Can be modified to generate different # of projections and use different
; a subset of them when reconstructing the model.  This can be used to show
; how smearing back projections and adding them can lead to a reconstruction.


; Useful constants that user should modify
; x11 = degrees between original projections. Safest to use factors of 180.
; x12 = dimensions of original image
; x13 = radius of model.  Keep it tight to the real object used.

x11=5
x12=500
x13=111

; constants used by program.  Generally should not be modified by user.
; x51 = # of projections to be generated (based on x11)
; x54 = radius of model (x13) with 25% padding
; x55 = width of a line to be chopped off from sinogram

x51=INT(180/x11)     ; # of projections that will be made, rounded down
x54=INT(x13*1.25)+1  ; padded radius of model
x55=1                ; width of line chopped from sinogram

; make the radon transform of the input image
RM 2DN
model    ; input image
x11      ; angular separation between projections
radon001 ; name of resulting radon image (sinogram)
x12      ; width of radon transform output
x54      ; radius of circle to transform
(0,0)    ; offset from center.
n

; the radon image *may* have very strange min and max values.
; This normalizes them to grey values of 0-255.

;FI X21,X22   ; file info. operation
;radon001     ; input file
;8,7          ; min and max location in header.  Applied to x21 and x22
;
;AR           ; arithmetic operation
;radon001     ; input file
;radon002     ; output file
;255*(P1-X21)/(X22-X21)

; Cut out each individual line of the sinogram and stick it into its 
; individual image file.  
; Interpolate the line image to stretch into a square image.
; Rotate the stretched image to it's corresponding angle of origin.

; The loop defaults to make each step in memory, but can easily be modified
; to save to a physical file for presentational purposes.

DO LB1 x60=1,x51

WI
radon001                  ; source sinogram
_1 ; line/line{***x60}    ; output 1D image
x12,x55                   ; cols,rows of windowed image
x55,x60                   ; across down coordinate for cut

IP
_1 ; line/line{***x60}    ; input 1D image
_2 ; smear/smear{***x60}  ; output 2D smear of 1D image
x12,x12                   ; width,height of 2D output

x56=x11*(x60-1)           ; find the angle for the current projection
RT
_2 ; smear/smear{***x60}  ; input 2D smear of 1D image
rot/rot{***x60}           ; rotated output
x56                       ; rotate by x11*(x60-1).  0deg is the first rotation.

CP TO TIFF
rot/rot{***x60}           ; rotated image
tif/tif{***x60}.tif       ; outputed tif

LB1

; Average the rotated images in various steps to generate reconstructions with 
; different number of projections
;DO LB2 x61=1,x51

AS R
rot/rot***
1-x51                     ; the individual images that will be averaged
a                         ; average (a)ll the images. Can also be (e)ven/(o)dd
avg001                    ; output location for average
_1                        ; output location for variance (not used here and discarded)

;LB2

; Make a montage of each step of the back projection step.  Each frame includes 
; a model made of one more projection than the frame before it.

; MN S
; montage is not yet implemented

EN D