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;
  
  public class BinaryWeights_erodeByOne extends Grid2D {
  
  	private final double focalLength;
  	private final double maxT;
  	private final double deltaT, dLambda;
  	private final int maxTIndex, maxLambdaIndex;
  	private Grid2D binaryMask;
  
  	public BinaryWeights_erodeByOne(final double focalLength, final double maxT,
  			final double deltaT, double maxLambda, double dLambda) {
  		// Call constructor from superclass
  		super(((int) Math.round(maxT / deltaT)), ((int) Math.round(maxLambda / dLambda) + 1));
  		
  		// Initialize parameters
  		this.focalLength = focalLength;
  		this.maxT = maxT;
  		this.deltaT = deltaT;
  		this.dLambda = dLambda;
  		this.maxLambdaIndex = (int) (Math.round(maxLambda / dLambda)) + 1;
  		this.maxTIndex = (int) Math.round(maxT / deltaT);
  
  		binaryMask = new Grid2D(this);
  		createMask();
  		createWeights();
  
  		// Correct for scaling due to varying angle lambda
  		NumericPointwiseOperators.multiplyBy(this, (float) (maxLambda / (Math.PI)));
  	}
  
  	private void createMask() {
  		double lambda, delta;
  		// iterate over the detector elements
  		for (int t = 0; t < maxTIndex; ++t) {
  			// compute delta of the current ray (detector element)
  			delta = Math.atan((t * deltaT - maxT / 2.d + 0.5 * deltaT)
  					/ focalLength);
  
  			// iterate over the projection angles
  			for (int b = 0; b < maxLambdaIndex; ++b) {
  				// compute the current lambda angle
  				lambda = b * dLambda;
  
  				// Default weight is 1
  				binaryMask.setAtIndex(maxTIndex - t - 1, b, 1.0f);
  				this.setAtIndex(maxTIndex - t - 1, b, 1.0f);
  
  				// First case: Handles weights for redundancies at the end of
  				// the scan
  				// Set weights to zero (Simulates the collimator)
  				if (lambda >= (Math.PI + 2 * delta)) {
  					binaryMask.setAtIndex(maxTIndex - t - 1, b, 0.f);
  					this.setAtIndex(maxTIndex - t - 1, b, 0.f);
  				}
  			}
  		}
  
  	}
  
  	private void createWeights() {
  		Grid2D tmp = new Grid2D(this);
  		NumericPointwiseOperators.fill(tmp, -100.f);
  
  		// basically an dilation step!
  		// iterate over the projection angles
  		for (int b = 0; b < maxLambdaIndex; ++b) {
  			// iterate over the detector elements
  			for (int t = 0; t < maxTIndex; ++t) {
  				if (t > 0) {
  					if (this.getAtIndex(t, b) == 1 || this.getAtIndex(t-1, b) == 1) {
  						tmp.setAtIndex(t, b, 1.f);
  					} 
  					else {
  						tmp.setAtIndex(t, b, 0.f);
  					}
  				}
  				if (t == 0)
  				{
  					tmp.setAtIndex(t, b, this.getAtIndex(t, b));
  				}
  			}
  		}
  		// iterate over the projection angles
  		for (int b = 0; b < maxLambdaIndex; ++b) {
  			// iterate over the detector elements
  			for (int t = 0; t < maxTIndex; ++t) {
  				this.setAtIndex(t, b, tmp.getAtIndex(t, b));
  			}
  		}
  
  	}
  
  	public Grid2D getBinaryMask() {
  		return binaryMask;
  	}
  
  	public void applyToGrid(Grid2D sino) {
  		NumericPointwiseOperators.multiplyBy(sino, this);
  	}
  
  	public static void main(String[] args) {
  		// fan beam bp parameters
  
  		double maxT = 400;
  		double deltaT = 1.d;
  		// set focal length according to the fan angle
  		double focalLength = (maxT / 2.0 - 0.5)
  				/ Math.tan(30.0 * Math.PI / 180.0);
  
  		new ImageJ();
  		Grid3D g = new Grid3D((int)maxT, 181, 360); 
  		for (int i = 0; i < 360; ++i) { 
  			double maxBeta = (double)(i+1) * Math.PI * 2.0 / 360.0;
  			double deltaBeta = maxBeta / 180;
  
  			BinaryWeights p = new BinaryWeights(focalLength, maxT, deltaT, maxBeta, deltaBeta); 
  			g.setSubGrid(i, p); 
  			//g.setSliceLabel("MaxBeta: " + Double.toString(maxBeta*180/Math.PI), i+1); }
  		}
  		g.show();
  
  
  		double maxBeta = Math.PI;
  		double deltaBeta = maxBeta / 180;
  
  		BinaryWeights_erodeByOne p = new BinaryWeights_erodeByOne(focalLength, maxT, deltaT, maxBeta,
  				deltaBeta);
  		p.show("Weights");
  		p.getBinaryMask().show("Mask");
  
  	}
  
  }
  /*
   * Copyright (C) 2010-2014 Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */