content: /*
   * Copyright (C) 2014 Madgalena Herbst
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
   */
  package edu.stanford.rsl.tutorial.noncircularfov;
  
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators;
  
  /**
   * Implementation of an algorithm that completes a given sinogram by using
   * alpha2 = -alpha1 and beta2 = beta1 + pi + 2 alpha1. One method completes the
   * sinogram and another method shows the double acquired data. Another method
   * finds the minimal complete set for a given sinogram and detector size.
   * 
   * @author Magdalena Herbst
   * 
   */
  
  public class Correspondences {
  
  	// parameters of the scanning geomtry
  	private double focalLength;
  	private double deltaBeta, deltaT;
  
  	int maxTIndex, maxBetaIndex;
  
  	/**
  	 * Focal length is defined in the constructor.
  	 * 
  	 * @param focalLength
  	 */
  
  	public Correspondences(double focalLength, Grid2D sino) {
  		this.focalLength = focalLength;
  
  		// set other parameters according to the sinogram.
  		maxTIndex = sino.getWidth();
  		maxBetaIndex = sino.getHeight();
  		deltaT = sino.getSpacing()[0];
  		deltaBeta = sino.getSpacing()[1];
  	}
  
  	/**
  	 * Method that displays the double acquired data. Double acquired data is
  	 * colored gray, missing areas or the result of zero-line-integrals are
  	 * colored black. The given sinogram has to be over a rotation of 360
  	 * degrees (it is the ground truth). The mask is a binary mask that is 1 at
  	 * the positions where data was acquired. The mask and the sinogram must
  	 * have the same dimensions.
  	 * 
  	 * @param mask
  	 * @param sino
  	 * @return the redundant areas
  	 * @throws IllegalArgumentException
  	 */
  	public Grid2D findRedundancies(Grid2D mask, Grid2D sino)
  			throws IllegalArgumentException {
  
  		if (mask.getHeight() != sino.getHeight()
  				|| mask.getWidth() != sino.getWidth()) {
  			throw new IllegalArgumentException("Dimensions must agree!");
  		}
  
  		double mid = maxTIndex / 2.0;
  		Grid2D erg = new Grid2D(maxTIndex, maxBetaIndex);
  		erg.setSpacing(deltaT, deltaBeta);
  
  		for (int y = 0; y < maxBetaIndex; y++) {
  			for (int x = 0; x < maxTIndex; x++) {
  
  				if (mask.getAtIndex(x, y) > 0 && sino.getAtIndex(x, y) > 0) {
  					int xc = (int) (maxTIndex - 1) - x;
  					double alpha = (Math.atan2((x - mid) * deltaT, focalLength));
  					double betac = y * deltaBeta + Math.PI + 2 * alpha;
  					betac = betac % (2 * Math.PI);
  					if (erg.getAtIndex(x, y) > 0) {
  						erg.setAtIndex(x, y, .5f);
  						erg.setAtIndex(xc, (int) (betac / deltaBeta), .5f);
  					} else {
  						erg.setAtIndex(x, y, 1.f);
  						erg.setAtIndex(xc, (int) (betac / deltaBeta), 1.f);
  					}
  
  				}
  			}
  		}
  		return erg;
  	}
  
  	/**
  	 * Methods that completes a given sinogram. The given sinogram has to be
  	 * over a rotation of 360 degrees (it is the ground truth). The mask is a
  	 * binary mask that is 1 at the positions where data was acquired. The mask
  	 * and the sinogram must have the same dimensions.
  	 * 
  	 * @param mask
  	 * @param sino
  	 * @return the completed sinogram
  	 * @throws IllegalArgumentException
  	 */
  	public Grid2D findCorrespondences(Grid2D mask, Grid2D sino)
  			throws IllegalArgumentException {
  
  		if ((mask.getHeight() != sino.getHeight())
  				|| (mask.getWidth() != sino.getWidth())) {
  			throw new IllegalArgumentException("Dimensions must agree!");
  		}
  
  		maxTIndex = sino.getWidth();
  		maxBetaIndex = sino.getHeight();
  		deltaT = sino.getSpacing()[0];
  		deltaBeta = sino.getSpacing()[1];
  
  		double mid = maxTIndex / 2.0;
  		Grid2D erg = new Grid2D(maxTIndex, maxBetaIndex);
  		erg.setSpacing(deltaT, deltaBeta);
  
  		// first set acquired values
  		for (int y = 0; y < maxBetaIndex; y++) {
  			for (int x = 0; x < maxTIndex; x++) {
  				if (mask.getAtIndex(x, y) > 0) {
  					erg.setAtIndex(x, y, sino.getAtIndex(x, y));
  				}
  			}
  		}
  
  		// then compute the missing values using the corresponding values
  		for (int y = 0; y < maxBetaIndex; y++) {
  			for (int x = 0; x < maxTIndex; x++) {
  				if (mask.getAtIndex(x, y) == 0) {
  					double xc = (maxTIndex - 1) - x;
  					double alpha = (Math.atan2((x - mid) * deltaT, focalLength));
  					double yc = y * deltaBeta + Math.PI + 2 * alpha;
  					yc = yc % (2 * Math.PI);
  
  					float val = InterpolationOperators.interpolateLinear(erg,
  							xc, yc / deltaBeta);
  					erg.setAtIndex(x, y, val);
  				}
  			}
  		}
  
  		return erg;
  	}
  
  	/**
  	 * Method that computes the missing pixels that result if segments of a
  	 * given sinogram are completed with findCorrespondences(...). The detector
  	 * follows the left boundary.
  	 * 
  	 * @param sino
  	 *            ground truth sinogram, has to be over a complete rotation.
  	 * @param start
  	 *            index of starting point
  	 * @param rotationLength
  	 *            in pixels
  	 * @param size
  	 *            the detector size in pixels
  	 * @return number of missing pixels
  	 */
  	public int missing(Grid2D sino, int start, int rotationLength, int size) {
  
  		int width = sino.getWidth();
  		int height = sino.getHeight();
  
  		// Create mask. Detector follows the left boundary.
  		Grid2D mask = new Grid2D(width, height);
  		for (int y = start; y < start + rotationLength; y++) {
  			for (int x = 0; x < width; x++) {
  
  				// y % height to simulate a never-ending sinogram
  				//if (sino.getAtIndex(x, y % height) > 0) {
  				if ((x>150)&&(x<150+size)) {
  					for (int i = 0; i < size; i++) {
  						if (x + i < width) {
  							mask.setAtIndex(x + i, y % height, 1.f);
  						}
  					}
  
  					break;
  				}
  
  			}
  		}
  
  		// Complete the sinogram.
  		Grid2D res = findCorrespondences(mask, sino);
  
  
  		int missing = 0;
  		for (int y = 0; y < height; y++) {
  			for (int x = 0; x < width; x++) {
  
  				// count the pixels that are zero in the completed image but not
  				// in the ground truth image
  				if (sino.getAtIndex(x, y) > 0 && res.getAtIndex(x, y) == 0) {
  					missing++;
  				}
  			}
  		}
  
  		return missing;
  	}
  
  	/**
  	 * Method that performs a grid search to find the minimal complete dataset.
  	 * A better way to choose the minimal set may be found, the current one is
  	 * not perfect. Note that this method may have a runtime up to 15 minutes
  	 * because it is not optimized.
  	 * 
  	 * @param sino
  	 * @param size
  	 *            size of the detector in pixel
  	 * @return result[0] = required rotation, result[1] = corresp. starting
  	 *         point, result[2] = corresp. missing pixel
  	 */
  	public int[] evaluateNumer(Grid2D sino, int size) {
  
  		int a = sino.getHeight() / 2; // start at 180 degrees
  		int b = sino.getHeight();
  
  		int[] werteStart = new int[b - a];
  		int[] werteMiss = new int[b - a];
  		int[] res = new int[3];
  
  		for (int rot = a; rot < b; rot++) {
  
  			// find best starting point for current rotationLength
  			int minMiss = Integer.MAX_VALUE;
  			int minMissStart = 0;
  
  			for (int start = 0; start < 360; start++) {
  				int miss = missing(sino, start, rot, size);
  
  				if (miss < minMiss) {
  					minMiss = miss;
  					minMissStart = start;
  				}
  			}
  
  			werteStart[rot - a] = minMissStart;
  			werteMiss[rot - a] = minMiss;
  
  		}
  
  		int min_idx = 0;
  
  		// uncomment the output to see the missing pixel for the other
  		// settings. The chosen minimal set might not be the best.
  
  		// System.out.println("rot\t" + "start\t" + "missing");
  		for (int i = 0; i < werteStart.length; i++) {
  			// System.out.println((i + a) + "\t" + werteStart[i] + "\t" +
  			// werteMiss[i]);
  
  			if (i != 0) {
  				// choose the setting where the number of missing pixel does not
  				// change any more (or only a little)
  				int change = 4; // or 0, 1, 2, 3,... did not find the perfect
  				// solution yet.
  				if (werteMiss[i] == werteMiss[i - 1] - change) {
  					min_idx = i;
  				}
  			}
  		}
  
  		res[0] = min_idx + a;
  		res[1] = werteStart[min_idx];
  		res[2] = werteMiss[min_idx];
  		return res;
  	}
  
  }
  