/*******************************************************************************
 * Copyright (c) 2010 Attila Kovacs <attila_kovacs@post.harvard.edu>.
 * All rights reserved. 
 * 
 * This file is part of crush.
 * 
 *     crush 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 3 of the License, or
 *     (at your option) any later version.
 * 
 *     crush 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 crush.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Contributors:
 *     Attila Kovacs <attila_kovacs@post.harvard.edu> - initial API and implementation
 ******************************************************************************/
// Copyright (c) 2009 Attila Kovacs 

package crush.sourcemodel;

import java.util.*;

import nom.tam.fits.*;
import nom.tam.util.*;

import util.*;

import java.io.*;

import crush.*;


public class Image implements Cloneable {
	public double[][] data;
	public int[][] flag;
	protected double[][] beam;
	
	public Unit unit = Unit.unity;
	public Instrument<?> instrument;

	public SphericalProjection projection;
	public Vector2D refIndex = new Vector2D();
	public double delta;
	
	public double smoothFWHM;  
	public double extFilterFWHM = Double.NaN;
	public double correctingFWHM = Double.NaN;	
	
	protected Unit dataUnit = new Unit("d.u.", 1.0);
	protected Unit dataUnitPerBeam = new Unit("d.u./beam");
	
	public Header header;
	public String fileName, sourceName;
	public String creator = "CRUSH", creatorVersion = CRUSH.version; 
	
	public boolean verbose = false;
	
	// 2 pi sigma^2 = a^2
	// a = sqrt(2 pi) sigma
	//   = sqrt(2 pi) fwhm / 2.35
	public static double fwhm2size = Math.sqrt(2.0 * Math.PI) / Util.sigmasInFWHM;
	
	public Image() {}
	
	public Image(int sizeX, int sizeY) {
		setSize(sizeX, sizeY);
	}
	
	public Image(double[][] data) {
		this.data = data;
		noFlag();
	}
	
	public Image(double[][] data, int[][] flag) {
		this.data = data;
		this.flag = flag;
	}
	
	@Override
	public Object clone() {
		try {return super.clone(); }
		catch(CloneNotSupportedException e) { return null; }
	}
	
	public Image copy() {
		Image copy = (Image) clone();

		// Make a copy of the fundamental data
		copy.data = (double[][]) copyOf(data);
		copy.beam = (double[][]) copyOf(beam);
		copy.flag = (int[][]) copyOf(flag);
		
		copy.dataUnit = (Unit) dataUnit.clone();
		copy.dataUnitPerBeam = (Unit) dataUnitPerBeam.clone();
		
		return copy;
	}

	
	public final int sizeX() { return data.length; }
	
	public final int sizeY() { return data[0].length; }
	
	public void setSize(int x, int y) {
		beam = new double[][] {{1.0}};
		data = new double[x][y];
		flag = new int[x][y];
		for(int i=0; i<x; i++) Arrays.fill(flag[i], 1); 
	}
	
	public void noFlag() {
		flag = new int[sizeX()][sizeY()];		
	}
	
	public boolean isHorizontal() {
		return projection.getReference() instanceof HorizontalCoordinates;
	}
	
	public boolean isEquatorial() {
		return projection.getReference() instanceof EquatorialCoordinates;
	}
	
	public boolean isEcliptic() {
		return projection.getReference() instanceof EclipticCoordinates;
	}
	
	public boolean isGalactic() {
		return projection.getReference() instanceof GalacticCoordinates;
	}
	
	public boolean isSuperGalactic() {
		return projection.getReference() instanceof SuperGalacticCoordinates;
	}
	
	public double getImageFWHM() {
		return Math.hypot(instrument.resolution, smoothFWHM);
	}
	
	protected void calcBeamUnits() {
		final double imageFWHM = getImageFWHM();
		dataUnitPerBeam.value = dataUnit.value * instrument.resolution*instrument.resolution / (imageFWHM * imageFWHM);	
	}
	
	public double getImageBeamArea() {
		final double A = fwhm2size * getImageFWHM();
		return A*A;		
	}
	
	public double getPointsPerSmoothingBeam() {
		final double nx = Math.max(1.0, fwhm2size * smoothFWHM / delta);
		return nx * nx;
	}
	
	public void setResolution(double value) { 
		delta = value; 
		smoothFWHM = Math.max(smoothFWHM, delta / fwhm2size);
		calcBeamUnits();
	}
	
	public double getResolution() {
		return delta;
	}

	public void setDataUnit(Unit unit) {
		dataUnit = unit;
		dataUnitPerBeam.name = dataUnit.name + "/beam";
		calcBeamUnits();
	}
	

	public double getExtFilterCorrectionFactor(double FWHM) {
		double effectiveFilterFWHM2 = FWHM*FWHM + extFilterFWHM*extFilterFWHM;
		double effectiveFWHM2 = FWHM*FWHM + smoothFWHM*smoothFWHM;
		if(Double.isNaN(extFilterFWHM)) return 1.0;
		return 1.0 / (1.0 - effectiveFWHM2/effectiveFilterFWHM2);
	}
	
	public void setImage(double[][] image) { for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) data[i][j] = image[i][j]; }

	public void addImage(double[][] image, double scale) {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) data[i][j] += scale * image[i][j];
	}

	public void addImage(double[][] image) { addImage(image, 1.0); }
	
	public void clear() {
		for(int i=0; i<sizeX(); i++) {
			Arrays.fill(data[i], 0.0);
			Arrays.fill(flag[i], 1);
		}
	}

	public void reset() {
		clear();
		smoothFWHM = delta / fwhm2size;
		extFilterFWHM = Double.NaN;
		beam = new double[][] {{ 1.0 }};
	}
	
	
	public final double valueAt(final MapIndex index) {
		return valueAtIndex(index.i, index.j);
	}
	
	public double valueAtIndex(int i, int j) { return data[i][j]; }
	
	public double valueAtIndex(Vector2D index) {
		return valueAtIndex(index.x, index.y);
	}
	
	// Interpolate (linear at edges, quadratic otherwise)	
	public double valueAtIndex(double ic, double jc) {		
		final int i = (int) Math.round(ic);
		final int j = (int) Math.round(jc);
		
		final int sizeX = sizeX();
		final int sizeY = sizeY();
		
		if(ic < 0.0 || ic >= sizeX || jc < 0.0 || jc >= sizeY) return Double.NaN;
		
		final double y0 = data[i][j];
		double a=0.0, b=0.0, c=0.0, d=0.0;
		
		if(i > 0 && i < sizeX-1) if((flag[i+1][j] | flag[i-1][j]) == 0) {
			a = 0.5 * (data[i+1][j] + data[i-1][j]) - y0;
			c = 0.5 * (data[i+1][j] - data[i-1][j]);
		}
		if(j > 0 && j < sizeY-1) if((flag[i][j+1] | flag[i][j-1]) == 0) {
			b = 0.5 * (data[i][j+1] + data[i][j-1]) - y0;
			d = 0.5 * (data[i][j+1] - data[i][j-1]);
		}
		
		ic -= i;
		jc -= j;
		
		return (a*ic+c)*ic + (b*jc+d)*jc + y0;
	}
	
	public double getWeight(int i, int j) { return 1.0; }
	
	public void scale(double value) {
		if(value == 1.0) return;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) data[i][j] *= value;
	}
	
	public void setUnit(Unit u) {
		unit = u;
	}
	
	public void setProjection(SphericalProjection projection) {
		this.projection = projection;
	}
	
	public double getMin() { 
		double min=Double.POSITIVE_INFINITY;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<data[i].length; j++) if(flag[i][j]==0) 
				if(data[i][j] < min) min = data[i][j];
		return min;
	}

	
	public double getMax() {
		double max=Double.NEGATIVE_INFINITY;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<data[i].length; j++) if(flag[i][j]==0) 
				if(data[i][j] > max) max = data[i][j];
		return max;
	}
	
	public Range getRange() {
		final Range range = new Range();
		range.empty();
		for(int i=0; i<sizeX(); i++) for(int j=0; j<data[i].length; j++) if(flag[i][j]==0) 
			range.include(data[i][j]);
		return range;
	}
	
	
	public MapIndex indexOfMax() {
		MapIndex index = new MapIndex();

		double peak = 0.0;

		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0)
			if(data[i][j] > peak) {
				peak = data[i][j];
				index.i = i;
				index.j = j;
			}

		return index;
	}
	
	public double mean() {
		double sum = 0.0;
		int n=0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			sum += data[i][j];
			n++;
		}
		return sum / n;
	}
	
	public double median() {
		double[] point = new double[countPoints()];

		if(point.length == 0) return 0.0;

		for(int i=0, k=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) point[k++] = data[i][j];

		Arrays.sort(point);
		int n = point.length;

		return n % 2 == 0 ? (point[n/2-1] + point[n/2]) / 2.0 : point[(n - 1) / 2]; 
	}
	
	public double getRMS() {
		double sum = 0.0;
		int n=0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			sum += data[i][j] * data[i][j];
			n++;
		}
		return n < 0 ? Double.NaN : Math.sqrt(sum / n);
	}
	
	public double getRobustRMS() {
		double[] chi2 = new double[countPoints()];
		if(chi2.length == 0) return 0.0;

		for(int i=0, k=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			final double value = data[i][j] * data[i][j];
			chi2[k++] = value * value;
		}
		// median(x^2) = 0.454937 * sigma^2 
		return Math.sqrt(ArrayUtil.median(chi2) / 0.454937);	
	}
	
	public void fluxClip(double level) {
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0)
			if(data[i][j] < level) flag[i][j] = 1;
	}
	
	public int countPoints() {
		int points = 0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) points++;
		return points;
	}
	
	public void crop(double dXmin, double dYmin, double dXmax, double dYmax) {		
		if(dXmin > dXmax) { double temp = dXmin; dXmin = dXmax; dXmax=temp; }
		if(dYmin > dYmax) { double temp = dYmin; dYmin = dYmax; dYmax=temp; }
		
		int imin = Math.max(0, (int)Math.floor(fracIndexOfdX(dXmin)));
		int jmin = Math.max(0, (int)Math.floor(fracIndexOfdY(dYmin)));
		int imax = Math.min(sizeX()-1, (int)Math.ceil(fracIndexOfdX(dXmax)));
		int jmax = Math.min(sizeY()-1, (int)Math.ceil(fracIndexOfdY(dYmax)));
	
		crop(imin, jmin, imax, jmax);
	}

	protected void crop(int imin, int jmin, int imax, int jmax) {
		if(verbose) System.err.println("Cropping.");

		double[][] olddata = data;
		int[][] oldflag = flag;
		
		setSize(imax-imin+1, jmax-jmin+1);
		
		for(int i=imin, i1=0; i<=imax; i++, i1++) for(int j=jmin, j1=0; j<=jmax; j++, j1++) {
			data[i1][j1] = olddata[i][j];
			flag[i1][j1] = oldflag[i][j];
		}

		refIndex.x -= imin;
		refIndex.y -= jmin;
	}
	
	public int[] getHorizontalIndexRange() {
		int min = sizeX()-1, max = 0;
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			if(i < min) min = i;
			if(i > max) max = i;
			break;
		}
		return max > min ? new int[] { min, max } : null;
	}
	
	public int[] getVerticalIndexRange() {
		int min = sizeY()-1, max = 0;
		for(int j=0; j<sizeY(); j++) for(int i=0; i<sizeX(); i++) if(flag[i][j] == 0) {
			if(j < min) min = j;
			if(j > max) max = j;
			break;
		}
		return max > min ? new int[] { min, max } : null;
	}
	
	public void autoCrop() {
		if(verbose) System.err.print("Auto-cropping. ");
		int[] hRange = getHorizontalIndexRange();
		int[] vRange = getVerticalIndexRange();
		if(verbose) System.err.println((hRange[1] - hRange[0] + 1) + "x" + (vRange[1] - vRange[0] + 1));
		this.crop(hRange[0], vRange[0], hRange[1], vRange[1]);
	}
	
	public boolean[][] getBooleanFlag() {
		boolean[][] mask = new boolean[sizeX()][sizeY()];
		for(int j=0; j<sizeY(); j++) for(int i=0; i<sizeX(); i++) mask[i][j] = flag[i][j] > 0;
		return mask;
	}
	
	public void flag(boolean[][] mask, int pattern) {
		for(int j=0; j<sizeY(); j++) for(int i=0; i<sizeX(); i++) if(mask[i][j]) flag[i][j] |= pattern;	 
	}
	
	public void unflag(boolean[][] mask, int pattern) {
		pattern = ~pattern;
		for(int j=0; j<sizeY(); j++) for(int i=0; i<sizeX(); i++) if(mask[i][j]) flag[i][j] &= pattern;	 
	}
	
	public double level(boolean robust) {
		double level = robust ? median() : mean();
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) data[i][j] -= level;
		return level;
	}
	
	
	public double countBeams() { return getArea() / getImageBeamArea(); }

	public double getArea() { return countPoints() * delta * delta; }

	
	// Convolves image to the specified beam resolution
	// by a properly chosen convolving beam...
	public void convolveTo(double FWHM) {
		if(smoothFWHM >= FWHM) return;
		convolve(Math.sqrt(FWHM * FWHM - smoothFWHM * smoothFWHM));
	}	

	public double[][] getConvolvedTo(double[][] image, double FWHM) {
		if(smoothFWHM >= FWHM) return image; 
		return getConvolved(image, Math.sqrt((FWHM * FWHM - smoothFWHM * smoothFWHM)));
	}

	public void convolve(double FWHM) {
		int size = (int)Math.ceil(2.0 * FWHM/delta);
		int step = (int)Math.ceil(FWHM/(5.0 * delta));
		size = 2*size + 1;	
		fastConvolve(getGaussian(FWHM, size), step);
		smoothFWHM = Math.hypot(smoothFWHM, FWHM);
		calcBeamUnits();
	}

	public double[][] getConvolved(double[][] image, double FWHM) {
		int size = (int)Math.ceil(2.0*FWHM/delta);
		int step = (int)Math.ceil(FWHM/(5.0 * delta));
		size = 2*size + 1;
		return getFastConvolved(image, getGaussian(FWHM, size), null, step);
	}   

	public void convolve(double[][] beam) {
		double[][] beamw = new double[sizeX()][sizeY()];
		data = getConvolved(data, beam, beamw);
		this.beam = beam;
	}
	
	public void fastConvolve(double[][] beam, int step) {
		double[][] beamw = new double[sizeX()][sizeY()];
		data = getFastConvolved(data, beam, beamw, step);
		this.beam = beam;
	}
	
	public double[][] getConvolved(double[][] image, double[][] beam, double[][] beamw) {
		double[][] convolved = new double[sizeX()][sizeY()];
		final int ic = (beam.length-1) / 2;
		final int jc = (beam[0].length-1) / 2;
		
		final WeightedPoint result = new WeightedPoint();
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			getConvolved(image, i, j, beam, ic, jc, result);
			convolved[i][j] = result.value;
			if(beamw != null) beamw[i][j] = result.weight;
		}
		
		return convolved;
	}
	
	
	// Do the convolution proper at the specified intervals (step) only, and interpolate (quadratic) inbetween
	public double[][] getFastConvolved(double[][] image, double[][] beam, double[][] beamw, int step) {
		if(step < 2) return getConvolved(image, beam, beamw);
		
		final int ic = (beam.length-1) / 2;
		final int jc = (beam[0].length-1) / 2;
		
		final WeightedPoint result = new WeightedPoint();
		
		final int nx = sizeX()/step + 1;
		final int ny = sizeY()/step + 1;
 
		final Image signalImage = new Image(nx, ny);
		Image weightImage = null;
		
		if(beamw != null) {
			weightImage = (Image) signalImage.clone();
			weightImage.data = new double[nx][ny];
		}
			
		for(int i=0, i1=0; i<sizeX(); i+=step, i1++) for(int j=0, j1=0; j<sizeY(); j+=step, j1++) {
			getConvolved(image, i, j, beam, ic, jc, result);
			signalImage.data[i1][j1] = result.value;
			if(beamw != null) weightImage.data[i1][j1] = result.weight;
			signalImage.flag[i1][j1] = result.weight == 0.0 ? 1 : 0;
		}
		
		final double[][] convolved = new double[sizeX()][sizeY()];
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			final double i1 = (double) i / step;
			final double j1 = (double) j / step;
			convolved[i][j] = signalImage.valueAtIndex(i1, j1);
			if(beamw != null) beamw[i][j] = weightImage.valueAtIndex(i1, j1);
		}
		
		return convolved;
	}
	
	
	// Beam fitting: I' = C * sum(wBI) / sum(wB2)
	// I(x) = I -> I' = I -> C = sum(wB2) / sum(wB)
	// I' = sum(wBI)/sum(wB)
	// rms = Math.sqrt(1 / sum(wB))
	public void getConvolved(final double[][] image, final int i, final int j, final double[][] beam, int ic, int jc, WeightedPoint result) {
		result.noData();
		
		ic += i;
		jc += j;
		
		// ib = i + ic - i1
		final int fromib = Math.min(beam.length-1, ic);
		final int fromjb = Math.min(beam[0].length-1, jc);
		final int toib = Math.max(0, ic - sizeX());
		final int tojb = Math.max(0, jc - sizeY());
			
		for(int ib=fromib, i1=ic-fromib; ib>toib; ib--,i1++) for(int jb=fromjb, j1=jc-fromjb; jb>tojb; jb--,j1++) if(flag[i1][j1] == 0) {
			final double B = beam[ib][jb];
			if(B != 0.0) {
				result.value += getWeight(i1, j1) * B * image[i1][j1];
				result.weight += getWeight(i1, j1) * Math.abs(B);		    
			}
		}

		result.value /= result.weight;
	}
	
	
	
	public void filterAbove(double FWHM) { filterAbove(FWHM, flag); }

	public void filterAbove(double FWHM, int[][] skip) {
		Image extended = copy();
		extended.flag = skip;
		extended.convolveTo(FWHM);
	
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) data[i][j] -= extended.data[i][j];		

		if(Double.isNaN(extFilterFWHM)) extFilterFWHM = FWHM;
		else extFilterFWHM = 1.0/Math.sqrt(1.0/(extFilterFWHM * extFilterFWHM) + 1.0/(FWHM*FWHM));

		filterCorrect(instrument.resolution, skip);
	}
	
	public void filterCorrect(double FWHM, int[][] skip) {
		if(!Double.isNaN(correctingFWHM)) return;
		
		final double filterC = getExtFilterCorrectionFactor(FWHM);
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(skip[i][j] != 0) scale(i, j, filterC);
		
		correctingFWHM = FWHM;
	}
	
	public void scale(int i, int j, double factor) {
		data[i][j] *= factor;
	}

	public double[][] getBeam(double FWHM) {
		int size = 3 * (int)Math.ceil(FWHM/delta) + 1;
		return getGaussian(FWHM, size);
	}
	
	public double[][] getGaussian(double FWHM, int size) {
		final double[][] beam = new double[size][size];
		final double sigma = FWHM / Util.sigmasInFWHM;
		final double A = -0.5 * (delta * delta) / (sigma * sigma);
		final double center = (size-1) / 2.0;

		for(int i=0; i<size; i++) for(int j=0; j<size; j++) {
			double dx = i - center;
			double dy = j - center;

			beam[i][j] = Math.exp(A*(dx*dx+dy*dy));
		}
		return beam;
	}	
	
	public void addGaussianAt(double flux, double dx, double dy, double FWHM, double beamRadius, float[][] renorm) {
		FWHM /= delta;
		
		final double sigma = FWHM / Util.sigmasInFWHM;
		final double A = -0.5 / (sigma*sigma);
		final double i0 = indexOfdX(dx);
		final double j0 = indexOfdY(dy);
		final int di = (int)Math.ceil(beamRadius * FWHM);
		
		final int imin = Math.max(0, (int)Math.floor(i0-di));
		final int imax = Math.min(sizeX()-1, (int)Math.ceil(i0+di));
		final int jmin = Math.max(0, (int)Math.floor(j0-di));
		final int jmax = Math.min(sizeY()-1, (int)Math.ceil(j0+di));
		
		for(int i=imin; i<=imax; i++) {
			double di2 = i - i0;
			di2 *= di2;
			
			for(int j=jmin; j<=jmax; j++) {
				double dj2 = j - j0;
				dj2 *= dj2;
				
				final double a = Math.exp(A * (di2 + dj2));
				
				data[i][j] += a * flux;
				renorm[i][j] += (float) a;
			}
		}
	}
	
	public void addGaussianAt(Image from, int i, int j, double FWHM, double beamRadius, float[][] renorm) {
		addGaussianAt(from.data[i][j], from.dXofIndex(i), from.dYofIndex(j), FWHM, 1.0, renorm);
	}
	
	public void regridTo(final SkyMap regrid) {		
		// replace pixels with Gaussians whose FWHM is the diagonal of the pixel (i.e. enclosing)...
		double FWHM = Math.sqrt(2.0 * delta * delta);
		
		// Remove the effect of pixellization from the smoothing to calculate the explicit smoothing...
		smoothFWHM = smoothFWHM * smoothFWHM - delta * delta / (fwhm2size * fwhm2size);
		if(smoothFWHM < 0.0) smoothFWHM = 0.0;
		smoothFWHM = Math.sqrt(smoothFWHM);
		
		// Add in the smoothing of the regrid algorithm...
		smoothFWHM = Math.hypot(smoothFWHM, FWHM);
		
		if(regrid.delta > delta) smoothFWHM = Math.hypot(smoothFWHM, regrid.delta / fwhm2size);
		
		regrid.clear();
		float renorm[][] = new float[regrid.sizeX()][regrid.sizeY()];
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j]==0) 
			regrid.addGaussianAt(this, i, j, FWHM, 1.0, renorm);
			
		regrid.normalize();
		
		// renormalize the weights with the taper...
		for(int i=0; i<regrid.sizeX(); i++) {
			Arrays.fill(regrid.flag[i], 1);
			for(int j=0; j<regrid.sizeY(); j++) regrid.weight[i][j] /= renorm[i][j];
		}
		
		// recalculate the flags...	
		for(int i=0; i<sizeX(); i++) {
			final int mini = Math.max(0, regrid.indexOfdX(dXofIndex(i-0.5)));
			final int maxi = Math.min(regrid.sizeX()-1, regrid.indexOfdX(dXofIndex(i+0.5)));
			
			for(int j=0; j<sizeY(); j++) if(flag[i][j]==0) {	
				final int minj = Math.max(0, regrid.indexOfdY(dYofIndex(j-0.5)));
				final int maxj = Math.min(regrid.sizeY()-1, regrid.indexOfdY(dYofIndex(j+0.5)));
			
				for(int i1=mini; i1<=maxi; i1++) for(int j1=minj; j1<=maxj; j1++) regrid.flag[i1][j1] = 0;
			}
			
		}
		
		regrid.sanitize();
	}
	
	
	public Vector2D[] getHistogram(double[][] image, double binSize) {
		Range range = getRange();
		
		int bins = 1 + (int)Math.round(range.max / binSize) - (int)Math.round(range.min / binSize);
		Vector2D[] bin = new Vector2D[bins];
		for(int i=0; i<bins; i++) bin[i] = new Vector2D(i*binSize, 0.0);
		
		for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) if(flag[i][j] == 0) {
			bin[(int)Math.round(image[i][j] / binSize)].y++;
		}
		
		return bin;
	}
	
	
	public void read(BasicHDU HDU) throws FitsException {		
		Object image = HDU.getData().getData();
			
		try {
			final float[][] fdata = (float[][]) image;
			for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) {
				if(!Float.isNaN(fdata[j][i])) data[i][j] = fdata[j][i] * unit.value;	    
				else flag[i][j] |= 1;
			}
		}
		catch(ClassCastException e) {
			final double[][] ddata = (double[][]) image;
			for(int i=0; i<sizeX(); i++) for(int j=0; j<sizeY(); j++) {
				if(!Double.isNaN(ddata[j][i])) data[i][j] = ddata[j][i] * unit.value;	    
				else flag[i][j] |= 1;
			}				
		}
	}
	
	public ImageHDU makeHDU(double unitValue, String unitName, String name) throws HeaderCardException, FitsException {
		float[][] fitsImage = new float[sizeY()][sizeX()];
		for(int i=0; i<sizeX(); i++) for(int j=0; j<data[0].length; j++) {
			if(flag[i][j] == 0.0) fitsImage[j][i] = (float) (data[i][j] / unitValue);
			else fitsImage[j][i] = Float.NaN;
		}
		ImageHDU hdu = (ImageHDU)Fits.makeHDU(fitsImage);

		hdu.addValue("EXTNAME", name, "The type of data contained in this HDU");
		addCoordinateInfo(hdu);

		Range range = getRange();
		
		hdu.addValue("DATAMIN", range.min / unitValue, "");
		hdu.addValue("DATAMAX", range.max / unitValue, "");

		hdu.addValue("BZERO", 0.0, "Zeroing level of the image data");
		hdu.addValue("BSCALE", 1.0, "Scaling of the image data");
		hdu.addValue("BUNIT", unitName, "The image data unit.");

		return hdu;
	}

	public void addCoordinateInfo(BasicHDU hdu) throws HeaderCardException {
		SphericalCoordinates reference = projection.getReference();
		
		hdu.addValue("CTYPE1", reference.coordinateSystem.get(0).wcsName + "-" + projection.getFitsID(), "Longitude Coordinate Type and Projection");
		hdu.addValue("CRPIX1", refIndex.x + 1, "Reference grid position");
		hdu.addValue("CRVAL1", reference.longitude()/Unit.deg, "Longitude at reference position (deg)");
		hdu.addValue("CDELT1", (reference.isReverseLongitude() ? -1 : 1) * delta/Unit.deg, "Grid spacing (deg)");	
		hdu.addValue("CROTA1", 0.0, "Axis rotation (deg)");
		
		hdu.addValue("CTYPE2", reference.coordinateSystem.get(1).wcsName + "-" + projection.getFitsID(), "Latitude Coordinate Type and Projection");
		hdu.addValue("CRPIX2", refIndex.y + 1, "Reference grid position");
		hdu.addValue("CRVAL2", reference.latitude()/Unit.deg, "Latitude at reference position (deg)");
		hdu.addValue("CDELT2", (reference.isReverseLatitude() ? -1 : 1) * delta/Unit.deg, "Grid spacing (deg)");
		hdu.addValue("CROTA2", 0.0, "Axis rotation (deg)");

		if(reference instanceof EquatorialCoordinates) {
			EquatorialCoordinates equatorial = (EquatorialCoordinates) reference;
			hdu.addValue("EQUINOX", equatorial.epoch.getYear(), "The catalog coordinate epoch.");
			hdu.addValue("RADESYS", equatorial.epoch instanceof JulianEpoch ? "FK5" : "FK4", "The Coordinate System Convention.");
		}	
		
		if(reference instanceof EclipticCoordinates) {
			EclipticCoordinates ecliptic = (EclipticCoordinates) reference;
			hdu.addValue("EQUINOX", ecliptic.epoch.getYear(), "The catalog coordinate epoch.");
			hdu.addValue("RADESYS", ecliptic.epoch instanceof JulianEpoch ? "FK5" : "FK4", "The Coordinate System Convention.");			
		}	
	}
	
	public void addLongHierarchKey(Cursor cursor, String key, String value) throws FitsException, HeaderCardException {
		addLongHierarchKey(cursor, key, 0, value);
	}

	public void addLongHierarchKey(Cursor cursor, String key, int part, String value) throws FitsException, HeaderCardException {
		if(value.length() == 0) value = "true";
		
		String alt = part > 0 ? "." + part : "";

		int available = 69 - (key.length() + alt.length() + 3);

		if(value.length() < available) cursor.add(new HeaderCard("HIERARCH." + key + alt, value, null));
		else { 
			if(alt.length() == 0) {
				part = 1;
				alt = "." + part;
				available -= 2;
			}

			cursor.add(new HeaderCard("HIERARCH." + key + alt, value.substring(0, available), null));
			addLongHierarchKey(cursor, key, (char)(part+1), value.substring(available)); 
		}
	}
	
	
	public void getCoordinateInfo(Header header) throws HeaderCardException {
		String type = header.getStringValue("CTYPE1");
	
		int sizeX = header.getIntValue("NAXIS1");
		int sizeY = header.getIntValue("NAXIS2");
		
		setSize(sizeX, sizeY);
		
		try { projection = SphericalProjection.forName(type.substring(5, 8)); }
		catch(Exception e) {
			System.err.println("ERROR! Unknown projection " + type.substring(5, 8));
		}
		
		SphericalCoordinates reference = null;
		
		if(type.toUpperCase().startsWith("ALON")) {
			reference = new HorizontalCoordinates();
			reference.parse(header);
		}
		else if(type.toUpperCase().startsWith("RA--")) {
			reference = new EquatorialCoordinates();
			reference.parse(header);
		}	
		else System.err.println("ERROR! Unknown Coordinate Definition " + type);
		
		projection.setReference(reference);
		
		refIndex.x = header.getDoubleValue("CRPIX1") - 1;
		refIndex.y = header.getDoubleValue("CRPIX2") - 1;
		
		delta = header.getDoubleValue("CDELT2") * Unit.deg;
	
		// Adjust the index to the determined pointing
		double di0 = (isHorizontal() ? -1 : 1) * header.getDoubleValue("RAP") * Unit.arcsec / delta;
		double dj0 = header.getDoubleValue("DECP") * Unit.arcsec / delta;
	
		refIndex.x += di0;
		refIndex.y += dj0;
	
		if(di0 != 0.0 || dj0 != 0.0)
			System.out.println("Shifting -> (" + 
					header.getDoubleValue("RAP") * Unit.arcsec / delta +  "," +
					header.getDoubleValue("DECP") * Unit.arcsec / delta + ")" );
		
	}
	
	@Override
	public String toString() {	
		SphericalCoordinates reference = projection.getReference();
		String projectionName = reference.getClass().getSimpleName();
		projectionName = projectionName.substring(0, projectionName.length() - "Coordinates".length());
		
		String info =
			"  " + projectionName + ": " + reference.toString() + "\n" +
			"  Projection: " + projection.getFullName() + " (" + projection.getFitsID() + ")\n" + 
			"  Map Size: " + sizeX() + " x " + sizeY() + " pixels. (" 
			+ Util.f1.format(sizeX() * delta / Unit.arcmin) + " x " + Util.f1.format(sizeY() * delta / Unit.arcmin) + " arcmin)." + "\n" +
			"  Grid Spacing: " + Util.f2.format(delta / Unit.arcsec) + " arcsec.\n";
		
		return info;
	}
	
	public void info() {
		System.err.println("\n" + toString() + "\n");
	}
	
	
	public void write(String name, Unit unit) throws HeaderCardException, FitsException, IOException {
		Fits fits = new Fits();	

		fits.addHDU(makeHDU(unit.value, unit.name, "Image"));
		editHeader(fits);

		BufferedDataOutputStream file = new BufferedDataOutputStream(new FileOutputStream(name));

		fits.write(file);	
		System.out.println(" Written Image to " + name);
	}
	
	public void editHeader(Fits fits) throws HeaderCardException, FitsException, IOException {
		
		nom.tam.util.Cursor cursor = fits.getHDU(0).getHeader().iterator();

		// Go to the end of the header cards...
		while(cursor.hasNext()) cursor.next();

		double imageFWHM = getImageFWHM();
		
		cursor.add(new HeaderCard("CREATOR", creator, "The software that created the image."));
		cursor.add(new HeaderCard("CRUSHVER", creatorVersion, "Version information."));
		cursor.add(new HeaderCard("DATE", FitsDate.getFitsDateString(), "The date of creation"));
		cursor.add(new HeaderCard("OBJECT", sourceName, "Source name as it appear in the raw data."));	
		//cursor.add(new HeaderCard("ORIGIN", "Caltech", "California Institute of Technology"));
		cursor.add(new HeaderCard("TELESCOP", instrument.getTelescopeName(), "Telescope name."));
		cursor.add(new HeaderCard("INSTRUME", instrument.name, "The instrument used."));	
		
		// The data descriptors
		cursor.add(new HeaderCard("BEAM", instrument.resolution / Unit.arcsec, "The instrument FWHM (arcsec) of the beam."));
		cursor.add(new HeaderCard("BMAJ", imageFWHM, "The beam major axis (radians)"));
		cursor.add(new HeaderCard("BMIN", imageFWHM, "The beam minor axis (radians)."));
		cursor.add(new HeaderCard("BPA", 0.0, "The beam position angle (radians)."));

		cursor.add(new HeaderCard("SMOOTH", smoothFWHM / Unit.arcsec, "The FWHM (arcsec) of the smoothing applied."));	
		cursor.add(new HeaderCard("SMTHRMS", true, "Is the Noise (RMS) image smoothed?"));
		cursor.add(new HeaderCard("RESOLUTN", imageFWHM / Unit.arcsec, "The effective image FWHM (arcsec)."));	
		
		
		if(!Double.isNaN(extFilterFWHM)) 
			cursor.add(new HeaderCard("EXTFLTR", extFilterFWHM / Unit.arcsec, "Large-scale structure filtering FWHM."));
			
		// Add the system descriptors...
		cursor.add(new HeaderCard("JAVA", Util.getProperty("java.vendor"), "Java vendor name."));
		cursor.add(new HeaderCard("JAVAVER", Util.getProperty("java.version"), "The Java version."));
		cursor.add(new HeaderCard("JAVAHOME", Util.getProperty("java.home"), "Java location."));
		cursor.add(new HeaderCard("JRE", Util.getProperty("java.runtime.name"), "Java Runtime Environment."));
		cursor.add(new HeaderCard("JREVER", Util.getProperty("java.runtime.version"), "JRE version."));
		cursor.add(new HeaderCard("JVM", Util.getProperty("java.vm.name"), "Java Virtual Machine."));
		cursor.add(new HeaderCard("JVMVER", Util.getProperty("java.vm.version"), "JVM version."));
		
		cursor.add(new HeaderCard("OS", Util.getProperty("os.name"), "Operation System name."));
		cursor.add(new HeaderCard("OSVER", Util.getProperty("os.version"), "OS version."));
		cursor.add(new HeaderCard("OSARCH", Util.getProperty("os.arch"), "OS architecture."));
		
		cursor.add(new HeaderCard("CPUS", Runtime.getRuntime().availableProcessors(), "Number of CPU cores/threads available."));
		cursor.add(new HeaderCard("DMODEL", Util.getProperty("sun.arch.data.model"), "Bits in data model."));
		cursor.add(new HeaderCard("CPUEND", Util.getProperty("sun.cpu.endian"), "CPU Endianness."));
		cursor.add(new HeaderCard("MAXMEM", Runtime.getRuntime().maxMemory() / (1024 * 1024), "MB of available memory."));
			
		cursor.add(new HeaderCard("COUNTRY", Util.getProperty("user.country"), "The user country."));
		cursor.add(new HeaderCard("LANGUAGE", Util.getProperty("user.language"), "The user language."));
		
		if(instrument != null) if(instrument.options != null) instrument.options.editHeader(cursor);
	}
	

	public final double dXofIndex(final double i) { return (i-refIndex.x) * delta; }
	   
    public final double dYofIndex(final double j) { return (j-refIndex.y) * delta; }
    
    public final int indexOfdX(final double dX) { return (int) Math.round(fracIndexOfdX(dX)); }
    
    public final int indexOfdY(final double dY) { return (int) Math.round(fracIndexOfdY(dY)); }

    public final double fracIndexOfdX(final double dX) { return refIndex.x + dX/delta; }
    
    public final double fracIndexOfdY(final double dY) { return refIndex.y + dY/delta; }
	
    
    public final void getIndex(final Vector2D offset, final MapIndex index) {
    	index.i = indexOfdX(offset.x);
    	index.j = indexOfdY(offset.y);    	
    }
    
    public final void getIndex(final Vector2D offset, final Vector2D index) {
    	index.x = indexOfdX(offset.x);
    	index.y = indexOfdY(offset.y);    	
    }
    
    public final void getOffset(final MapIndex index, final Vector2D offset) {
    	offset.x = dXofIndex(index.i);
    	offset.y = dYofIndex(index.j);
    }
    
    public final void getOffset(final Vector2D index, final Vector2D offset) {
    	offset.x = dXofIndex(index.x);
    	offset.y = dYofIndex(index.y);
    }
    
    public boolean containsIndex(double i, double j) {
    	return i >= 0 && i < sizeX() && i >= 0 && j < sizeY();
    }
    	
    public static Object copyOf(final Object image) {

		if(image instanceof double[][]) {
			final double[][] orig = (double[][]) image;
			final double[][] copy = new double[orig.length][orig[0].length];
			for(int i=0; i<orig.length; i++) System.arraycopy(orig[i], 0, copy[i], 0, orig[0].length);	
			return copy;
		}	
		else if(image instanceof int[][]) {
			final int[][] orig = (int[][]) image;
			final int[][] copy = new int[orig.length][orig[0].length];
			for(int i=0; i<orig.length; i++) System.arraycopy(orig[i], 0, copy[i], 0, orig[0].length);	
			return copy;
		}

		return null;
	}
    
    
    
}

