content: package edu.stanford.rsl.tutorial.fan.redundancy;
  
  import ij.ImageJ;
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.Grid3D;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  
  // Variable according to Slaney 
  // Note: Parker weights are not defined for super-short scans (maxBeta < pi + 2*gammaM) !! Use Silver or Compensation weights instead!
  
  public class ParkerWeights extends Grid2D {
  
  	private final double gammaM;
  	private final int maxTIndex, maxBetaIndex;
  
  	public ParkerWeights(final double focalLength, final double maxT,
  			final double deltaT, double maxBeta, double deltaBeta) {
  		// Call constructor from superclass
  		super((int) Math.round(maxT / deltaT), (int)Math.round(maxBeta / deltaBeta) + 1);
  		
  		// Initialize parameters
  		this.maxBetaIndex = (int)(Math.round(maxBeta / deltaBeta)) + 1;
  		this.maxTIndex = (int) Math.round(maxT / deltaT);
  		this.gammaM = (double) (Math.atan((maxT / 2.d)/ focalLength) );
  		double beta, alpha;
  
  		// iterate over the detector elements
  		for (int t = 0; t < maxTIndex; ++t) {
  			// compute alpha of the current ray (detector element)
  			alpha = Math.atan((t * deltaT - maxT / 2.d + 0.5*deltaT) / focalLength);
  			
  			// iterate over the projection angles
  			for (int b = 0; b < maxBetaIndex; ++b) {
  				beta = b * deltaBeta;
  				
  				// Shift weights such that they are centered (Important for maxBeta < pi + 2 * gammaM)
  					beta += (Math.PI+2*gammaM-maxBeta)/2.0;
  				
  				// Adjust beta if out of range [0, 2*pi]
  				if (beta < 0) {
  					continue;
  				}
  				if (beta > Math.PI *2.d) {
  					continue;
  				}
  
  				// implement the conditions as described in Parker's paper
  				if (beta <= 2 * (gammaM - alpha)) {
  					double tmp = beta * Math.PI / 4.d / (gammaM - alpha);
  					float val = (float) Math.pow(Math.sin(tmp), 2.d);
  					
  					if (Double.isNaN(val)){
  						continue;
  					}
  					this.setAtIndex(t, b , val);
  
  				} else if (beta < Math.PI - 2.d * alpha) {
  					this.setAtIndex(t, b , 1);
  				}
  				else if (beta <= (Math.PI + 2.d * gammaM) + 1e-12) {
  					double tmp = (Math.PI / 4.d) * ( (Math.PI + 2.d*gammaM - beta) / (gammaM + alpha) );
  					float val = (float) Math.pow(Math.sin(tmp), 2.d);
  					if (Double.isNaN(val)){
  						continue;
  					}
  					this.setAtIndex(t, b , val);
  				}
  			}
  		}
  		
  		// Correct for scaling due to varying angle
  		NumericPointwiseOperators.multiplyBy(this, (float)( maxBeta / (Math.PI)));
  
  	}
  
  	public void applyToGrid(Grid2D sino) {
  		NumericPointwiseOperators.multiplyBy(sino, this);
  	}
  
  	public static void main (String [] args){
  		//fan beam bp parameters
  
  		double maxT = 285;
  		//double gammaM =Math.atan2(maxT / 2.f - 0.5, focalLength) * 2.0;
  		double fan = 10.0*Math.PI/180.0;
  		double focalLength = (maxT/2.0-0.5)/Math.tan(fan);
  		new ImageJ();
  		double deltaT = 1.d;
  		
  		Grid3D g = new Grid3D((int)maxT, 181, 360, false);
  		for (int i = 0; i < 360; ++i)
  		{
  			double maxBeta =  (double)(i+1) * Math.PI * 2.0 / 360.0;
  			double deltaBeta = maxBeta / 180;
  			
  			ParkerWeights p = new ParkerWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta);
  			//g.setSliceLabel("MaxBeta: " + Double.toString(maxBeta*180/Math.PI), i+1);
  			g.setSubGrid(i, p);
  			
  		}
  		
  		g.show();
  
  	}
  
  }
  /*
   * Copyright (C) 2010-2014 Martin Berger
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */