/*******************************************************************************
 * 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.text.*;
import java.util.*;

import util.*;


public class GaussianSource extends CircularRegion {
	public DataPoint peak, FWHM;
	public boolean isCorrected = false;
	
	public GaussianSource() { }

	public GaussianSource(String line, AstroImage forImage) throws ParseException {
		super(line, forImage);
	}
	
	public GaussianSource(AstroMap map, Vector2D offset, double r) {
		super(map, offset, r);
	}
	
	public GaussianSource(SphericalCoordinates coords, double r) {
		super(coords, r);
	}
	
	public Bounds getBounds(AstroImage image, double beams) {
		double origRadius = radius.value;
		radius.value = beams * image.getImageFWHM();
		Bounds bounds = getBounds(image);
		radius.value = origRadius;
		return bounds;
	}
	
	@Override
	public DataPoint finetunePeak(AstroMap map) {
		AstroImage rms = map.getRMSImage();
		peak = super.finetunePeak(map);
		Vector2D centerIndex = getIndex(map.grid);
		if(peak == null) {
			peak.value = map.valueAtIndex(centerIndex);
			peak.setRMS(rms.valueAtIndex(centerIndex));
		}
		return peak;
	}
	
	public void scale(double factor) {
		peak.scale(factor);		
	}
	
	public void getChi2(AstroMap map, WeightedPoint chi2, double level) {
		chi2.noData();
		Bounds bounds = getBounds(map);
		Vector2D centerIndex = getIndex(map.grid);
		final double sigma = radius.value / Util.sigmasInFWHM / map.getResolution();
		final double A = -0.5 / (sigma*sigma);

		for(int i=bounds.fromi; i<=bounds.toi; i++) for(int j=bounds.fromj; j<=bounds.toj; j++) if(map.flag[i][j] == 0) {
			final double di = i-centerIndex.x;
			final double dj = j-centerIndex.y;
			final double dev = (map.data[i][j] - level - peak.value * Math.exp(A*(di*di+dj*dj))) / map.getRMS(i, j);
			chi2.value += dev * dev;
			chi2.weight += 1.0;
		}
	}
	
	public void addGaussian(AstroImage image, double FWHM, double scaling) {
		Bounds bounds = getBounds(image, 3.0);
		Vector2D centerIndex = getIndex(image.grid);
		final double sigma = radius.value / Util.sigmasInFWHM / image.getResolution();
		final double A = -0.5 / (sigma*sigma);
			
		for(int i=bounds.fromi; i<=bounds.toi; i++) for(int j=bounds.fromj; j<=bounds.toj; j++) {
			final double di = i-centerIndex.x;
			final double dj = j-centerIndex.y;
			image.data[i][j] += scaling * peak.value * Math.exp(A*(di*di+dj*dj));
		}
		
	}
	
	public void add(AstroImage image) { add(image, null); }
	
	public void add(AstroImage image, Collection<Region> others) { add(image, radius.value, 1.0, others); }
	
	public void addPoint(AstroImage image) { addPoint(image, null); }
	
	public void addPoint(AstroImage image, Collection<Region> others) { add(image, image.getImageFWHM(), 1.0, others); }	
	
	public void subtract(AstroImage image) { subtract(image, null); }
	
	public void subtract(AstroImage image, Collection<Region> others) { add(image, radius.value, -1.0, others); }
	
	public void subtractPoint(AstroImage image) { subtractPoint(image, null); }
	
	public void subtractPoint(AstroImage image, Collection<Region> others) { add(image, image.getImageFWHM(), -1.0, others); }
	
	
	public void add(AstroImage image, final double FWHM, final double scaling) {
		add(image, FWHM, scaling, null);
	}
	
	public void add(AstroImage image, final double FWHM, final double scaling, final Collection<Region> others) {			
		double filterFWHM = Double.NaN;
		double imageFWHM = image.getImageFWHM();
	
		if(!Double.isNaN(image.extFilterFWHM)) 
			filterFWHM = Math.sqrt(
					image.extFilterFWHM * image.extFilterFWHM 
					+ FWHM * FWHM 
					- imageFWHM * imageFWHM
			);		
		
		// Remove source only after the filtering corrections are applied to its flux...
		addGaussian(image, FWHM, scaling);

		// Correct for filtering.
		// Consider that only the tip of the source might escape the filter...
		
		double filterBeamScale = Math.pow(FWHM / filterFWHM, 2.0);
		
		if(!Double.isNaN(image.extFilterFWHM)) if(image instanceof AstroMap) {
			AstroMap map = (AstroMap) image;
	
			filterBeamScale *= Double.isNaN(map.filterBlanking) ? 1.0 : Math.min(1.0, map.filterBlanking / peak.significance());
			addGaussian(image, filterFWHM, -scaling * filterBeamScale);
		}

		final double sigma = radius.value / Util.sigmasInFWHM / image.getResolution();
		final double A = -0.5 / (sigma*sigma);		

		Vector2D centerIndex = getIndex(image.grid);
		// Adjust prior detections.for the filtering around this one.
		if(others != null) for(Region region : others) if(region instanceof GaussianSource) if(region != this) {
			GaussianSource source = (GaussianSource) region;
			Vector2D sourceIndex = source.getIndex(image.grid);
			final double di = sourceIndex.x - centerIndex.x;
			final double dj = sourceIndex.y - centerIndex.y;
			
			if(!Double.isNaN(image.extFilterFWHM))
				source.peak.value -= scaling * filterBeamScale * peak.value * Math.exp(A*(di*di+dj*dj));
		}
		
	}
	
	public double getCorrectionFactor(AstroMap map, double FWHM) {	
		double correction = 1.0;	
		
		// Correct for filtering.
		// Consider that only the tip of the source might escape the filter...
		if(!Double.isNaN(map.extFilterFWHM)) {
			double filterFraction = Double.isNaN(map.filterBlanking) ? 1.0 : Math.min(1.0, map.filterBlanking / peak.significance());
			double filtering = 1.0 - 1.0 / map.getExtFilterCorrectionFactor(FWHM);;
			correction *= 1.0 / (1.0 - filtering * filterFraction);
		}
		
		return correction;
	}
	
	
	public void correct(AstroMap map, double FWHM) {	
		if(isCorrected) throw new IllegalStateException("Source is already corrected.");
		double correction = getCorrectionFactor(map, FWHM);
		peak.scale(correction);
		isCorrected = true;
	}
	

	public void uncorrect(AstroMap map, double FWHM) {
		if(!isCorrected) throw new IllegalStateException("Source is already uncorrected.");
		double correction = getCorrectionFactor(map, FWHM);
		peak.scale(1.0 / correction);
		isCorrected = false;
	}
	
	
	// Formula from Kovacs et al. (2006)
	public void setSearchRadius(AstroImage image, double pointingRMS) {
		double beamSigma = image.instrument.resolution / Util.sigmasInFWHM;
		radius.value = Math.sqrt(4.0 * pointingRMS * pointingRMS - 2.0 * beamSigma * beamSigma * Math.log(1.0 - 2.0 / peak.significance()));
	}

	
	@Override
	public void parse(String line, AstroImage forImage) throws ParseException {
		if(line == null) return;
		if(line.length() == 0) return;
		if(line.charAt(0) == '#' || line.charAt(0) == '!') return;

		StringTokenizer tokens = new StringTokenizer(line);

		if(tokens.countTokens() < 5) return;

		id = tokens.nextToken();

		EquatorialCoordinates coords = new EquatorialCoordinates();
		coords.parse(tokens.nextToken() + " " + tokens.nextToken() + " " + tokens.nextToken());
		//System.err.println(coords.toString());
		if(forImage.getReference() instanceof Precessing) {
			CoordinateEpoch epoch = ((Precessing) forImage.getReference()).getEpoch();
			coords.precess(epoch);
		}
			
		radius.value = Double.parseDouble(tokens.nextToken()) * Unit.arcsec;
		if(tokens.hasMoreTokens()) peak.value = Double.parseDouble(tokens.nextToken());
		
		String nextArg = tokens.hasMoreTokens() ? tokens.nextToken() : null;
		if(nextArg != null) {
			if(nextArg.equals("+-")) nextArg = tokens.nextToken();
			try { 
				peak.setRMS(Double.parseDouble(nextArg));
				nextArg = tokens.hasMoreTokens() ? tokens.nextToken() : null;
			}
			catch(NumberFormatException e) {}
		}
		
		Unit unit = nextArg == null ? new Unit("Jy/beam", forImage.getUnit("Jy/beam")) : new Unit(nextArg, forImage.getUnit(nextArg));
		while(tokens.hasMoreTokens()) comment += tokens.nextToken() + " ";
		comment.trim();
		
		peak.scale(unit.value);
	}

	@Override
	public String toCrushString(AstroImage image) {
		return id + "\t" + super.toCrushString(image) + "  " + getPeakString(image.unit);
	}
	
	@Override
	public String getComment() {
		return "s/n=" + Util.f2.format(peak.significance()) + " " + super.getComment();
	}
	
	@Override
	public StringTokenizer parseCrush(String line, AstroImage forImage) {
		StringTokenizer tokens = super.parseCrush(line, forImage);
		
		peak.value = Double.parseDouble(tokens.nextToken());
		String next = tokens.nextToken();
		if(next.equals("+-")) {
			peak.setRMS(Double.parseDouble(tokens.nextToken()));
			next = tokens.nextToken();
		}
		
		peak.value *= forImage.getUnit(next);
		
		return tokens;
	}
	
	public String getPeakString(Unit unit) {
		return DataPoint.toString(peak, unit);
	}
	
	public void centroid(AstroMap map) {
		Bounds bounds = getBounds(map, 2.0);
		Vector2D index = new Vector2D();
		double sumw = 0.0;
		
		for(int i=bounds.fromi; i<=bounds.toi; i++) for(int j=bounds.fromj; j<=bounds.toj; j++) if(map.flag[i][j] == 0) {
			double w = map.getS2N(i,j);
			index.x += w * i;
			index.y += w * j;
			sumw += Math.abs(w);
		}
		index.scale(1.0/sumw);
		map.grid.getCoords(index, coords);
	}
	
	public double spread(AstroMap map) {
		Bounds bounds = getBounds(map);
		Vector2D centerIndex = getIndex(map.grid);
		double sumd2 = 0.0, sumw = 0.0;
		
		for(int i=bounds.fromi; i<=bounds.toi; i++) for(int j=bounds.fromj; j<=bounds.toj; j++) if(map.flag[i][i] == 0) {
			double w = map.getS2N(i,j);
			double d = (i - centerIndex.x);
			sumd2 += w * d * d;
			d = (j - centerIndex.y);
			sumd2 += w * d * d;
			sumw += Math.abs(w);
		}
		
		return Math.sqrt(sumd2 / sumw);
	}
	
}
