/* 
% 
%   Set of utilities image processing functions for the ImageDecorrelationAnalysis plugin
%
% ---------------------------------------
%
% A detailled description of the method can be found in : 
% "Descloux, A., K. S. Grußmayer, and A. Radenovic. "Parameter-free image 
% resolution estimation based on decorrelation analysis."
% Nature methods (2019): 1-7."
%
%   Copyright © 2018 Adrien Descloux - adrien.descloux@epfl.ch, 
%   École Polytechnique Fédérale de Lausanne, LBEN/LOB,
%   BM 5.134, Station 17, 1015 Lausanne, Switzerland.
%
%  	This program 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.
%
%  	This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
*/

import ij.ImageStack;
import ij.process.FHT;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;

public class utils {
	
	
	// This class contains generic utility Image Processing functions required for the analysis

	
	static public FloatProcessor getPreprocessedImage(FloatProcessor im) {
		
		// make it square via zero-pad (mandatory for fht.transform)
		int newSize = Math.max(im.getWidth(), im.getHeight());

		newSize = (int)Math.pow(2, Math.ceil(Math.log((double)newSize)/Math.log(2.0)));
		
		float[] pixelsIn = (float[]) im.getPixels();
		
		// Since the image is already apodized, the padding value is given by any edge pixel
		float meanIm = pixelsIn[0]; 
		
		float[] pixelsOut = new float[newSize*newSize];
		for (int k = 0 ; k < newSize*newSize; k++) { pixelsOut[k] = meanIm; }
		
		int ox = (int) (newSize-im.getWidth())/2;
		int oy = (int) (newSize-im.getHeight())/2;

		int xIn = 0, xOut = ox;
		int yIn = 0, yOut = oy;
		for (xIn = 0, xOut = ox ; xIn <im.getWidth(); xIn++,xOut++)
			for (yIn = 0,yOut = oy; yIn < im.getHeight(); yIn++,yOut++) 
				pixelsOut[yOut*newSize +xOut] = pixelsIn[yIn*im.getWidth() + xIn];
			
		return new FloatProcessor(newSize,newSize,pixelsOut,null);
	}
	
	static public FloatProcessor apodizeEdges(FloatProcessor im) {

		float[] pin = new float[im.getWidth()*im.getHeight()];
		pin = (float[]) im.getPixels();
		
		double dist = 0;
		int off = 20; // number of pixels used for smooth apodization
		double r = im.getWidth()/2-off;
		
		float edgeMean = 0f;
		int count = 0;
		for (int x = 0; x < im.getWidth(); x++)
			for (int y=0;y < im.getHeight(); y++) {
				dist = (x-im.getWidth()/2)*(x-im.getWidth()/2) + (y-im.getHeight()/2)*(y-im.getHeight()/2);
				if (dist > r*r) {
					edgeMean += pin[y*im.getWidth() + x];
					count++;
				}
			}
		edgeMean = edgeMean/count;
		double x0 = 0;
		double y0 = 0;

		for (int x = 0 ; x < im.getWidth(); x++)
			for (int y = 0 ; y < im.getHeight() ; y++) {
			
				x0 = Math.abs(x-im.getWidth()/2);
				y0 = Math.abs(y-im.getHeight()/2);
				if (Math.abs(x0-im.getWidth()/2) <= off | Math.abs(y0-im.getHeight()/2) <= off) { 
					double d = Math.min(Math.abs(x0-im.getWidth()/2), Math.abs(y0-im.getHeight()/2));
					double c = (Math.cos(d*Math.PI/off-Math.PI)+1)/2; 
					
					pin[y*im.getWidth() + x] = (float) (c*(pin[y*im.getWidth() + x]-edgeMean))+edgeMean;
				}
				else if (Math.abs(x-im.getWidth()/2) > im.getWidth()/2 & Math.abs(y-im.getHeight()/2) > im.getHeight()/2) {
					pin[y*im.getWidth() + x] = (float) edgeMean;
				}
			}
		
		FloatProcessor out = new FloatProcessor(im.getWidth(),im.getHeight(),pin);
		
		return out;
	}
	
	static public FloatProcessor substract(FloatProcessor I1, FloatProcessor I2) {
		float[] pixelsI1 = (float[]) I1.getPixels();
		float[] pixelsI2 = (float[]) I2.getPixels();

		float[] pixelsOut = new float[pixelsI1.length];

		for (int k = 0; k < pixelsI1.length; k++)
			pixelsOut[k] = pixelsI1[k]-pixelsI2[k];
		
		return new  FloatProcessor(I1.getWidth(),I1.getHeight(),pixelsOut);
	}
	
	static public FloatProcessor[] normalizeFFT(FloatProcessor[] Ir) {
		
		float[] pixelsInR = (float[]) Ir[0].getPixels();
		float[] pixelsInI = (float[]) Ir[1].getPixels();
		
		FloatProcessor[] out = new FloatProcessor[2];
		out[0] = new FloatProcessor(Ir[0].getWidth(),Ir[0].getHeight());
		out[1] = new FloatProcessor(Ir[0].getWidth(),Ir[0].getHeight());
				
		float[] pixelsOutR = new float[pixelsInR.length];
		float[] pixelsOutI = new float[pixelsInR.length];
		
		for (int k = 0 ; k < pixelsInR.length; k++) {
			float mag = (float) Math.sqrt(pixelsInR[k]*pixelsInR[k] + pixelsInI[k]*pixelsInI[k]);
			if (mag != 0) { // avoid dividing by zero
				pixelsOutR[k] = pixelsInR[k]/mag;
				pixelsOutI[k] = pixelsInI[k]/mag;
			}
		}
		
		out[0].setPixels(pixelsOutR);
		out[1].setPixels(pixelsOutI);
		return out;
	}
	
	static public FloatProcessor[] getComplexFFT(ImageProcessor ip) {
		FHT fht = new FHT(ip);
		fht.transform();

		FloatProcessor[] ret = new FloatProcessor[2];

		ImageStack stack1 = fht.getComplexTransform();
		ret[0] = ((FloatProcessor) stack1.getProcessor(1));
		ret[1] = ((FloatProcessor) stack1.getProcessor(2));

		return ret;
	}
	
	// map the input value val, which belongs to [valmin,valmax] to the range [mapmin,mapmax]
	static public double linmap(double val, double valmin, double valmax, double mapmin, double mapmax) {
		double out = 0;
		// range is now [0,1]
		out = (val-valmin)/(valmax-valmin);
		// range is now [mapmin,mapmax]
		out = out*(mapmax-mapmin) + mapmin;
		return out;
	}
	
	static public double clamp(double val, double valmin, double valmax) {
		if (val < valmin)
			return valmin;
		else if (val > valmax)
			return valmax;
		else
			return val;
	}
	
}
