(* ::Package:: *)

(* ::Code:: *)
(*(************************************************************************)*)
(*(*								FMSfit									*)*)
(*(*			 Program for fluorescence spectral simulations				*)*)
(*(*																		 *)*)
(*(* 					Iztok Urbancic & Ajasja Ljubetic					*)*)
(*(*							  January 2013			                   *)*)
(*(*                                                                          *)*)
(*(*		                                                                  *)*)
(*(*		Modification for additional hemoglobin oxygenation fitting		*)*)
(*(*                                                                           *)*)
(*(*                               Rok Podlipec                                 *)*)
(*(*                                2018-2019                                    *)*)
(*(************************************************************************)*)
(**)


(* ::Section:: *)
(*(* functions for computing and debugging *)*)


(* *** functions to distribute compiled functions between kernels *** *)

ClearAll[CompiledFunctionNames];
CompiledFunctionNames[pattern_] :=
    Select[Names["Global`*"],MatchQ[Evaluate[Symbol@#],_CompiledFunction]&];

ClearAll[DistributeCompiledFunction];
DistributeCompiledFunction[name_] :=
    With[ {copy = Symbol[name]},
        ParallelEvaluate[If[ MatchQ[Evaluate[Symbol@name],_Symbol],(*Print["Setting ",name];*)
                             Evaluate[Symbol[name]] = copy
                         ]]
    ];

ClearAll[DistributeCompiledFunctions];
DistributeCompiledFunctions[pattern_] :=
    DistributeCompiledFunction/@CompiledFunctionNames[pattern];


(* *** functions for debugging *** *)

ClearAll[Puts];
SetAttributes[Puts,HoldAll];
$Verbose = True;
$VerboseLevel = 3;
Puts[str_,n_:1] :=
    If[ $Verbose &&($VerboseLevel>=n),
        Print[ToString@StandardForm@str]
    ];
Puts[str__,n_:1] :=
    If[ $Verbose&&($VerboseLevel>=n),
        Print[ToString@StandardForm@Row[{str}," "]]
    ];

ClearAll[PutsVar];
SetAttributes[PutsVar,{HoldAll,Listable}]
PutsVar[var_Symbol,n_Integer:1] :=
    Puts[Row@{SymbolName@Unevaluated@var," = ",Short[Evaluate@var,10]},n];

ClearAll[debug];
SetAttributes[debug,HoldAll];
debug[code_] :=
    Internal`InheritedBlock[{Message},Module[ {inMessage},
                                          Unprotect[Message];
                                          Message[args___]/;!MatchQ[First[Hold[args]],_$Off] :=
                                              Block[ {inMessage = True},
                                                  Print[{Shallow/@Replace[#,HoldForm[f_[___]]:>HoldForm[f],1],Style[Map[Short,Last[#],{2}],Red]}&@Drop[Drop[Stack[_],-7],4]];
                                                  Message[args];
                                                  Throw[$Failed,Message];
                                              ]/;!TrueQ[inMessage];
                                          Protect[Message];
                                      ];
                                      Catch[StackComplete[code],Message]];
                                      
$Line



(* ::InheritFromParent:: *)
(**)


(* ::Section:: *)
(*(* IMPORT *)*)


(* 
ImportImages: import all necessary data about the experiment: name, images, reference images, lambda range, exposure time ...

input: 	folder 					.. full path where to set the file-opening 
		lctfCorrectionFileName	.. full name of the file with LCTF correction curve (parameter can be omitted)
		
output:	experimentName 			.. name of the .dsc file without file extension  
		images					.. 3D matrix with 16-bit data from all images in the spectral series; sorted by wavelength, as saved
		intensityData			.. data of all spectral images superimposed, normalized to the maximal value
		lambdaSample			.. list of wavelengths as sampled in the experiment
		timeExposure			.. exposure time of each image in seconds
		randomSampling			.. True for experiment with stochastic wavelength sampling, False for linear
		experimentData			.. a table of some info about experiments - currently not used
		importFolder			.. full path of the folder with the experiment; returned to be used as input at next ImportImages
		bleachRefMode			.. True for experiment with a reference wavelength for bleaching correction (used with linear sampling)
		bleachRefImages			.. 3D matrix with 16-bit data from images in the bleaching reference series; sorted by time taken
		bleachRefLambda			.. wavelength at which bleaching reference images were recorded
		bleachRefFrameNumbers	.. consecutive numbers of bleaching reference images in the whole series, e.g. {1,4,7,11,16,23} (not always equally spaced)
		"BatchModeQ"                  .. option: True/False to control behavior in batch analysis mode
 
*)

ClearAll[ImportImages];
Options[ImportImages] = {"BatchModeQ"->False};
ImportImages[folder_,filename_,lctfCorrectionFileName_:FileNameJoin[{NotebookDirectory[],"lib","LctfCamera.dat"}],opts:OptionsPattern[]] := 
   Block[ {importPath,importFolder,experimentName,dsc,expType,dataColumn,
   	lambda1,lambda2,lambdaStep,randomSampling,lambdaSample,dark,timeExposure,
   	imageNames,images,intensityData,lctfCorrection,correctionFactors,
   	bleachRefMode,bleachRefLambda,bleachRefNames,bleachRefImages,bleachRefLambdas,bleachRefFrameNumbers,bleachRefCorrectionFactor,
   	experimentData},
		
		(*if not given, ask for filename*)
		If[filename=="",		
			importPath = SystemDialogInput["FileOpen",{folder,{"Data files"->{"*.dsc"}}},WindowTitle->"Select .dsc file"];
			If[ToString@importPath=="$Canceled", Return[]];
			,
			importPath = FileNameJoin[{folder,filename}];
		];

		importFolder = FileNameTake[importPath, {1, -2}];
		experimentName = StringDrop[FileNameTake[importPath, -1], -4];
		(*Print[experimentName];*)
		(*name = FileNameTake[path,-1];*)
        
        (*import .dsc file with info about experiment*)
        dsc = Import[importPath,"Table"];
        
		(*check if the loaded file belongs to a spectral experiment; if not, warn the user (unless in batch analysis mode) and end the importing routine*)
		expType=With[{dataRow=Position[dsc,"ExpTyp"][[1,1]]}, dsc[[dataRow,3;;]]];
		If[Flatten[StringPosition[expType,"spectrum",IgnoreCase->True]]=={},
			If[!OptionValue@"BatchModeQ",CreateDialog[{TextCell["This is not a spectral experiment and cannot be fitted."],DefaultButton[]}]];
			Return[]
		];

		(*column in which the data appear; changed in FMS+ 2.8 when units were appended to each row*)
		dataColumn=If[dsc[[7,-1]]<=2.7,-1,-2];
		
		lambda1=With[{dataRow=Position[dsc,"ExpStr"][[1,1]]}, dsc[[dataRow,dataColumn]]];
		lambda2=With[{dataRow=Position[dsc,"ExpEnd"][[1,1]]}, dsc[[dataRow,dataColumn]]];
		lambdaStep=With[{dataRow=Position[dsc,"ExpStp"][[1,1]]}, dsc[[dataRow,dataColumn]]];
		bleachRefLambda=With[{dataRow=Position[dsc,"WvlRef"][[1,1]]}, dsc[[dataRow,-1]]];				
		timeExposure=With[{dataRow=Position[dsc,"ExpTim"][[1,1]]}, dsc[[dataRow,-1]]];		
		dark=With[{dataRow=Position[dsc,"PrABkg"][[1,1]]}, dsc[[dataRow,-1]]];		
            (*tweezers can add considerable unhomogenious background, which cannot be subtracted throughout the image;
                       Andor camera clamps the offset; *)
           
		If[dsc[[7,-1]]<=2.7(*file version*),
        	(*old dsc files*)
            randomSampling = If[ DeleteDuplicates[NumberQ/@Flatten[dsc[[13]]]]=={True},
                                 True,
                                 False
                             ];
            lambdaSample = If[ randomSampling,
                               Flatten[dsc[[13]]],
                               Range[lambda1,lambda2,lambdaStep]
                           ];
            ,
   		(*new dsc files*)
            randomSampling = With[{dataRow=Position[dsc,"ExpTyp"][[1,1]]}, dsc[[dataRow,3]]=="Stochastic"];
            lambdaSample = If[ randomSampling,
							With[{dataRow=Position[dsc,"RndWvl"][[1,1]]}, Flatten[dsc[[dataRow,3;;]]]],
                            Range[lambda1,lambda2,lambdaStep]
                           ];
		];
		

		(*load the file with LCTF transmittance to correct the data from images*)
        lctfCorrection = Import[lctfCorrectionFileName, "Table"];
        lctfCorrection = Transpose[{lctfCorrection[[All, 1]], N[MeanFilter[lctfCorrection[[All, 2]], 5]]}];
        (*select the values at wavelengths that appear in the experiment*)
        correctionFactors = lctfCorrection[[Sort[lambdaSample] - 399, 2]];
        
        (*import images, convert them into a numerical 3D matrix, subtract offset (dark) and correct for LCTF transmittance*)
        imageNames = FileNames[experimentName~~NumberString~~".tif",{importFolder}];
        imageNames = Select[imageNames, StringLength[FileNameTake[#,-1]] == StringLength[experimentName]+7 &];
		images = 2^16 (ImageData[Import[#]]&/@imageNames)-dark;
        images /= correctionFactors;
        (*get intensity image as superposition of all spectral images in the series*)
        intensityData = Total[images]/Length[images]/2^16;
        (*intensityData -= Min[Flatten[intensityData]];
        intensityData /= Max[Flatten[intensityData]];*)
        (*Print[GraphicsRow[{Show[ArrayPlot[Total[images],ColorFunction->GrayLevel],Frame->True,FrameLabel->{"x [pix]","y [pix]"},LabelStyle->{FontSize->14},FrameTicks->{{All,None},{All,None}},ImageSize->400],ListLinePlot[Transpose[{Sort[lambdaSample],images[[All,250,250]]}],PlotRange->{0,All},Frame->True,FrameLabel->{"\[Lambda] [nm]","I  [a.u.]"},LabelStyle->{FontSize->14}]},ImageSize->600]];*)
        
        (*linear bleaching reference:*)
        bleachRefMode = False;
		{bleachRefImages,bleachRefFrameNumbers}={{0},{0}};
        If[randomSampling == False,
        	bleachRefNames = FileNames[experimentName <> "_blc_" ~~ NumberString ~~ ".tif",{importFolder}];
        	If[bleachRefNames != {},
				bleachRefMode = True;
				bleachRefImages = 2^16 (ImageData[Import[#]] & /@ bleachRefNames) - dark;
				bleachRefLambdas = ToExpression[StringTake[#, {-7, -5}]] & /@ bleachRefNames;
				bleachRefFrameNumbers = Flatten[Position[lambdaSample, #] & /@ bleachRefLambdas] + 1;
				bleachRefFrameNumbers += Range[0, 5];
				bleachRefFrameNumbers[[1]] -= 1;
				bleachRefCorrectionFactor = lctfCorrection[[bleachRefLambda - 399, 2]];
				bleachRefImages /= bleachRefCorrectionFactor;
			];
		];
      
		(*Print[Timing[Print[{experimentName,images[[1;;2]]}]]];*)
        
 		Return[{experimentName,images,intensityData,lambdaSample,timeExposure,randomSampling,experimentData,importFolder,
        	{bleachRefMode,bleachRefImages,bleachRefLambda,bleachRefFrameNumbers}}];
    ];




(* ::Section:: *)
(*(* BACKGROUND *)*)


(* 
BackgroundSubtractFromImage: subtracts the average background spectrum from the selected part of the image (blue rectangle)

input:  images 					.. 3D matrix with 16-bit data from all images in the spectral series; sorted by wavelength, as saved
		bleachRefImages			.. 3D matrix with 16-bit data from images in the bleaching reference series; sorted by time taken
		cropBackground              .. {{xmin,ymin},{xmax,ymax}} of the region to calculate noise from [pix]
		bleachRefLambdaPosition     .. position of bleachRefLambda in the lambda range
		
output: imagesBackgroundSubtracted            .. as images (above), with background signal (average spectrum) subtracted
		bleachRefImagesBackgroundSubtracted   .. as bleachRefImages (above), with background signal (average intensity at reference lambda) subtracted
		backgroundSpectrum                    .. list of intensities over lambda, averaged over the selected part of the image
*)
ClearAll[BackgroundSubtractFromImage]
BackgroundSubtractFromImage[images_,bleachRefImages_,cropBackground_,bleachRefLambdaPosition_]:=Block[{backgroundSpectrum,imagesBackgroundSubtracted,backgroundBleachRefLambda,bleachRefImagesBackgroundSubtracted},

	(*get spectrum from the blue rectangle and subtract it from all images*)
	backgroundSpectrum = SpectrumAvg[images,cropBackground];(*ordered by wavelength*)
	imagesBackgroundSubtracted = images - backgroundSpectrum;
	
	(*subtract the appropriate value of the background spectrum from bleaching reference images*)
	backgroundBleachRefLambda = backgroundSpectrum[[bleachRefLambdaPosition]];
	bleachRefImagesBackgroundSubtracted = bleachRefImages - backgroundBleachRefLambda;

	Return[{imagesBackgroundSubtracted,bleachRefImagesBackgroundSubtracted,backgroundSpectrum}];
];


(* 
BackgroundReset: reset images with subtracted background to original values

input:  imagesBackgroundSubtracted            .. as images (above), with background signal (average spectrum) subtracted
		bleachRefImagesBackgroundSubtracted   .. as bleachRefImages (above), with background signal (average intensity at reference lambda) subtracted
		backgroundSpectrum                    .. list of intensities over lambda, averaged over the selected part of the image
		bleachRefLambdaPosition               .. position of bleachRefLambda in the lambda range
		
output: images 					.. *original* 3D matrix with 16-bit data from all images in the spectral series; sorted by wavelength, as saved
		bleachRefImages			.. *original* 3D matrix with 16-bit data from images in the bleaching reference series; sorted by time taken
		backgroundSpectrumReset    .. 0 
*)
ClearAll[BackgroundReset]
BackgroundReset[imagesBackgroundSubtracted_,bleachRefImagesBackgroundSubtracted_,backgroundSpectrum_,bleachRefLambdaPosition_]:=Block[{images,bleachRefImages,backgroundBleachRefLambda,backgroundSpectrumReset},

	(*add the background spectrum back to the data*)
	images = imagesBackgroundSubtracted + backgroundSpectrum; (*ordered by wavelength*)
	
	(*add the appropriate value of the background spectrum to bleaching reference images*)
	backgroundBleachRefLambda = backgroundSpectrum[[bleachRefLambdaPosition]];
	bleachRefImages = bleachRefImagesBackgroundSubtracted + backgroundBleachRefLambda;
	
	(*reset the background spectrum*)
	backgroundSpectrumReset = 0;

	Return[{images,bleachRefImages,backgroundSpectrumReset}];
];





(* ::Section:: *)
(*(* IMAGE MANIPULATION *)*)


(*
CropImages: crop 512x512 images to a subsection

input: 	images 			.. 3D matrix of data from a set of 512x512 images
		{xmin,ymin} 	.. bottom left corner coordinates [pix]; y-axis reversed because ArrayPlot plots rows top-down
		{xmax,ymax} 	.. top right corner coordinates [pix]; 

output:	images[[...]]	.. 3D submatrix		 
*)
ClearAll[CropImages]
CropImages[images_,{{xmin_,ymin_},{xmax_,ymax_}}] := images[[All,512-ymax+1;;512-ymin+1,xmin;;xmax]];


(*
SNRfromImages: calculate mean signal-to-noise (SNR) ratio from a cropped part of a series of images

input:	images 			.. 3D matrix of data from a set of 512x512 images
		crop 			.. {{xmin,ymin},{xmax,ymax}} of the region to calculate signal from [pix]
		cropBackground 	.. {{xmin,ymin},{xmax,ymax}} of the region to calculate noise from [pix]

output: noise			.. standard deviation of signal from cropBackground, avereged over all images in the series
		signal/noise	.. average of max signals from points in cropped region / noise
*)
ClearAll[SNRfromImages]
SNRfromImages[images_,crop_,cropBackground_] :=
    Block[ {imagesCrop,backgroundCrop,signal,noise},
        imagesCrop = CropImages[images,crop];
        backgroundCrop = CropImages[images,cropBackground];
        noise = Mean[StandardDeviation[Flatten[#]]&/@backgroundCrop];

		(*if image doesn't have dark background spots noise is rather stdev of the spectra tails*)
		
        signal = Mean[Max[Flatten[#]]&/@imagesCrop];
        Return[{signal/noise,noise}];
    ];

(*
SpectrumAvg: get average spectrum from a cropped subregion of images

input: 	images 			.. 3D matrix of data from a set of 512x512 images
		crop 			.. {{xmin,ymin},{xmax,ymax}} of the region to calculate signal from [pix]

output: list of intensities 
*)
ClearAll[SpectrumAvg]
SpectrumAvg[images_,crop_] := Mean[Flatten[#]]&/@CropImages[images,crop]
    

(*
SpectraFromImages: get a list of spectra from all {x,y} points in the image series (resolution reduced) 
				with some additional info, needed for optimization

input: 	images						.. 3D matrix of data from a set of 512x512 or cropped images
		cropXY						.. {{xmin,ymin},{xmax,ymax}} of the region to be analized [pix]
		lambdaShuffle				.. list of numbers that converts ordered wavelength list into sampled row: lambdaSample = Sort[lambdaSample][[lambdaShuffle]]  
		nAvg						.. averaging for nAvg pix in every direction from the central point -> over (2 nAvg + 1) x (2 nAvg + 1) pix 
		noise						.. standard deviation of background signal
		signalThreshold				.. intensity threshold to filter points with signal < signalThreshold * noise
		bleachRefMode				.. True for experiment with a reference wavelength for bleaching correction (used with linear sampling)
		bleachRefImages				.. 3D matrix with 16-bit data from images in the bleaching reference series; sorted by time taken
		numberOfVariableParameters	.. number of parameters to be optimized in fitting routine

output:	out 						.. a list of {{x,y,0,0,0,0,...}, spectrum, SNR, fitSwitch, bleachSeries, variableParameterStartingPoint} from all points
		 {x,y,0,0,0,0,...}			.. coordinates and places for other parameters
		 spectrum					.. list of intensities through images at position {x,y} - IN THE ORDER AS ACQUIRED!!!
		 SNR						.. signal-to-noise in that point
		 fitSwitch					.. switch to perform fit or not: True if signal above signalThreshold*noise, False if below
		 bleachSeries				.. list of intensities through bleachRefImages
		 variableParameterStartingPoint .. list of places for starting values for variable parameters
*)
ClearAll[SpectraFromImages];
SpectraFromImages[images_,cropXY_,lambdaShuffle_,nAvg_,noise_,signalThreshold_,{bleachRefMode_,bleachRefImages_},numberOfVariableParameters_]:=
   Block[{xy,points,imagesAvg,spectra,snr,out,fitSwitch,bleachSeries,variableParameterStartingPoint},
   		(*define points where to take spectra from*)
   		points = Tuples[Range[nAvg+1,#-(nAvg),2 nAvg+1]&/@(Dimensions[images[[1]]])];
   		(*get their xy coordinates*)
   		xy = {cropXY[[1,1]]-1,cropXY[[2,2]]+1}+{1,-1}Reverse[#1]&/@points;
   		xy = PadRight[xy, {Length[xy],16}, 0.];

   		(*extract spectra*)
   		spectra = (Mean[Flatten[#]]&/@images[[lambdaShuffle,#1[[1]]-nAvg;;#1[[1]]+nAvg,#1[[2]]-nAvg;;#1[[2]]+nAvg]]&)/@points;
   		(*calculate SNR for every spectrum*)
      	snr = (Mean[Sort[#1][[-3;;All]]]/noise&)/@spectra;

   		(*check if intensity above threshold*)
       	fitSwitch = #[[3]]>signalThreshold&/@spectra;
       	(*bleaching reference series*)
   		bleachSeries = If[bleachRefMode, 
   			(Mean[Flatten[#]]&/@bleachRefImages[[All,#1[[1]]-nAvg;;#1[[1]]+nAvg,#1[[2]]-nAvg;;#1[[2]]+nAvg]]&)/@points,
   			(*else*)ConstantArray[Null,Length[points]]];
   		(*prepare places for variableParameterStartingPoint*)
   		variableParameterStartingPoint = ConstantArray[0.,{Length[spectra],numberOfVariableParameters}];
   		(*output*)
   		out = Transpose[{xy,spectra,snr,fitSwitch,bleachSeries,variableParameterStartingPoint}];
   		(*order by rows and columns*)
        out = SortBy[out,{#[[1,2]]&,#[[1,1]]&}(*Total[#1[[2]]]>Total[#2[[2]]]&*)];

   		Return[out];
	];





(*
AlignImages: align images in the stack to compensate for particle movement

input: 	images			.. 3D matrix of data from a set of images (arbitrary dimensions)
		nAvg			.. averaging for nAvg pix in every direction from the central point -> over (2 nAvg + 1) x (2 nAvg + 1) pix 
		"AlignMethod"	.. option: built-in numerical method for image alignment: "Translation", "Rigid" (translation + rotation), ... 

output: imagesAlign		.. 3D matrix of aligned images (smoothened!) 
		imageIntCompareAlign .. aligned intensity image(*column of two small intensity images to compare before/after alignment*)
*)
ClearAll[AlignImages];
Options[AlignImages]={"AlignMethod"->"Translation"};
AlignImages[images_,nAvg_,opts:OptionsPattern[]]:=Block[{imageIntOrig,imagesAlign,imageIntAlign,imageIntCompareAlign,imageIntensityMax,imageBrightestIndex},
	(*generate intensity image from original slices*)
	imageIntOrig=Show[ArrayPlot[Total[images],ColorFunction->GrayLevel,PlotRange->All],(*PlotLabel->Style["original",14],*)Frame->False,ImageSize->150];
	(*smoothen images*)
	imagesAlign=ParallelMap[MeanFilter[#,nAvg]&,images];
	(*find the brightest image - average of top 4% of pixels, neglecting 1% (potential hot pixels)*)
	imageIntensityMax=With[{numberOfPixels=Times@@Dimensions[imagesAlign[[1]]]},
			ParallelMap[Mean[Sort[Flatten[#]][[Round[.95 numberOfPixels];;Round[.99 numberOfPixels]]]]&,imagesAlign]
		];
	imageBrightestIndex=Position[imageIntensityMax,Max[imageIntensityMax]][[1,1]];
	(*convert from ImageData back to images*)
	imagesAlign=Image[#,"Bit16"]&/@imagesAlign;
	(*align all with first*)
	imagesAlign=With[{imageBrightest=imagesAlign[[imageBrightestIndex]]},
					ParallelMap[Check[ImageAlign[imageBrightest,#,"Transformation"->OptionValue@"AlignMethod",Background->Black],#]&,imagesAlign]];
	(*convert back to ImageData*)
	imagesAlign=2^16ImageData/@imagesAlign;
	(*fill cropped areas of aligned images with some background signal*)
	(*imagesAlign=FillCroppedAreas[imagesAlign];*)
	(*generate intensity image from aligned slices*)
	imageIntAlign=Show[ArrayPlot[Total[imagesAlign],ColorFunction->GrayLevel,PlotRange->All],(*PlotLabel->Style["aligned",14],*)Frame->False,ImageSize->{300}(*150*)];
	(*compare intensity images - original and aligned*)
	imageIntCompareAlign=imageIntAlign;(*GraphicsColumn[{imageIntOrig,imageIntAlign},Spacings->1,ImageSize->{300}];*)
	Return[{imagesAlign,imageIntCompareAlign}];
];


(* ::Section:: *)
(*(* BLEACHING *)*)


(*
BleachFitLinearRef: calculates bleaching rates from reference data series: y/y0 = Exp[-rate x] --> rate = -Log[y/y0]/x ; 
				averaged over all reference points
				
input:	bleachSeries 			.. list of bleaching series (time series of intensity at every {x,y} at reference wavelength)
		bleachRefFrameNumbers 	.. consecutive numbers of bleaching reference images in the whole series, e.g. {1,4,7,11,16,23} (not always equally spaced)

output:	bleachRates				.. list of bleaching rates in all points [frames^-1] 			
*)
ClearAll[BleachFitLinearRef]
BleachFitLinearRef[bleachSeries_,bleachRefFrameNumbers_] :=
    Block[ {bleachSeriesNorm, bleachRates, bleachRateAvg, out},
        bleachSeriesNorm = bleachSeries;
        bleachSeriesNorm /= bleachSeriesNorm[[All,1]];
        bleachRates = (# / (bleachRefFrameNumbers[[2;;]]-1)) & /@ (-Log[bleachSeriesNorm[[All,2 ;;]]]);
        bleachRates = Re[bleachRates];
        bleachRates = Clip[bleachRates,{0,Infinity}];
        bleachRates = Mean /@ bleachRates;
        Return[bleachRates];
    ]


(* ::Subsection:: *)
(*Spectral data of hemoglobin absorbance*)


(*ClearAll[hemSpectra]
hemSpectra[lambda_]:=Block[{fileName,data,dataY},
SetDirectory[NotebookDirectory[]];
fileName=FileNames[];
(*data=FileNames[__~~".dat","Data"];*)
data=Import[fileName[[1]]];
dataY=data[[Flatten[Position[data,lambda]][[1]],2]]

Return[dataY]
]
*)


(* ::Section:: *)
(*(* MODEL SPECTRUM FUNCTIONS *)*)


(*
SpectrumDescriptiveC: compiled function with spectral model (for derivation see Iztok's diary IU@FMS_01 p .108): 
					I[lambda] = Exp[ -Log[2] (Log[4 asymmetry/width (lambda-lambdaMax) + 1] / asymmetry /2 )^2 ] or 
								0 if lambda < (lambdaMax - width/(4 asymmetry))

input: 	par 		.. spectral parameters {lambdaMax, width, asymmetry} (positions 5-7 and 10-12 in the complete set of parameters)
		 lambdaMax	.. peak position [nm]
		 width		.. FWHM [nm]
		 asimmetry	.. ((part of FWHM right of lambdaMax) - (part of FWHM left of lambdaMax)) / FWHM
		 lambda		.. wavelength [nm] to calculate I[lambda] 

output: intensity (a number)
*)
ClearAll[SpectrumDescriptiveC]
SpectrumDescriptiveC = Compile[{{par,_Real,1},{lambda,_Real,0}},
	If[ lambda > par[[1]]-par[[2]]/(4. par[[3]]),
			If[par[[1]]<600,
				Exp[-Log[2.](Log[4.par[[3]]/par[[2]] (lambda-par[[1]])+1.]/par[[3]]/2.)^2],
				
				(*next functions represent oxygenated and de-oxygenated hemoglobin absorbance spectra*)
				(*the oxy-hemoglobin absorbance curve is composed of two Gauss functions with Lorenz tails and additional linear correction*)

				N[par[[4]]*(0.93*(Exp[-((lambda-543)/100.)^2/0.034]*1./(((lambda-543)/100.)^2+1.)+ Exp[-((lambda-576)/100.)^2/0.02]*1./(((lambda-576)/100.)^2+1.)) + (1.-0.93)*(-1.5*(lambda-500.)/100.+1.09)) + 
				(1.-par[[4]])*(Exp[-((lambda-559.1)/100.)^2/0.143]*1./(((lambda-559.1)/100.)^2+1.) + (-0.93*(lambda-500.)/100.+0.52))]
				
]	,		

        0
    	],
    RuntimeAttributes->{Listable}
    ];
DistributeDefinitions[SpectrumDescriptiveC];

(*
SpectrumDescriptiveTableC: calculates model spectrum intensities over a list of wavelengths for given parameters 
						(SpectrumDescriptiveC is Listable, but somehow does not work without this )

input: 	param 		.. spectral parameters {lambdaMax, width, asymmetry} (positions 6-8 and 12-14 in the complete set of parameters) - as above
		lambdaRange	.. list of wavelengths
		
output: list of intensities in the same order as input wavelengths
*)
ClearAll[SpectrumDescriptiveTableC]
SpectrumDescriptiveTableC = 
        Compile[{{param, _Real, 1}, {lambdaRange, _Real, 1}},
	          (SpectrumDescriptiveC[param, #]) & /@ lambdaRange
       ,CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}];
       
DistributeDefinitions[SpectrumDescriptiveTableC];

(*
SpectraSum: calculates model spectrum as a sum of two components;
In case of characterization of oxygenated hemoglobin, second component is subtracted from the first component

input: 	parameters 		.. full list for 2 components 
							{x,y,intensity,intensity1,portion1,lambdaMax1,width1,asimmetry1,bleachingRate1,intensity2,p2,lambdaMax2,w2,a2,b2,HemOxy,chi2};
							bleaching rates in [frames^-1]
		lambdaSample	    .. list of wavelengths, ordered by TIME OF ACQUISITION (shuffled in the case of stochastic sampling)
		
output:	spectrum		.. list of intensities (only) in the same order as lambdaSample 
*)
ClearAll[SpectraSum];
SpectraSum[parameters_,lambdaSample_] :=
    Block[ {spectrumSimComponents,bleachingConstants,bleachingFactors,spectrum},
        (*calculate the two spectra*)

		If[Length[parameters]==16,parameters=Append[parameters,0],None];

        spectrumSimComponents = Chop[SpectrumDescriptiveC@@{#,lambdaSample}&/@{parameters[[6;;9]],Pick[parameters,{0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,0},1](*parameters[[12;;14]]*)}];
        (*prepare bleaching correction for both components*)
        bleachingFactors = Exp[- Range[0,Length[lambdaSample]-1] #]&/@parameters[[{9,15}]];
        (*weight the two spectra by bleaching factors and spectral portions and then sum them together*)
        spectrumSimComponents = MapThread[Times,{spectrumSimComponents,bleachingFactors,parameters[[3]]*parameters[[{5,11}]]}];
        spectrum = Total[spectrumSimComponents];
        Return[spectrum];
    ];
DistributeDefinitions[SpectraSum];

(*compiled version of SpectraSum*) 
ClearAll[SpectraSumC];
SpectraSumC = Compile[{{par, _Real,1}, {lambdaSample,_Real,1}(*,{hemSpectra,_Real,1}*)}, 
				 Block[ {spectrumSimComponents={{0.},{0.}},bleachingConstants,bleachingFactors,par2},


                       If[Length[par]==16,
							par2=Append[par,0];
						(*calculate the two spectra*)
                        spectrumSimComponents = 
							SpectrumDescriptiveTableC@@ {#, lambdaSample} &/@ {par2[[6;;9]],Pick[par2,{0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,0},1](*par[[12;;14]]*)},
						spectrumSimComponents = 
							SpectrumDescriptiveTableC@@ {#, lambdaSample} &/@ {par[[6;;9]],Pick[par,{0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,0},1](*par[[12;;14]]*)};
                            (*SpectrumDescriptiveTableC @@ {#, lambdaSample} & /@ {par[[6;;8]], par[[12;;14]]};*)];

			
				        (*prepare bleaching correction for both components*)
				        bleachingFactors = Exp[- Range[0,Length[lambdaSample]-1] #]&/@par[[{9,15}]];
				        (*weight the two spectra by bleaching factors and spectral portions and then sum them together*)
				        spectrumSimComponents = MapThread[Times,{spectrumSimComponents,bleachingFactors,par[[3]]*par[[{5,11}]]}];
			       
						spectrumSimComponents[[1]]-spectrumSimComponents[[2]]
				        (*Total[spectrumSimComponents]*)
				       ]                      
                      ,CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}];



(* ::Section:: *)
(*(* OPTIMIZATION *)*)


(*
ParametersOutsideLimitsQ: check if any of the parameters is outside their limits; used during optimization

input: 	par 		.. list of all variable paramters
		parLimits	.. list of all corresponding limits {{par1min,par1max}, ...}

output: result		.. True if any of par is outside, False if all within limits
*)                                      
ClearAll[ParametersOutsideLimitsQ]
ParametersOutsideLimitsQ = 
  Compile[{{par, _Real, 1}, {parLimits, _Real, 2}},
   Block[{i = 1, result = False},
    For[i = 1, i <= Length@par, i++,
     If[Not[parLimits[[i, 1]] <= par[[i]] <= parLimits[[i, 2]]], 
      result = True; Break[]]];
    result
    ]
   ];                                      


(*
CompiledChi2: compute the simulated spectrum and compare it to the experimental in terms of chi2; lower intensities filtered not to contribute much 

input:	par				.. list of all parameters; bleaching in [frames^-1]
		spectrumExp		.. list of experimental spectral intensities, ordered as lambdaSample
		lambdaSample	.. list of wavelengths, ordered BY TIME OF ACQUISITION (shuffled in the case of stochastic sampling)
		noise			.. standard deviation of background signal

output:	chi2			.. normalized (/noise^2) and reduced (/N spectral points) chi2
*)
ClearAll[CompiledChi2];
With[ {infinity = $MaxMachineNumber},
    CompiledChi2 = 
      Compile[{{par, _Real, 1}, {spectrumExp, _Real, 1}, {lambdaSample, _Real, 1}, {noise, _Real, 0}},
       Block[ {spectrumSim, filterWeights, totalFilterWeights, intensitySim, chi2, noise2},
           (*calculate simulated spectrum*)
           spectrumSim = SpectraSumC[par,lambdaSample];

		   (*noise2 = Sort[spectrumExp];*)
		   (*noise is optionally set to 1. if there is no dark spots of background in the image*)
			noise2=30.;
           (*calculate chi2; points weighted by intensities - a Gaussian filter damps the influence of data with signal close to noise*)
           filterWeights = 1. - Exp[-spectrumExp^2/2./(2. noise2)^2];
           totalFilterWeights = Total[filterWeights spectrumSim^2];
           intensitySim = 
            If[ totalFilterWeights > 0., (*then*)
                Total[filterWeights spectrumExp spectrumSim]/totalFilterWeights,
                (*else*)
                infinity
              ];
           chi2 = Total[filterWeights (intensitySim spectrumSim - spectrumExp)^2]/(*Length[spectrumExp]*)Total[filterWeights]/noise2^2;
           (*Print[{par[[{5,8}]],chi2}];*)
           chi2
       ]
       , {{filterWeights, _Real, 1}, {spectrumSim, _Real, 1}, {intensitySim, _Real, 1}, {totalFilterWeights, _Real, 1}}
       , CompilationOptions -> {"InlineCompiledFunctions" -> True, "ExpressionOptimization" -> True, "InlineExternalDefinitions" -> True}
       , RuntimeOptions -> {"CompareWithTolerance" -> False, 
         "EvaluateSymbolically" -> False}];
];
 
 
 
(*
CompiledLimitedChi2: returns large chi2 if parameters get out of the specified ranges; else calculates chi2 by calling CompiledChi2
					bleaching in [frames^-1] 

input:	par							.. list of all parameters
		spectrumExp					.. list of experimental spectral intensities, ordered as lambdaSample
		lambdaSample				.. list of wavelengths, ordered BY TIME OF ACQUISITION (shuffled in the case of stochastic sampling)
		noise						.. standard deviation of background signal
		variableParamtersPositions 	.. indices of variable parameters in the list of all parameters; e.g. {6,9} for variable {lambdaMax1,bleach1}
		variableParametersLimits	.. list of limits for variable parameters {{par1min,par1max},..}
		chi2Outside					.. a large value for chi2, if parameters get out of their limits, e.g. 100.

output:	chi2						.. normalized (/noise^2) and reduced (/N spectral points) chi2
*)  
ClearAll[CompiledLimitedChi2];
CompiledLimitedChi2 = 
      Compile[{{par, _Real, 1}, {spectrumExp, _Real, 1}, {lambdaSample, _Real, 1}, {noise, _Real, 0},
      	       {variableParametersPositions, _Integer, 1}, {variableParametersLimits, _Real, 2}, {chi2Outside, _Real, 0}},
       Block[ {chi2=chi2Outside},
(*       	 Print["par = ", par];
         Print["variableParametersPositions = ", variableParametersPositions];
         Print["variableParametersLimits = ", variableParametersLimits];
*)         If[	Not@ParametersOutsideLimitsQ[par[[variableParametersPositions]], variableParametersLimits],(*then*)
           		chi2=CompiledChi2[par,spectrumExp, lambdaSample, noise]; 	
         ];
         chi2	
       ]
       , CompilationOptions -> {"InlineCompiledFunctions" -> True, "ExpressionOptimization" -> True, "InlineExternalDefinitions" -> True}
       , RuntimeOptions -> {"CompareWithTolerance" -> False, "EvaluateSymbolically" -> False}];

CompiledChi2w[parameters:{_?NumericQ..},spectrumExp:{_?NumericQ..},lambdaSample:{_?NumericQ..},noise_?NumericQ] :=
    CompiledChi2w[parameters,spectrumExp,lambdaSample,noise] = CompiledChi2[parameters,spectrumExp,lambdaSample,noise];
DistributeDefinitions[CompiledChi2w];


(*
SpectrumBleachFit: runs the simulation for one spectrum to get all parameter values; 
				uncompiled version writes results (parameters,chi2) into a temporary file;
				bleaching in [frames^-1] 

input:	spectrum						.. list of experimental spectral intensities, ordered as lambdaSample
		lambdaSample					.. list of wavelengths, ordered BY TIME OF ACQUISITION (shuffled in the case of stochastic sampling)
		noise							.. standard deviation of background signal
		initialPoint					.. starting values for variable parameters
		fixedOrVariableParameters		.. list of all parameters; fixed with numbers, variable with symbols, e.g. {1,1,1,1,1,lm1,75,.24,b1,0,0,550,75,.24,0,HemOxy}
		variableParamtersPositions 		.. indices of variable parameters in the list of all parameters, e.g. {6,9} for variable {lm1,b1}
		variableParameterSymbols		.. list of symbols for variable parameters, e.g. {lm1,b1} 
		variableParametersLimits		.. list of limits for variable parameters {{par1min,par1max},..}, e.g. {{400,700},{0,1}};
		variableParameterRangeSizes		.. list of spans between the limits for variable parameters {par1max-par1min,..}, e.g. {300,1}
		variableParameterConstraints	.. list of formated constraints for variable parameters {{par1min < par1 < par1max},..}, e.g. {400<=lm1<=700,0<=b1<=1}
		performFit						.. fit switch; True if optimization is to be performed, False if not

output: appends a row with all parameter values and chi2 to a file $stream (global variable!) 
*)
ClearAll[SpectrumBleachFit];
SpectrumBleachFit[spectrum_, lambdaSample_, noise_, 
	initialPoint_, fixedOrVariableParameters_, variableParametersPositions_, variableParameterSymbols_, 
	variableParameterLimits_, variableParameterRangeSizes_, variableParameterConstraints_,
	performFit_:True] :=
    Module[ {initialPoints, chi2norm, solution, parameters},
        If[performFit, (*then*)
	        (*minimization of n parameters needs n+1 starting points;
	        one is the initialPoint, others are generated by varying only one parameter for 10% of the range (one along each parameter, +/- random); 
	        clipped to constraints*)
	        initialPoints = 
	        	Join[
	        		{initialPoint},
					Table[
						MapThread[
							Clip,
								{initialPoint+0.1 variableParameterRangeSizes[[i]] Sign[RandomReal[{-1.,1.}]] UnitVector[Length[variableParametersPositions],i]
								,variableParameterLimits}],
					{i,Length[variableParametersPositions]}]
					];
			
	        (* results in the form {chi2, param1, param2, ..} *)
	        Check[
	        	
				{chi2norm,solution} = NMinimize[
	            	Join[{CompiledChi2w@@{fixedOrVariableParameters,spectrum,lambdaSample,noise}},variableParameterConstraints],
	            	variableParameterSymbols,
	            	Method->{"NelderMead","PostProcess"->False,"InitialPoints"->initialPoints},
	            	MaxIterations->300,PrecisionGoal->2,AccuracyGoal->3
	            	];
	            parameters = fixedOrVariableParameters/. solution;

	       		,(*nekaj je slo narobe*)
			    Print["Something went wrong with fitting at point ", fixedOrVariableParameters[[;;2]]];
			    Print["spectrum=", spectrum];
			    Print["initialPoints=", initialPoints];
			    Print["fixedOrVariableParameters=", fixedOrVariableParameters];
			    Print["variableParametersPositions=", variableParametersPositions];
			    Print["variableParameterLimits=", variableParameterLimits];
			    
			    parameters = fixedOrVariableParameters;
			    parameters[[variableParametersPositions]] = initialPoint;
			    chi2norm = Infinity;
	        	];,
	   
	   	(*else (performFit==False)*)
            parameters = fixedOrVariableParameters;
			parameters[[variableParametersPositions]] = initialPoint;
			chi2norm = Infinity;     
        ]; 
        
(*Print[{"parameters",parameters}];*)

		Write[$stream,OutputForm[Row[Riffle[Flatten[{parameters,chi2norm}],"\t"]]]];
        ClearSystemCache[];
    	(*Return[{parameters,chi2norm}];*)
    ];


(*
CompiledSpectrumBleachFit: runs the simulation for one spectrum to get all parameter values; 
						compiled version returns a joined list of all parameters and chi2;
						bleaching in [frames^-1] 

input:	spectrum						.. list of experimental spectral intensities, ordered as lambdaSample
		initialPoint					.. starting values for variable parameters
		fixedOrVariableParameters		.. list of all parameters; fixed with numbers, variable with symbols, e.g. {1,1,1,1,1,lm1,75,.24,b1,0,0,550,75,.24,0}
		variableParamtersPositions 		.. indices of variable parameters in the list of all parameters, e.g. {6,9} for variable {lm1,b1}
		variableParametersLimits		.. list of limits for variable parameters {{par1min,par1max},..}, e.g. {{400,700},{0,1}}
		variableParameterRangeSizes		.. list of spans between the limits for variable parameters {par1max-par1min,..}, e.g. {300,1}
		performFit						.. fit switch; True if optimization is to be performed, False if not
		$$Fitter					.. compiled minimization function, passed as GLOBAL VARIABLE, containing all the info about lambdaSample, noise, ..	
		
output: list of all parameter values with chi2 appended as last 
*)

ClearAll[CompiledSpectrumBleachFit];
CompiledSpectrumBleachFit[spectrum_, initialPoint_, fixedOrVariableParameters_, 
	variableParametersPositions_, variableParameterLimits_, variableParameterRangeSizes_, 
	performFit_:True] :=
    Module[ {initialPoints, result, out},
        
If[performFit, (*then*)
	        (*minimization of n parameters needs n+1 starting points;
	        one is the initialPoint, others are generated by varying only one parameter for 10% of the range (one along each parameter, +/- random); 
	        clipped to constraints*)
	        initialPoints = 
	        	Join[
	        		{initialPoint},
					Table[
						MapThread[
							Clip,
								{initialPoint + 0.1 variableParameterRangeSizes[[i]] Sign[RandomReal[{-1.,1.}]] UnitVector[Length[variableParametersPositions],i]
								,variableParameterLimits}],
					{i,Length[variableParametersPositions]}]
					];
			initialPoints=Clip[initialPoints,{10^-6,Infinity}];

(*Print[{"initialPoints",initialPoints}];
Print[{"spectrum",spectrum}];*)

	        (* results in the form {chi2, param1, param2, ..} *)
	        Check[
	        	result = $$Fitter[spectrum, initialPoints];(*$$Fitter is a GLOBAL VARIABLE!*);
	       		,(*if any message was generated during fitting, report:*)

			    Print["Something went wrong with fitting at point ", fixedOrVariableParameters[[;;2]]];
			    Print["spectrum=", spectrum];
			    Print["initialPoints=", initialPoints];
			    Print["fixedOrVariableParameters=", fixedOrVariableParameters];
			    Print["variableParametersPositions=", variableParametersPositions];
			    Print["variableParameterLimits=", variableParameterLimits];
			    result={Infinity}~Join~initialPoint;
	        	];,
	   
	   	(*else (performFit==False)*)
           result={Infinity}~Join~initialPoint;     
        ]; 
        
		(*result[[2;;]]=MapThread[Clip[#1,#2]&,{result[[2;;]],variableParameterLimits}];*)

        (*return complete vector of parameters with chi2 added in the end*) 
        out = fixedOrVariableParameters;
        out[[variableParametersPositions]] = Drop[result, 1];
        AppendTo[out, First@result];

(*Print[out];*)

        out
    ];



(*
ReconstructSpectrum: calculates optimal intensity of the fitted spectrum and returns appropriately scaled spectrum with both components;
					intensity determined by dchi2/dI = 0 ... -> I = Sum_i (ySim_i yExp_i) / Sum_i(ySim_i) 		
					(derivation in Iztok's diary IU@FMS_01 p .126)
					
input:	parameters		.. list of all parameters; intensity not yet optimized; bleaching in [frames^-1] 
		spectrumExp		.. list of experimental spectral intensities, ordered as lambdaSample
		lambdaSample	.. list of wavelengths, ordered BY TIME OF ACQUISITION (shuffled in the case of stochastic sampling)				
					
output:	{{spec,specComponent1,specComponent2},intensity}
			spec			.. list of simulated spectral intensities, ordered as lambdaSample
			specComponent1	.. list of simulated spectral intensities of the first spectral component, ordered as lambdaSample
			specComponent2	.. list of simulated spectral intensities of the second spectral component, ordered as lambdaSample
			intensity		.. parameter to scale model lineshape -> replaces the value at 3rd position in the complete list
*)
ClearAll[ReconstructSpectrum];
ReconstructSpectrum[parameters_,spectrumExp_,lambdaSample_] :=
    Block[ {spec,intensity,specComponent1,specComponent2,totalSpec},
        spec = SpectraSum[parameters,lambdaSample];
        totalSpec=Total[spec^2];
        intensity =If[totalSpec>0,(*then*)
        	         Total[spec*spectrumExp]/totalSpec
        	         ,(*else*)
        	         0];
        (*spec = spec*intensity;*)
        spec = SpectraSum[ReplacePart[parameters,3->intensity],lambdaSample];
        specComponent1 = SpectraSum[ReplacePart[parameters,{3->intensity parameters[[5]], 5->1.,11->0.}],lambdaSample];
        specComponent2 = SpectraSum[ReplacePart[parameters,{3->intensity parameters[[11]],5->0.,11->1.}],lambdaSample];
        Return[{{spec,specComponent1,specComponent2},intensity}];
    ];
(*DistributeDefinitions[ReconstructSpectrum];*)

(*compiled version of ReconstructSpectrum*)
ClearAll[ReconstructSpectrumC];
ReconstructSpectrumC =  
    Compile[{{par, _Real, 1}, {spectrumExp, _Real, 1}, {lambdaSample, _Real, 1}},
      Block[ {spec={0.}, intensity=0., specComponent1={0.}, specComponent2={0.}, compPar={0.}, totalSpec=0.},

(*Print[{"par",par}];*)
            spec = SpectraSumC[par, lambdaSample];
(*Print["mu"];*)
            totalSpec=Total[spec^2];   
            intensity =If[totalSpec>0.,(*then*) 
                         Total[spec*spectrumExp]/totalSpec
                         ,(*else*) 
                         0.];
            spec = spec*intensity;
            (*spec = SpectraSum[ReplacePart[par,3->intensity],lambdaSample];*) 
            
            If[par[[5]]<1,
            	compPar=par;
            	compPar[[{3,5,11}]]={intensity*par[[5]],1.,0.};
            	specComponent1 = SpectraSumC[compPar,lambdaSample];
            	compPar[[{3,5,11}]]={intensity*par[[11]],0.,1.};
            	specComponent2 = SpectraSumC[compPar,lambdaSample];
            	,
            	specComponent1 = ConstantArray[0,Length[spectrumExp]];
            	specComponent2 = ConstantArray[0,Length[spectrumExp]];
            ];
            Join[spec,specComponent1,specComponent2,{intensity}]
           ]
           ,{{SpectraSumC[_,_], _Real, 1}}              
           ,CompilationOptions -> {"InlineCompiledFunctions" -> True, "InlineExternalDefinitions" -> True}
]; 

(*
ReconstructSpectrumCWrapper: calls ReconstructSpectrumC and reformats its output which comes as a flat list (specspecComponent1specComponent2intensity);
							
input: 	spectraBlock	.. {parameters,spectrumExp,performFit}
		 parameters		.. list of all parameters; intensity not yet optimized; bleaching in units of [frames^-1]
		 spectrumExp	.. list of experimental spectral intensities, ordered as lambdaSample
		 performFit		.. fit switch; True if optimization was be performed, False if not

output:	{{spec,specComponent1,specComponent2},intensity}
			spec			.. list of simulated spectral intensities, ordered as lambdaSample
			specComponent1	.. list of simulated spectral intensities of the first spectral component, ordered as lambdaSample
			specComponent2	.. list of simulated spectral intensities of the second spectral component, ordered as lambdaSample
			intensity		.. parameter to scale model lineshape -> replaces the value at 3rd position in the complete list
*)      
ReconstructSpectrumCWrapper[spectraBlock_,lambdaSample_]:=
    With[{spectrum=spectraBlock[[2]],parameters=spectraBlock[[1]],performFit=spectraBlock[[3]]},
    Block[{rawReconstruc},  
      If[!performFit, 
      	(*if the spectrum was not optimized, return zeros*) 
      	 {ConstantArray[0,{3,Length@spectrum}],0.}
       ,(*else*)
         (*vrne tri spektre in na koncu intenziteto v enem flat listu*)
         rawReconstruc=ReconstructSpectrumC[parameters,spectrum,lambdaSample];
         {Partition[rawReconstruc,Length@spectrum],Last@rawReconstruc} 
        ]
    ]];
  



(* ::Section:: *)
(*(* MAIN FITTING ROUTINE *)*)
(**)


(*
AnalizeSpectra: take a list of spectra from all points, prepare everything for fitting, run optimization, and draw some fits from randomly chosen points;
					input bleaching rates in [frames^-1], output in [s^-1]
				
input: 	listOfSpectra			.. list of spectra from all points in the image, ordered by TIME OF ACQUISITION
		lambdaSample			.. list of wavelengths, ordered BY TIME OF ACQUISITION (shuffled in the case of stochastic sampling)
		lambdaOrder				.. list of numbers to order lambdaSample by wavelength: lambdaSample[[lambdaOrder]] = Sort[lambdaSample]
		noise					.. standard deviation of background signal
		timeExposure			.. exposure time of each image in seconds
		parameterSwitches		.. list of True/False for all parameters; True for variable ones (to be optimized), False for fixed		
		parameterFixedValues	.. list of values for all parameters; used as initial points for optimization of variable parameters, kept fixed for the rest;  bleachRates in [frames^-1]
		parameterLimits			.. list of limits to keep parameters within during optimization: {{par1min,par1max},..} for all 15 parameters;  bleachRates in [frames^-1]
		bleachRefMode			.. True for experiment with a reference wavelength for bleaching correction (used with linear sampling)
		bleachRefLambda			.. wavelength at which bleaching reference images were recorded [nm]
		bleachRefFrameNumbers	.. consecutive numbers of bleaching reference images in the whole series, e.g. {1,4,7,11,16,23} (not necessarily equally spaced)  
		"UseCompiled"			.. True/False whether to use compiled functions 
		"UseParallel"			.. True/False whether to run optimization on parallel kernels
		"GeneratePlots"			.. True/False whether to generate graphical output

output:	parameters				.. list of all parameters+chi2 for all points (dimensions Npoints x 16); bleaching rates in [s^-1]!
		plots					.. plots of experimental and fitted spectra from 6 points with different intensities
		selectedPointsToLabel	.. list of {x,y} coordinates where the 6 spectra were taken from
		intMaxAbs				.. max intensity of all data points (cut-off at 0.1%)
		{intMinRel,intMaxRel}	.. limits for intensity coding, relative to max intensity (filters 0.5% of highest and lowest values);
*)
ClearAll[AnalizeSpectra];
Options[AnalizeSpectra] = {"UseCompiled" -> True, "UseParallel" -> True, "GeneratePlots"->True};
AnalizeSpectra[listOfSpectra_,lambdaSample_,lambdaOrder_,noise_,timeExposure_,
	parameterSwitches_,parameterFixedValues_,parameterLimits_,
	{bleachRefMode_,bleachRefLambda_,bleachRefFrameNumbers_},
	opts:OptionsPattern[]] :=
    Block[ {t,spectra,
    parameterSymbols,variableParameterSymbols,variableParameterPositions,fixedOrVariableParameters,
    variableParameterRangeSizes,variableParameterStartingPoints,variableParameterConstraints,variableParameterLimits,
    bleachRates,bleachPlacesToReplace,
    kernels,tmpOutFilenames,parameters,intensities,spectraFit,chi2OutsideLimits=100.,
    plots={},selectedPointsToLabel,intFitted,intMaxAbs=1.,intMinRel=0.,intMaxRel=1.,histograms,imagesCC},
		
		
		(************* PREPARE *************)

		spectra = listOfSpectra;

		(*positions of parameters to fit*)
        variableParameterPositions = Flatten[Position[parameterSwitches,True]];
		(*symbols*)
		parameterSymbols = {x,y,i0,i1,p1,lm1,w1,a1,b1,i2,p2,lm2,w2,a2,b2,hemOxy,chi2};
        variableParameterSymbols = parameterSymbols[[variableParameterPositions]];
 
        (*get a vector of all parameters with variable ones replaced by symbols, e.g. {1,1,1,1,1,lm1,w1,.21,b1,0,0,550,60,.2,0.,1.}*)
        fixedOrVariableParameters = parameterFixedValues;
        fixedOrVariableParameters[[variableParameterPositions]] = variableParameterSymbols;
        
        (*if fitted with 2 spectral components with variable portions, replace value for portion2 with 1-portion1*)
        If[ parameterSwitches[[5]]||(parameterFixedValues[[5]]<1),
            fixedOrVariableParameters[[11]] = 1.-p1;
        ];
        fixedOrVariableParameters = ConstantArray[fixedOrVariableParameters,Length[spectra]];

	
        (***** FIT LINEAR BLEACHING *****)
		
		If[bleachRefMode && Or@@parameterSwitches[[{9,15}]],
			t = AbsoluteTiming[
			(*fit bleaching series*)
			bleachRates = BleachFitLinearRef[spectra[[All,5]],bleachRefFrameNumbers];
			bleachRates = Clip[bleachRates,parameterLimits[[9]]];
			(*replace in fixedOrVariableParameters*)
			bleachPlacesToReplace = Select[{9,15},MemberQ[variableParameterPositions,#]&];
			fixedOrVariableParameters[[All,bleachPlacesToReplace]] = Transpose@ConstantArray[bleachRates,Length[bleachPlacesToReplace]];
			(*drop bleach from fitting*)
			variableParameterPositions = Select[variableParameterPositions,!MemberQ[bleachPlacesToReplace,#]&];
			variableParameterSymbols = parameterSymbols[[variableParameterPositions]];
			][[1]];
			PrintTemporary[Row[{"BleachFitLinearRef .. ",t, " s"}]];
		];

 
		(************* PREPARE *************)
       
        (*define ranges and limits for fitting*)
        variableParameterLimits =  parameterLimits[[variableParameterPositions]];
        variableParameterConstraints = (parameterLimits[[#,1]]<=parameterSymbols[[#1]]<=parameterLimits[[#,2]]&)/@variableParameterPositions;
        variableParameterRangeSizes = Flatten[Differences/@variableParameterLimits];


        (*initial values for variable parameters*)
        variableParameterStartingPoints = ConstantArray[parameterFixedValues[[variableParameterPositions]],Length[spectra]];
   
     
       (* If[ optimizeStartingPoint,
            t = AbsoluteTiming[
            variableParameterStartingPoints = OptimizeStartingPoints[spectra,noise,lambdaSample,{fixedOrVariableParameters,variableParameterSymbols,variableParameterStartingPoints[[1]],variableParameterConstraints,variableParameterRangeSizes,parameterConstraints[[variableParameters]]}][[All,variableParameters]];
            ][[1]];
            PrintTemporary[Row[{"OptimizeStartingPoints .. ",t, " s"}]];
        ];*)



		(*add new info to spectra multilist: {{x,y,...},{spectrum},signalToNoise,performFit,{bleachSeries},{variableParameterStartingPoint}}*)
		spectra[[All,1,3;;]] = fixedOrVariableParameters[[All,3;;]];
		spectra[[All,6]] = variableParameterStartingPoints;
		Puts["spectra=",Short@spectra];		
		
		
		(************* MODEL FIT *************)
		
        If[ OptionValue@"UseCompiled", (*then*)

        (**************Compiled************)
            (*create compiled fitting function $$Fitter*)
            Quiet[With[{Minimizer = NelderMeadMinimize`Dump`CompiledNelderMead[
            	                              Function[Evaluate@variableParameterSymbols,
            	                              	       Evaluate@CompiledLimitedChi2[fixedOrVariableParameters[[1]],spectr,lambdaSample,noise,
            	                              	                                   variableParameterPositions,variableParameterLimits,chi2OutsideLimits]]
            	                              ,variableParameterSymbols,"ReturnValues"->"AugmentedOptimizedParameters"]},
                      $$Fitter = Compile[{{spectr,_Real,1},{initialParameterValues,_Real,2}},
                      	                 Minimizer[initialParameterValues,0.001,300]
                      	                 ,CompilationOptions->{"InlineCompiledFunctions"->True,"InlineExternalDefinitions"->True,"ExpressionOptimization"->True}
                      	                 ,RuntimeOptions->{"Speed","EvaluateSymbolically"->False}
                      	                 ,CompilationTarget->"MVM"
                      	                ];
                      ];
            ,{CompiledFunction::cfta}];
        

            If[ OptionValue@"UseParallel", (*then*)
            
            (*parallel*)
                DistributeCompiledFunction["$$Fitter"];
                DistributeDefinitions[CompiledSpectrumBleachFit];
                t = AbsoluteTiming[
                	parameters = ParallelMap[
			                		With[ {spectrum = #[[2]], initialPoint = #[[6]], fixOrVarPar = #[[1]], performFit=#[[4]]}, 
			                			CompiledSpectrumBleachFit[spectrum,initialPoint,fixOrVarPar,
			                				variableParameterPositions,variableParameterLimits,variableParameterRangeSizes,performFit]
			                        ]&,spectra]
                ][[1]];
                PrintTemporary[Row[{"Parallel CompiledSpectrumBleachFit .. ",t, " s"}]];
 
 ,
   
 
            (*else not parallel*)
				t = AbsoluteTiming[
                	parameters = Map[
			                		With[ {spectrum = #[[2]], initialPoint = #[[6]], fixOrVarPar = #[[1]], performFit=#[[4]]}, 
			                			CompiledSpectrumBleachFit[spectrum,initialPoint,fixOrVarPar,
			                				variableParameterPositions,variableParameterLimits,variableParameterRangeSizes,performFit]
			                        ]&,spectra]
                ][[1]];
                PrintTemporary[Row[{"CompiledSpectrumBleachFit .. ",t, " s"}]];


            ];
            (*Puts["parameters=",parameters];*),

        (*ELSE*)    
        (**************Uncompiled************)

           (*prepare temp output log files*)
	        kernels = ParallelEvaluate[$KernelID];
	        tmpOutFilenames = ("log-tmp_kernel-"<>ToString[#]<>".txt")&/@kernels;
            With[ {f = NotebookDirectory[],ps = parameterSymbols},
                ParallelEvaluate[
                	SetDirectory[f];
                	$stream = OpenAppend["log-tmp_kernel-"<>ToString[$KernelID]<>".txt",PageWidth->Infinity];
                	Write[$stream,OutputForm[Row[Riffle[ps,"\t"]]]];
                ]
            ];

	        SetSharedVariable[spectra,lambdaSample,noise,fixedOrVariableParameters,variableParameterSymbols,variableParameterStartingPoints,variableParameterConstraints,variableParameterRangeSizes];
	        t = AbsoluteTiming[
		        	ParallelMap[
			        	With[ {spectrum = #[[2]], initialPoint = #[[6]], fixOrVarPar = #[[1]], performFit=#[[4]]}, 
		               		SpectrumBleachFit[spectrum, lambdaSample, noise, 
		               			initialPoint, fixOrVarPar, variableParameterPositions, variableParameterSymbols, 
								variableParameterLimits, variableParameterRangeSizes, variableParameterConstraints,
								performFit]
			        	]&,spectra]
	        ][[1]];
	        PrintTemporary[Row[{"SpectrumBleachFit .. ",t, " s"}]];
	        ParallelEvaluate[Close[$stream]];
            Clear[$stream];


            (*read the temp log files*)
            parameters = Flatten[Import[#,"Table"][[2;;]]&/@tmpOutFilenames,1];
            parameters = SortBy[parameters,{#[[2]]&,#[[1]]&}];
            DeleteFile/@tmpOutFilenames;
        ];
        
         
        (*correct the portions of the second spectral component and calculate both partial intensities*)
        parameters[[All,11]] = If[ parameterSwitches[[5]]||(parameterFixedValues[[5]]<1),
        	1.-parameters[[All,5]],
        	0
        ];
  
     
(*If[Length[parameters[[1]]]==18, parameters={Drop[parameters[[1]],{17}]},None];*)
	
		(****** RECONSTRUCT SPECTRA ******)

        If[ OptionValue@"UseCompiled", (*then*)
	        (*mySpectra[[All,3]]=parameters;*)
	        spectra[[All,1]]=parameters;

            If[ OptionValue@"UseParallel", (*then*)
		        DistributeCompiledFunction["ReconstructSpectrumC"];
		        t = AbsoluteTiming[	        
			        With[{lambdaSamp=lambdaSample},
			        {spectraFit,intensities} = Transpose@ParallelMap[
			            ReconstructSpectrumCWrapper[#,lambdaSamp]
			          &,spectra[[All,{1,2,4}]]];
			        ]  
	            ][[1]];

      	       PrintTemporary[Row[{"Parallel ReconstructSpectrum Compiled .. ",t, " s"}]];	
           ,(*else not parallel*)
		        t = AbsoluteTiming[	        
			        With[{lambdaSamp=lambdaSample},
			        {spectraFit,intensities} = Transpose@Map[
			            ReconstructSpectrumCWrapper[#,lambdaSamp]
			          &,spectra[[All,{1,2,4}]]];
			        ]  
	            ][[1]];
	            PrintTemporary[Row[{"ReconstructSpectrum Compiled .. ",t, " s"}]];	

          
            ];
        ,(*else not compiled*)
        
	        t = AbsoluteTiming[
		        {spectraFit,intensities} = Transpose[(*Parallel*)Map[
		        If[ spectra[[#,4]],
		            ReconstructSpectrum[parameters[[#]],spectra[[#,2]],lambdaSample],
		            {ConstantArray[0,{3,Length[lambdaSample]}],0.}
		        ]&,Range[Length[spectra]]]];
	        ][[1]];
	        PrintTemporary[Row[{"ReconstructSpectrum .. ",t, " s"}]];
        
        ];

        (*add intensity info to fit parameters and calculate partial intensities of spectral components*)

		parameters[[All,3]] = Clip[intensities,{0,Infinity},{0,Infinity}];
        parameters[[All,{4,10}]] = intensities*parameters[[All,{5,11}]];
	    


        (************* PRESENT FITTING RESULTS *************)
		(*If[Length[parameters]==18,parameters=Drop[parameters,{17}],None];*)


		(*clip parameters to their ranges and round the numbers*)
		parameters[[All,variableParameterPositions]] = Transpose@MapThread[Clip[#1,#2]&,{Transpose@parameters[[All,variableParameterPositions]],parameterLimits[[variableParameterPositions]]}];
		parameters = Transpose@MapThread[Round[#1,#2]&,{Transpose@parameters,{1,1,1,1,.01,.1,.1,.001,.0001(*(timeExposure Length[lambdaSample])^-1/100.*),1,.01,.1,.1,.001,.0001, .01(*(timeExposure Length[lambdaSample])^-1/100.*),.1}}];
	    

(*convert bleching rate units from frames^-1 into s^-1*)
		parameters[[All,{9,15}]]/=(timeExposure+.2);
	     		
		(*in case of fitting linear bleaching data, add back the eliminated variableParameterPositions for bleaching rates*)
		If[bleachRefMode && Or@@parameterSwitches[[{9,15}]],
			variableParameterPositions = Sort[Join[variableParameterPositions,bleachPlacesToReplace]];
		];
		(*if portions were variable, draw plots for both*)
		If[parameterSwitches[[5]],
			variableParameterPositions = Sort[Join[variableParameterPositions,{4,10,11}]];,
			variableParameterPositions = Sort[Join[variableParameterPositions,{4}]];
		];


        If[ OptionValue@"GeneratePlots", (*then*)
	       	
	       	t = AbsoluteTiming[
	       		
	       		(*get useful intensity ranges; 0.5% of extreme values (usually bad fits with high or low - even negative - intensity) ignored*)
	    		intFitted = Sort[Select[parameters[[All,3]],0<#(*<2^16*)&]];
	    		intMaxAbs = intFitted[[Round[.999 Length[intFitted]]]];
	    		intMinRel = Round[intFitted[[Round[.005 Length[intFitted]]+1]]/intMaxAbs,.01];
	    		intMaxRel = Round[intFitted[[Round[.995 Length[intFitted]]]]/intMaxAbs,.01];


					       		
	        	(*generate plots of fitted spectra from selected points - evenly distributed over intensities*)
	        	{plots,selectedPointsToLabel} = PlotFittedSpectra[parameters,spectra[[All,2]],spectraFit,lambdaSample,lambdaOrder,
	        		{bleachRefMode,(*bleachSeries*)spectra[[All,5]],bleachRefFrameNumbers,
	        		(*position of fitted b*)If[Length@bleachPlacesToReplace>0,bleachPlacesToReplace[[1]],9],timeExposure},
	        		intMaxAbs {intMinRel,intMaxRel}];
	        ][[1]];

	        PrintTemporary[Row[{"Generating spectra plots .. ",t, " s"}]];
	        
	   ];
		(*export also spectraFit, which are data of fit, 1.spec Comp, 2. spec Comp*)
       Return[{parameters,plots,selectedPointsToLabel,{Round[intMaxAbs],intMinRel,intMaxRel},spectraFit}];
    ];
 


(* ::Section:: *)
(*(* IMAGE ANALYSIS ROUTINE *)*)
(**)


(*
AnalizeImages: prepares images (crop, alignment), extracts spectra and sends them to optimization;
			input bleaching rates in [frames^-1], output in [s^-1]

input:	images					.. 3D matrix with 16-bit data from all images in the spectral series; sorted by wavelength, as saved
		cropXY					.. {{xmin,ymin},{xmax,ymax}} of the region to be analized [pix]
		lambdaSample			.. list of wavelengths as sampled in the experiment
		timeExposure			.. exposure time of each image in seconds
		avg						.. how many pixels to average in each direction; avg x avg
		noise					.. standard deviation of background signal
		signalThreshold			.. intensity threshold to filter points with signal < signalThreshold * noise
		parameterSwitches		.. list of True/False for all parameters; True for variable ones (those to be optimized), False for fixed
		parameterFixedValues	.. list of values for all parameters; used as initial points for optimization of variable parameters, kept fixed for the rest;  bleachRates in [frames^-1]
		parameterLimits			.. list of limits to keep parameters within during optimization: {{par1min,par1max},..} for all 15 parameters;  bleachRates in [frames^-1]
		bleachRefMode			.. True for experiment with a reference wavelength for bleaching correction (used with linear sampling)
		bleachRefImages			.. 3D matrix with 16-bit data from images in the bleaching reference series; sorted by time taken
		bleachRefLambda			.. wavelength at which bleaching reference images were recorded
		bleachRefFrameNumbers	.. consecutive numbers of bleaching reference images in the whole series, e.g. {1,4,7,11,16,23} (not always equally spaced)
		"UseCompiled"			.. option: True/False whether to use compiled functions 
		"UseParallel"			.. option: True/False whether to run optimization on parallel kernels
		"GeneratePlots"			.. option: True/False whether to generate graphical output
		"AlignImages"			.. option: True/False whether to align images to compensate for particle movement before extracting spectra
		"AlignMethod"			.. option: built-in numerical method for image alignment: "Translation", "Rigid" (translation + rotation), ... 

output:	parameters				.. list of all parameters+chi2 for all points (dimensions Npoints x 16); bleaching rates in [s^-1]!
		plots					.. plots of experimental and fitted spectra from 6 points with different intensities
		selectedPointsToLabel	.. list of {x,y} coordinates where the 6 spectra were taken from
		intMaxAbs				.. max intensity of all data points (cut-off at 0.1%)
		{intMinRel,intMaxRel}	.. limits for intensity coding, relative to max intensity (filters 0.5% of highest and lowest values)
		imageIntensityCompareAlign	.. aligned intensity image (*column of two small intensity images to compare before/after alignment*)
		spectra[[All,2]]              ..spectral points for each image pixel
*)
ClearAll[AnalizeImages];
Options[AnalizeImages] = {"UseCompiled" -> True, "UseParallel" -> True, "GeneratePlots"->True, "AlignImages"->False, "AlignMethod"->"Translation"};
AnalizeImages[images_,cropXY_,lambdaSample_,timeExposure_,avg_,noise_, signalThreshold_,
	          parameterSwitches_,parameterFixedValues_,parameterLimits_, 
	          {bleachRefMode_, bleachRefImages_, bleachRefLambda_, bleachRefFrameNumbers_},
	          opts:OptionsPattern[]] :=
    Block[{t,lambdaOrder,lambdaShuffle,imagesCrop,bleachRefImagesCrop,imageIntensityCompareAlign,spectra,spectraNoise,signal,nAvg,(*noise,snr,*)imagesAvg,
    	referenceSpectraIndex,nSpectra,colors,parameters,imagesCC,histograms,plots,selectedPointsToLabel,intMaxAbs,intMinRel,intMaxRel,specFit},
        
        ParallelEvaluate[$HistoryLength=0];

        lambdaOrder = Ordering[lambdaSample];  (*list of indices to sort stochastic lambdas*)
        lambdaShuffle = Flatten[Position[Sort[lambdaSample],#]&/@lambdaSample];  (*list of indices to shuffle ordered lambdas into original sampling list*)
        imagesCrop = images[[All,512-cropXY[[2,2]]+1;;512-cropXY[[1,2]]+1,cropXY[[1,1]];;cropXY[[2,1]]]];
		
		If[bleachRefMode,
			bleachRefImagesCrop = bleachRefImages[[All,512-cropXY[[2,2]]+1;;512-cropXY[[1,2]]+1,cropXY[[1,1]];;cropXY[[2,1]]]];
		];
		
		nAvg = Round[(avg-1)/2]; (*how many pixels to average in each direction; avg x avg -> (2 nAvg + 1) x (2 nAvg + 1)*)

		If[OptionValue@"AlignImages",
			t=AbsoluteTiming[
				{imagesCrop,imageIntensityCompareAlign} = AlignImages[imagesCrop,nAvg,"AlignMethod"->OptionValue@"AlignMethod"];
				If[bleachRefMode, 
					bleachRefImagesCrop = AlignImages[bleachRefImagesCrop,nAvg,"AlignMethod"->OptionValue@"AlignMethod"][[1]];
				];
			][[1]];
			PrintTemporary[Row[{"Aligning images .. ",t, " s"}]];
			,
			imageIntensityCompareAlign="";
		];

		(*get spectra from images in the form {{x,y}, spectrum, SNR, fitSwitch, bleachSeries (for bleachRefMode), place for variableParameterStartingPoint}*)
        spectra = SpectraFromImages[imagesCrop,cropXY,lambdaShuffle,nAvg,noise,signalThreshold,{bleachRefMode,bleachRefImagesCrop},Count[parameterSwitches,True]];
        PrintTemporary[Row[{"N of points above threshold: ",Count[spectra[[All,4]],True]}]];

        (*fit all selected spectra*)
        {parameters,plots,selectedPointsToLabel,{intMaxAbs,intMinRel,intMaxRel},specFit} = 
            AnalizeSpectra[spectra,lambdaSample,lambdaOrder,noise,timeExposure,
            		parameterSwitches,parameterFixedValues,parameterLimits,
            		{bleachRefMode,bleachRefLambda,bleachRefFrameNumbers}
            	,"UseCompiled"->OptionValue@"UseCompiled"
            	,"UseParallel"->OptionValue@"UseParallel"
            	,"GeneratePlots"->OptionValue@"GeneratePlots"];

        (*Export[filenameCrop<>"_I-image.png",imageIntensity];*)
        (*second last parameter to export are all spectral points from each image pixel*)
        (*last parameter to export fit of the first component*)
        Return[{parameters,plots,selectedPointsToLabel,intMaxAbs(*{intMaxAbs,intMinRel,intMaxRel}*),imageIntensityCompareAlign,spectra[[All,2]],specFit[[All,2]]}];
    ];


(* ::Section:: *)
(*(* PRESENTATION OF RESULTS *)*)


(*
IntensityImage: creates an intensity image by superimposing all spectral images

input:  intensityData			.. data of all spectral images superimposed, normalized to the maximal value
		adjustContrastMinMax     .. list of two values between 0 and 1 to scale contrast for intensityImage: {minIntensity,maxInensity} 

output:  imageIntensity         	.. ArrayPlot of all spectral images superimposed; intensity scaled from min to max
*)
ClearAll[IntensityImage]
IntensityImage[intensityData_,adjustContrastMinMax_]:= 
	ArrayPlot[intensityData, PlotRange -> {{0, 512}, {0, 512}, adjustContrastMinMax}, ColorFunction->GrayLevel, ClippingStyle->{Black,White}, 
		Frame->False, ImageSize->{300}];


(*
IntensityHistogram: plots a distribution of pixel brightness values

input:	intensityData			.. data of all spectral images superimposed, normalized to the maximal value
		adjustContrastMinMax	.. list of two values between 0 and 1 to scale contrast for intensityImage: {minIntensity,maxInensity} 

output: intensityHistogram		.. Histogram of brightness values across the image
*)
ClearAll[IntensityHistogram]
IntensityHistogram[intensityData_,adjustContrastMinMax_]:=
	Histogram[Flatten[intensityData], {(adjustContrastMinMax[[2]] - adjustContrastMinMax[[1]])/100.}, {"Log", "Probability"},
		PlotRange -> {{.5,1.5} adjustContrastMinMax, {-13.5,0}}, PlotRangePadding -> None, 
		Axes -> None, ChartStyle -> Black, 
		Frame -> True, FrameLabel -> {"I [au]", "dP/dI"}, LabelStyle -> {FontSize -> 14}, LabelingFunction -> None, FrameTicks -> {{None,None},{Automatic,Automatic}},
		PerformanceGoal -> "Speed",	ImageSize -> 300, ImagePadding->{{60,5},{60,5}}];


(*
PlotFittedSpectra: generates plots of fitted spectra from 6 points that are approx. evenly distributed over intensity range

input:	parameters				.. list of all parameters+chi2 for all points (dimensions Npoints x 14); bleaching rates in [s^-1]
		spectraExp				.. list of experimental spectra from all points, ordered by TIME OF ACQUISITION 
		spectraFit				.. list of simulated spectra from all points, ordered by TIME OF ACQUISITION
		lambdaSample			.. list of wavelengths as sampled in the experiment
		lambdaOrder				.. list of numbers to order lambdaSample by wavelength: lambdaSample[[lambdaOrder]] = Sort[lambdaSample]
		bleachRefMode			.. True for experiment with a reference wavelength for bleaching correction (used with linear sampling)
		bleachSeries			.. list of intensities through bleachRefImages from all points
		bleachRefFrameNumbers	.. consecutive numbers of bleaching reference images in the whole series, e.g. {1,4,7,11,16,23} (not always equally spaced)
		fittedBleachPosition	.. position of the one bleaching parameter that was optimized (9 or 15); 8 if both fixed
		timeExposure			.. exposure time of each image in seconds
		{intMin,intMax}			.. limits for intensity
		
output: plots					.. plots of experimental and fitted spectra from 6 points with different intensities
		selectedPointsToLabel	.. list of {x,y} coordinates where the 6 spectra were taken from		
*)
ClearAll[PlotFittedSpectra]
PlotFittedSpectra[parameters_,spectraExp_,spectraFit_,lambdaSample_,lambdaOrder_,
	{bleachRefMode_,bleachSeries_,bleachRefFrameNumbers_,fittedBleachPosition_,timeExposure_},{intMin_,intMax_}] :=
    Block[ {selectedParameters,intMaxAbs=1.,(*intMax,intMin,*)selectedIntensities,selection,selectedPointsToLabel,parametersNoBleach,spectraFitNoBleachSelected,labels,tooltips,plots,bleachPlots},
    	(*take only fitted ones*)
        selectedParameters = Select[parameters,0<#[[3]](*<2^16*)&];
        (*take 6 intensities evenly distributed over the range*)
		selectedIntensities = Range[intMax,intMin+(intMax-intMin)/6.,-(intMax-intMin)/6.];
        (*find indexes of points in the original dataset that are close enough to these intensities*)
        selection = Flatten[
            ParallelTable[With[{int=selectedIntensities[[i]]},
                Position[
                	parameters[[All,3]],
                    (*Select[selectedParameters[[All,3]],Abs[1-#/selectedIntensities[[i]]]<.05&,1][[1]],*)
                    Select[RandomSample[selectedParameters[[All,3]]],Abs[1-#/int]<.05&,1][[1]],
                    {1},1]],
                {i,Length[selectedIntensities]}]
                ];
		selectedPointsToLabel = parameters[[selection,;;2]];
		(*spectra without bleaching*)
		parametersNoBleach = parameters[[selection,;;-2]];


		parametersNoBleach[[All,{9,15}]] = 0.;
		spectraFitNoBleachSelected = ParallelMap[
			With[{par=#},SpectraSumC[par,lambdaSample]]&,
			parametersNoBleach];
		(*plot labels with {x,y} info and mouse-over tooltips with parameter values*)
		labels = Style["point "<>ToString[Position[selection,#][[1,1]]]<>": {x,y} = "<>ToString[parameters[[#,;;2]]],9]&/@selection;
		tooltips = {$parameterSymbolsForm,parameters[[#]]} &/@selection;
		If[DeleteDuplicates[parameters[[All,5]]]=={1.},
			tooltips[[All,2,10;;17]]="-"
		];


		(*spectra plots*)			
		plots =(*Grid[Partition*) Row[
            ParallelMap[
            	With[{specExp=spectraExp[[selection[[#]],lambdaOrder]],label=labels[[#]],tooltip=tooltips[[#]],
	            	specFit=spectraFit[[selection[[#]],All,lambdaOrder]],specFitNoBleach=spectraFitNoBleachSelected[[#,lambdaOrder]]},
	            	Tooltip[
		            	Show[
		                	ListPlot[specExp,DataRange->{Min[lambdaSample],Max[lambdaSample]},PlotRange->{-0.5,1.1 intMax},
		                		Frame->True,FrameLabel->{"\[Lambda] [nm]",None},LabelStyle->{FontSize->9},
		                		PlotLabel->label,ImageSize->180],
		                	ListLinePlot[specFit,DataRange->{Min[lambdaSample],Max[lambdaSample]},PlotRange->{-0.5,1.1 intMax},
		                		PlotStyle->(Directive[Darker[Blue],#1]&)/@{Dashing[{}],Dotted,Dotted}],
		                	ListLinePlot[specFitNoBleach,DataRange->{Min[lambdaSample],Max[lambdaSample]},PlotRange->{-0.5,1.1 intMax},
		                		PlotStyle->Directive[Gray,Dashed]]
		                	],
		                TableForm[tooltip]
		                ]
	            	]&,
            	Range[Length[selection]]
            	],
            "     "(*5]*)];
        (*bleachSeries plots*)
        If[bleachRefMode,
			bleachPlots = Row[
	           		ParallelMap[
		           		With[{bleachData=bleachSeries[[selection[[#]]]],bleachPar=parameters[[selection[[#]],fittedBleachPosition]]},
			           		Show[
			                	ListPlot[Transpose[{bleachRefFrameNumbers,bleachData}],
			                		PlotRange->{{0,All},{0,1.1 intMax}},PlotStyle->Red,
			                		Frame->True,FrameLabel->{"frame N",None},LabelStyle->{FontSize->9},
			                		PlotLabel->labels[[#]],ImageSize->180],
			                	Plot[bleachData[[1]] Exp[-bleachPar*(timeExposure+.2) (n-1)], 
			                		{n, 1, bleachRefFrameNumbers[[-1]]},
			                		PlotStyle->Red]
		                		]
		           		]&,
	            	Range[Length[selection]]
	           		],
            	"     "];
			plots = Column[{plots,bleachPlots},Spacings->2];     
 
        ];
        Return[{plots, selectedPointsToLabel}]
    ]


(*
WeightedHistogramData: generates points for histogram-like plot where contributions of points are weighted by their intensity (first parameter)

input:	data			.. list of pairs {intensity,parameterValue} for all points
		binSize			.. width of bins to group points into

output: dataWeighted 	.. list of pairs {binPosition,binIntensity} for histogram-like plot
*)
ClearAll[WeightedHistogramData]
WeightedHistogramData[data_(*as {intensity,parValue}*), binSize_] := 
  Block[{dataGather, binIntensities, binPositions,dataWeighted, histogram},
	(*take only fitted ones*)
   dataGather = Select[data,0<#[[1]](*<2^16*)&];								
   (*sort by parameterValue to put in the right order*)
   dataGather = SortBy[dataGather, Last];							
   (*round parameterValue to use later as binPosition*)
   dataGather[[All, 2]] = Round[dataGather[[All, 2]], binSize];		
   (*group by parameterValue*)
   dataGather = GatherBy[dataGather, Last];							
   
   (*center bin positions*)
   binPositions = dataGather[[All, 1, 2]]-binSize/2.;
   
   (*sum up all intensities in each bin*)
   binIntensities = Total /@ dataGather[[All, All, 1]];				
   (*convert into "probability"*)
   binIntensities /= (Total[binIntensities] binSize);
   
   (*combine positions and intensities*)				
   dataWeighted = Transpose[{binPositions, binIntensities}];
   (*add zeros before first and after last point - needed for nicer plots*)
   dataWeighted = Append[Prepend[dataWeighted,{dataWeighted[[1,1]]-binSize,0}],{dataWeighted[[-1,1]]+binSize,0}];
   (*Print[Show[{Histogram[Select[data, #[[1]] > 0 &][[All, 2]], {binSize},"PDF"],ListLinePlot[dataWeighted]}]];*)
	Return[dataWeighted];
   ];
 
 (*
 WeightedHistograms: creates histograms for all variable parameters and chi2; contributions of points are weighted by spectral intensity and portion of the component
 
 input:	data						.. list of all parameters+chi2 from all points; bleaching rates in [s^-1]
 		variableParameterPositions	.. positions of variable parameters in each parameter list, e.g. {6,9} for {lm1,b1} 
 		colorCodeLimits				.. list of limits to scale colors: {{par1min,par1max},..} for all 15 parameters
 
 output: histograms					.. list of histogram-like plots for chi2 and variable parameters
 *)  
ClearAll[WeightedHistograms]
WeightedHistograms[data_,variableParameterPositions_,colorCodeLimits_] :=
    Block[ {data2show,intMaxAbs,intMin,intMax,binSizes,intensities,histogramData,histograms},
    	binSizes = {1,1,1,10,.01,.5,1,.01,.0001,10,.01,.5,1,.01,.0001,.01,.1};
		(*binSizes[[{4,10}]] = Round[(colorCodeLimits[[{4,10},2]]-colorCodeLimits[[{4,10},1]])/20.];*)
		binSizes = (colorCodeLimits[[All,2]]-colorCodeLimits[[All,1]])/25.;
  	
		(*take only fitted ones*)
        data2show = SortBy[Select[data,0<#[[3]](*<2^16*)&],#[[3]]&];
		
		(*clip intensities to the specified ranges and normalize them*)
        intensities = MapThread[(Clip[data2show[[All,#1]],{#2,#3}]-#2)/(#3-#2)&,
						{{4,10},colorCodeLimits[[{4,10},1]],colorCodeLimits[[{4,10},2]]}];
		(*convert data pairs {intensity*portion, parameterValue} into lists of pairs of {binPosition,binIntensity} - for each variable parameter and chi2*)

histogramData = WeightedHistogramData[(*data2show[[All,{3,#}]]*)
        									Transpose[{
        										intensities[[If[#<=9,1,2]]],(*partial intensities*)
        										data2show[[All,#]](*parameter values*)
        										}],										
        									binSizes[[#]]]&/@variableParameterPositions(*Join[variableParameterPositions,{-1}]*);
       (*draw plots*)
        histograms = Join[
            (*(*chi2 histogram*)
            {ListLinePlot[histogramData[[-1]],
            	InterpolationOrder->0,PlotRange->{{0,Automatic},{0,All}},AxesOrigin->{0,0},Filling -> Axis,
            	Frame->True,FrameTicks->{{None,None},{Automatic,Automatic}},FrameLabel->{Style[$parameterSymbolsForm[[-1]]],None
            	(*	Style["dP/d"<>StringCases[$parameterSymbolsForm[[-1]],x__~~" ["~~__->x][[1]]]*)},
            	LabelStyle->{FontSize->16},ImagePadding->{{20,20},{40,10}},ImageSize->300]},*)
            (*histograms of fitted parameters*)
            Table[
            	ListLinePlot[histogramData[[i]],
	            	InterpolationOrder->0,PlotRange->{colorCodeLimits[[variableParameterPositions[[i]]]],{0,All}},AxesOrigin->{0,0},Filling -> Axis,
	            	Frame->True,FrameTicks->{{None,None},{Automatic,Automatic}},
	            	FrameLabel->{Style[$parameterSymbolsForm[[variableParameterPositions[[i]]]]],None
	            		(*Style["dP/d"<>StringCases[$parameterSymbolsForm[[variableParameterPositions[[i]]]],x__~~" ["~~__->x][[1]]]*)},
	            	LabelStyle->{FontSize->16},ImagePadding->{{20,20},{50,10}},ImageSize->300],
            	{i,Length[variableParameterPositions]}]
            ];
        Return[histograms]
    ];


(*
ColorCodedImage: creates a color-coded image according to values of a parameter and corresponding color legend; data as {{x,y,partialIntensity,value},..}

input:	data			.. list of {x,y,partialIntensity,parameterValue} for all points 
		nAvg			.. averaging for nAvg pix in every direction from central point -> over (2 nAvg + 1) x (2 nAvg + 1) pix 
		{parMin,parMax} .. parameter limits for color scaling
		{iMin,iMax}		.. limits for intensity coding (absolute signal values)
		parameterName	.. a string to label legend axis
		hemRegionDensityThres  ..threshold value of Hem spectral portion (above which Hem is 100% present in the tissue). The value defines also the attenuantion of the blood regions in \[Lambda]max coded image.
	    parameterP2       ..relative portion of second spectral component

output:	Row[{colorCodedImage, legend},"  "]
*)
ClearAll[ColorCodedImage]
ColorCodedImage[data_,nAvg_,{parMin_,parMax_},{iMin_,iMax_},parameterName_,hemRegionDensityThresh_,parameterP2_] :=
    Block[ {Color,ColorHSB,int,intMaxAbs,colors,dataXYHSB ,intMin,intMax,dataScaled,intensitiesScaled,imageCC,legendScale,parameterLabel,legend,g,dataPop},
        (*definition of color scheme: parameterValue defines hue, intensity defines brightness, saturation set to max (1.)*)
        Color[{intensity_,parameterValue_}] := Hue[.8(1-parameterValue^1.7(*1.3*))^1.7(*1.5*),1.,intensity]; (*returns Hue[H,S,V], needed for legend*)
        ColorHSB[{intensity_,parameterValue_}] := {.8(1-parameterValue^1.7(*1.3*))^1.7(*1.5*),1.,intensity}; (*only returns {H,S,V} values, needed for Raster*)
        
        (*scale intensities and parameterValues to interval [0-1]*)
		dataPop=ConstantArray[0,Length[data]];

(*Print[{"parMin",parMin}];
Print[{{"iMin","iMax"},{iMin,iMax}}];
*)
        Which[500 < parMin < 600, 				
(*		If[parameterP2[[#]] >= hemRegionDensityThresh, 
			dataPop[[#]] = N[data[[#,3]]*Exp[-(parameterP2[[#]]-hemRegionDensityThresh)/0.01]], 
			dataPop[[#]]=data[[#,3]]]
			&/@ Range[Length[data]];*)
				dataPop = N[data[[#,3]]*Exp[-Power[Abs[parameterP2[[#]]-0.03]/(hemRegionDensityThresh*1.5 + 0.)]]]&/@ Range[Length[data]];
				dataScaled = Transpose[{Abs[(dataPop-iMin)]/(iMax-iMin),((data[[All,-1]]-parMin)/(parMax-parMin))}],
			0.03 < parMin < hemRegionDensityThresh,
				dataPop = N[data[[#,3]]*Exp[-Power[hemRegionDensityThresh/(parameterP2[[#]] + 0.001),2]]]&/@ Range[Length[data]];
				dataScaled = Transpose[{Abs[(dataPop-iMin)]/(iMax-iMin),((data[[All,-1]]-parMin)/(parMax-parMin))}],
			hemRegionDensityThresh < parMin < 0.7,
				dataPop = N[data[[#,3]]*Exp[-Power[hemRegionDensityThresh/(parameterP2[[#]] - 0.001),2]]]&/@ Range[Length[data]];
				dataScaled = Transpose[{Abs[(dataPop-iMin)]/(iMax-iMin),((data[[All,-1]]-parMin)/(parMax-parMin))}],
			Not[500 < parMin < 600] || Not[0.03 < parMin < 0.7],
			    dataScaled = Transpose[{(data[[All,3]]-iMin)/(iMax-iMin),((data[[All,-1]]-parMin)/(parMax-parMin))}];    		
];
			

		dataScaled = Clip[Chop[dataScaled],{0,1}];
        (*define colors*)
        colors = ColorHSB/@dataScaled;
        (*combine {x,y} and {H,S,B} into a list of {{x,y},{H,S,B}}*)
        dataXYHSB = Transpose[{data[[All,;;2]],colors}];
        (*split list by rows*)
        dataXYHSB = GatherBy[dataXYHSB ,#[[1,2]]&]; 

        (*create color-coded image*)

        imageCC = Graphics[{Raster[dataXYHSB [[All,All,-1]],ColorFunction->Hue]},
        	(*Frame->True,(*FrameLabel->{"x","y"},*)FrameTicks->{{All,None},{All,None}},LabelStyle->{FontSize->12},*)ImageSize->{{300},{300}}];
        (*create a color legend*)
		legendScale = Log10[parMax];
		legendScale = If[-3 < legendScale < 3, 1, Floor[legendScale]];
		parameterLabel = If[legendScale==1, parameterName, StringReplace[Evaluate[parameterName]," ["->StringJoin[" [10e",ToString[legendScale],"  "]]];
		legend = With[{min = If[legendScale==1, parMin, parMin/10^legendScale], max = If[legendScale==1, parMax, parMax/10^legendScale]},
			ListDensityPlot[
        		Flatten[Table[{x,y,(y-min)/(max-min)},{x,0,1,.5},{y,min,max,(max-min)/10.}],1],
        		PlotRange->{{0,1},{min,max},{0,1}},ColorFunction->Function[{value},Color[{1.,value}]],ColorFunctionScaling->False,
        		FrameTicks->{{All,None},{None,None}},FrameLabel->{None,TraditionalForm[parameterLabel]},LabelStyle->{FontSize->18},
				AspectRatio->8.,ImageSize->{300},ImagePadding->{{70,1},{10,10}}]
			];
        g = Row[{imageCC,legend},"  "];
        Return[g];
    ];



(*
ColorCodedImagesPointLabel: adds small white numbers to all color-coded images at points where plotted spectra where extracted from

input: 	imagesCC		.. list of elements Row[unlabelled color-coded image, corresponding legend] for all fitted parameters
		pointsToLabel	.. list of {x,y} coordinates where to put labels	
		dataXY			.. list of {x,y} coordinates from all points to get exact lower-left corner from
		nAvg			.. averaging for nAvg pix in every direction from central point -> over (2 nAvg + 1) x (2 nAvg + 1) pix 

output: imagesCClabeled .. list of elements Row[{LABELLED color-coded image, corresponding legend},"  "] for all fitted parameters
*)
ClearAll[ColorCodedImagesPointLabel]
ColorCodedImagesPointLabel[imagesCC_,pointsToLabel_,xy_,nAvg_]:=
	Block[{pointPositions, pointLabels, imagesCClabeled,(*extract lower left corner:*)xMin=Min[xy[[All,1]]],yMin=Min[xy[[All,2]]]},
		(*get positions of selected points where spectra were taken from; in relative coordinates of Raster as half-integers: {{.5,.5},{1.5,.5},...}*)
		pointPositions = ((# - {xMin,yMin})/ (2 nAvg + 1) +.5) &/@pointsToLabel;
        (*create small white numbers as labels*)
        pointLabels = MapThread[												
        	Inset[Style[ToString[#2],White,10],#1,Alignment->Center]&,
        	{pointPositions, Range[Length[pointsToLabel]]
        	}];
        (*label the images*)
        imagesCClabeled = imagesCC[[All,1]]; (*get rid of Row*)
        imagesCClabeled[[All,1]] = Show[#,Graphics[pointLabels]]&/@imagesCClabeled[[All,1]]; (*add labels*)
        imagesCClabeled = Row[#,"  "]&/@imagesCClabeled; (*group back into rows[image,legend]*)
        Return[imagesCClabeled];
	];


(*
ImagesParametersTable: sorts images, resulting from parameter analysis (color-coded images or histograms), into a 2D matrix, 
						aligning the parameters of the two spectral components
						
input:  images						.. a list of created images
		variableParameterPositions	.. positions of variable parameters in each parameter list, e.g. {4,6,9} for {i1,lm1,b1}, including partial intensities  

output: imagesTable					.. 2D matrix of images; rows for spectral components, columns for parameters 
*)
ClearAll[ImagesParametersTable]
ImagesParametersTable[images_,variableParameterPositions_]:=Block[{imagesTable},
	If[MemberQ[variableParameterPositions,5],(*2 spectral components fitted*)
		(*create empty table for all possible parameters*)
		imagesTable = ConstantArray[Null,{2,6}];
		With[{(*define columns in the table for all fitted parameters of each spectral component (1,2)*)
			columns1 = Mod[Select[variableParameterPositions,#<=9&]-4,6]+1, 
			columns2 = Mod[Select[variableParameterPositions,10<=#<=16&]-4,6]+1},
			(*fill the table*)
			imagesTable[[1,columns1]] = images[[;;Length[columns1]]];
			imagesTable[[2,columns2]] = images[[Length[columns1]+1;;]];
			(*delete empty columns*)
			imagesTable = imagesTable[[All,Sort[DeleteDuplicates[Join[columns1,columns2]]]]]
		]
		,
		(*1 spectral component fitted*)
		imagesTable = {images}	
	];
	
	Return[imagesTable]
];


(*
AnalizeFitParameters: draw color-coded images and histograms of fitted parameter values  

input:	parameters				.. list of all parameters+chi2 for all points (dimensions Npoints x 16); bleaching rates in [s^-1]!
		parameterSwitches		.. list of True/False for all parameters; True for variable ones (to be optimized), False for fixed		
		colorCodeLimits			.. list of limits to scale colors: {{par1min,par1max},..} for all 15 parameters
		avg						.. averaging for avg x avg pix in every direction
		selectedPointsToLabel	.. list of {x,y} coordinates where the 6 spectra were taken from
        hemRegionDensityThresh    ..
  		
output:	histograms				.. 2D matrix of histograms for chi2 and all variable parameters (2 rows for 2 spectral components)
		imagesCC				.. 2D matrix of Row[{color-coded image, corresponding legend}] for all fitted parameters
		imagesCCshow			.. imagesCC with labels drawn over the pixels where the plotted spectra were extracted from
*)
ClearAll[AnalizeFitParameters]
AnalizeFitParameters[parameters_,parameterSwitches_,colorCodeLimits_,avg_,selectedPointsToLabel_,hemRegionDensityThresh_]:=
	Block[{nAvg, variableParameterPositions, t, histograms, imagesCC,imagesCCshow},
		
		(*averaging over nAvg poins in each direction from the central point (avg = 2 nAvg + 1)*)	
		nAvg = Round[(avg-1)/2];
		
		(*indexes of variable parameters in the list of parameters*)
		variableParameterPositions = Flatten[Position[parameterSwitches,True]];

		
		(*include intensity data; if 2 components were used, draw plots for both portions and intensities*)
		If[parameterSwitches[[5]],
			variableParameterPositions = Sort[Join[variableParameterPositions,{4,10,11}]];,
			variableParameterPositions = Sort[Join[variableParameterPositions,{4}]];

		];


		(*generate histograms, weighted by intensities*)
        t = AbsoluteTiming[
        	histograms = WeightedHistograms[parameters,variableParameterPositions,colorCodeLimits];
          	(*ordinary histograms (not intensity-weighted):*)
        	(*histograms = Row[Join[
        			{Histogram[Select[parameters,#[[3]]>0&][[All,-1]],PlotRange->{{0,All},{0,All}},AxesOrigin->{0,0},Frame->True,FrameLabel->{"\!\(\*SuperscriptBox[\(\[Chi]\), \(2\)]\)","N"},LabelStyle->{FontSize->12},ImageSize->250]},
        			Table[Histogram[Select[parameters,#[[3]]>0&][[All,i]],PlotRange->{{All,All},{0,All}},AxesOrigin->Automatic,Frame->True,FrameLabel->{Style[$parameterSymbolsForm[[i]]],"N"},LabelStyle->{FontSize->12},ImageSize->250],{i,variableParameters}]
        		],"\t"];*)
    	 ][[1]];
        PrintTemporary[Row[{"Generating histograms .. ",t, " s"}]];


		(*generate color-coded images*)
        t = AbsoluteTiming[
             imagesCC = ParallelMap[
							With[{
								parCC = parameters[[All,{1,2,If[#<=9,4,10],#}]], 
								limits = colorCodeLimits[[#]],
								symbol = $parameterSymbolsForm[[#]]
								},
								ColorCodedImage[parCC,nAvg,limits,colorCodeLimits[[If[#<=9,4,10]]],symbol,hemRegionDensityThresh,parameters[[All,11]]]
							]&,
							variableParameterPositions];
        ][[1]];
        PrintTemporary[Row[{"Generating color-coded images .. ",t, " s"}]];

		
		(*draw labels over color-coded images from where the plotted spectra were extracted from*)
		imagesCCshow = ColorCodedImagesPointLabel[imagesCC, selectedPointsToLabel, parameters[[All, ;; 2]], nAvg];
		
		(*sort the images into a grid according to fitted parameters*)
		imagesCC = ImagesParametersTable[imagesCC, variableParameterPositions];
		imagesCCshow = ImagesParametersTable[imagesCCshow, variableParameterPositions];
		histograms = ImagesParametersTable[histograms, variableParameterPositions];
		

		Return[{histograms, imagesCC, imagesCCshow}];
];




(* ::Section:: *)
(*(* EXPORT *)*)


(*
ExportSpectra: exports spectra from the plot in manipulate and fitted parameters

input:	experimentName			.. string for name basis
		cropBackground			.. {{xmin,ymin},{xmax,ymax}} of the region where background was extracted from [pix]
		avg						.. averaging avg x avg [pix]
		parametersPoint			.. fitted parameters for the fit from the green point
		spectraPoint			.. list of {lambdaSample, experimental spectrum intensities from the green point (spec), 
											fitted spec, fitted spec component 1, fitted spec component 2,
											fitted bleaching corrected (BC) spec, fitted BC spec component 1, fitted BC spec component 2
											background spec}, 
									all ordered by lambda
		fitSwitches				.. list of True/False for all parameters; True for variable ones (to be optimized), False for fixed	
		backgroundSubtraction	.. True/False whether background was subtracted before fitting or not

output: .dat file with all the input spectra (components omitted if only one was fitted) and info about the fit (parameters, averaging, background region, etc.)
*)
ClearAll[ExportSpectra]
ExportSpectra[folder_, experimentName_, cropBackground_, avg_, parametersPoint_, spectraPoint_, fitSwitches_, backgroundSubtraction_]:=
	Block[{parameterSymbols, fit2components, table, tableSpectra, exportName, exportPath, exportFolder},

		(*add xy of the point to the name to prevent overwritting*)        
        exportName = experimentName <> "_pointXY-" <> ToString@parametersPoint[[1]] <> "-" <> ToString@parametersPoint[[2]];
                 
		(*ask for folder to export, name can be modified*)
        exportPath = SystemDialogInput["FileSave",FileNameJoin[{folder,exportName}]];
        If[ToString@exportPath=="$Canceled",Return[folder]];
        {exportFolder,exportName}=FileNameTake[exportPath,#]&/@{{1,-2},-1};

		
		(*prepare output table*)
		fit2components = fitSwitches[[5]] || parametersPoint[[5]]!=1.;
		
		parameterSymbols = {"x[px]","y[px]","i0[au]","i1[au]","p1[]","lm1[nm]","w1[nm]","a1[]","b1[/s]","i2[au]","p2[]","lm2[nm]","w2[nm]","a2[]","b2[/s]","HemOxy","chi2[N sigma^2]"};
        
		table = {
			{"Data:"},
			{experimentName},
			{If[StringLength[exportName] > StringLength[experimentName]+1,
				StringTake[exportName,{StringLength[experimentName]+1,-1}]]},
			{"\nFit:"},
			{"averaging [pix x pix]: ", avg},
			Join[{"parameters: "}, parameterSymbols], 
			Join[{""}, If[fit2components, parametersPoint, ReplacePart[parametersPoint, Partition[Range[10,15],1] -> "-"]]],
			{"\nBackground: "},
			Flatten[{"subtracted: ", backgroundSubtraction}], 
			{"area: ", "xMin","yMin","xMax","yMax"}, 
			Join[{""}, Flatten@cropBackground]
			};
			
		tableSpectra = Transpose[Round[spectraPoint]];
		tableSpectra = Join[{{"lambda[nm]", "experiment", "fit", "fitComp1", "fitComp2", "fitBC", "fitBCcomp1", "fitBCcomp2", "background"}}, tableSpectra];
		If[!fit2components, tableSpectra = tableSpectra[[All,{1,2,3,6,9}]]];

		table = Join[table, {{"\nSpectra:"}}, tableSpectra];
		
		Export[exportPath<>".xls", table, PageWidth->Infinity];
		
		Return[exportFolder];
	];




(*
ExportFileName: create a filename with an additional string as prefix/sufix and file format extension in the end

input:	name		.. string for name basis (experiment name)
		addToName	.. string to be added for diferent exported items (parameters, color-coded images, spectra, histograms, ...)
		prefixSufix	.. "Prefix" or "Sufix" - where to place the additional string
		fileType	.. string with format extension

output: all input strings appropriately joined
*)
ClearAll[ExportFileName]
ExportFileName[name_,addToName_,prefixSufix_,fileType_]:=
	Which[prefixSufix == "Prefix",addToName<>"_"<>name<>fileType,
			prefixSufix == "Sufix",name<>"_"<>addToName<>fileType
	];

(*
ExportParameters: export table of fitted parameters, preceded by rows of parameter names, fitSwitches and initialValues

input:	exportFolder			.. path of the directory to export to
		exportName				.. experiment name to include in the filename
		parameters				.. list of all parameters+chi2 from all points; bleaching rates in [s^-1]
		spectraPoints             ..spectral points for all image pixels
		specFit                   ..spectral points of 1. coomponent fit
		initialValues			.. list of initial paramter values - starting values for variable parameters or fixed for the rest; bleaching rates in [frames^-1]! 
		fitSwitches				.. list of True/False for all parameters; True for variable ones (to be optimized), False for fixed	
		timeExposure			.. exposure time of each image in seconds
		"SaveNamePrefixSufix"	.. option: "Prefix" or "Sufix" - where to place the additional string
		"SaveNameParameters"	.. option: string to be added to filename to export parameters

output: .txt file
*)
ClearAll[ExportParameters]
Options[ExportParameters] = {"SaveNamePrefixSufix" -> "Prefix", "SaveNameParameters" -> "fitParameters"};
ExportParameters[exportFolder_,exportName_,parameters_,spectraPoints_,specFit_,initialValues_,fitSwitches_,timeExposure_,opts:OptionsPattern[]] :=
    Block[ {initialValuesUnits, header, out, outName,parameters2,parameters2Pop,
        (*at the end added spectral points*)
    	parameterSymbols = Flatten[{"x[px]","y[px]","i[au]","i1[au]","p1[]","lm1[nm]","w1[nm]","a1[]","b1[/s]","i2[au]","p2[]","lm2[nm]","w2[nm]","a2[]","b2[/s]","OxyHem","chi2[N sigma^2]",Table["Ch"<>ToString[i],{i,1,Length[spectraPoints[[1]]]}],Table["1. comp fit"<>ToString[i],{i,1,Length[specFit[[1]]]}]}]},
        (*convert initial bleaching rates into s^-1*)
        initialValuesUnits = initialValues;
        initialValuesUnits[[{9,15}]]/=(timeExposure+0.2);
        (*first rows with general info*)
        header = {{"parameters:"},parameterSymbols,{"\n fitted:"},fitSwitches,{"\n starting values:"},initialValuesUnits,{"\n results:"}};
        
        (*add spectral point at the end of the table*)

        parameters2 = Flatten[Append[parameters[[#]],spectraPoints[[#]]]]&/@Range[Length[parameters]];
        parameters2Pop = Flatten[Append[parameters2[[#]],specFit[[#]]]]&/@Range[Length[parameters]];
        
Print[{"parameters2Pop[[1;;5]]",parameters2Pop[[1;;5]]}];

        (*everything together*)
        out = Join[header,parameters2Pop];
		
        (*export*)
        SetDirectory[exportFolder];
        outName = ExportFileName[exportName,OptionValue@"SaveNameParameters",OptionValue@"SaveNamePrefixSufix",".txt"];
        $stream = OpenWrite[outName,PageWidth->Infinity];
        Write[$stream,
            OutputForm[Column[
                Row[Riffle[#,"\t"]]&/@out
                ]]
            ];
        Close[$stream];
        Clear[$stream];
    ];

(*
IntensityImageCropRectangle: if only a cropped region is analized, add a red rectangle around cropped area in intensity image

input:	imageIntensity	.. ArrayPlot of all spectral images superimposed; intensity scaled from min to max
		crop			.. {{xmin,ymin},{xmax,ymax}} of the region analized [pix] 

output:	imageIntensity with red rectangle
*)
ClearAll[IntensityImageCropRectangle]
IntensityImageCropRectangle[imageIntensity_,crop_]:=
	If[crop == {{1,1},{512,512}},
		imageIntensity,
		Show[imageIntensity,Epilog->{EdgeForm[Red],Transparent,Rectangle@@crop}]
	];

(*
ImagesCCTogether: joins all color-coded images; if selected, includes intensity image and alignment images

input: 	imageIntensity				.. ArrayPlot of all spectral images superimposed; intensity scaled from min to max
		imagesCC					.. 2D matrix of Row[{color-coded image, corresponding legend}] for all fitted parameters (2 rows for 2 spectral components)
		imageIntensityCompareAlign	.. aligned intensity image of cropped area					
		"SaveImagesCCIntensity "		.. option: True/False whether to include intensity image at the beginning of the panel with color-coded images
		"SaveImagesCCAlign "			.. option: True/False whether to include alignment comparison image at the end of the panel with color-coded images	

output:	grid of all images
*)
ClearAll[ImagesCCTogether]
Options[ImagesCCTogether] = {"SaveImagesCCIntensity " -> True, "SaveImagesCCAlign " -> False};
ImagesCCTogether[imageIntensity_,imagesCC_,imageIntensityCompareAlign_,opts:OptionsPattern[]] := 
	Block[{dimensions,imagesOut},
		dimensions = Dimensions[imagesCC];
		imagesOut = imagesCC;

		If[dimensions[[1]] == 1,
			(*1 spectral component used*)
			If[OptionValue@"SaveImagesCCAlign ",
				imagesOut[[1]] = Prepend[imagesOut[[1]],imageIntensityCompareAlign];
			];
			If[OptionValue@"SaveImagesCCIntensity ",
				imagesOut[[1]] = Prepend[imagesOut[[1]],imageIntensity];
			];
			,
			(*2 spectral components used*)
			If[OptionValue@"SaveImagesCCIntensity " || OptionValue@"SaveImagesCCAlign ",
				imagesOut = Prepend[#,""]&/@imagesOut;
				If[OptionValue@"SaveImagesCCIntensity ",
					imagesOut[[1,1]] = imageIntensity;
				];
				If[OptionValue@"SaveImagesCCAlign ",
					imagesOut[[2,1]] = imageIntensityCompareAlign;
				];
			];
		];
		
		Return[imagesOut];
	];

(*
ColumnOfImagesWithFilename: put experiment name at the top of the image if selected

input: 	listOfImages	.. list of all graphics to appear in the column
		name			.. name to put at the top
		showName		.. True/False whether to include the name

output:	column of all graphics
*)
ClearAll[ColumnOfImagesWithFilename]
ColumnOfImagesWithFilename[listOfImages_,name_,showName_]:=
	Block[{listOut, out},
		listOut = listOfImages;
		If[showName,
			listOut = Prepend[listOut,Style[name,10]];
		];
		out = Column[listOut,Center,3,ItemSize->Full];
		Return[out]; 
	];


(*
ExportImages: create and export graphical output as user defined

input:	exportFolder			.. path of the directory to export to
		exportName				.. experiment name to include in the filename
		crop					.. {{xmin,ymin},{xmax,ymax}} of the region analized [pix] 
		imageIntensity         	.. ArrayPlot of all spectral images superimposed; intensity scaled from min to max
		imagesCC				.. list of Row[{color-coded image, corresponding legend}] for all fitted parameters
		histograms				.. list of histograms for chi2 and all variable parameters
		plots					.. plots of experimental and fitted spectra from 6 points with different intensities
		imageIntensityCompareAlign	.. aligned intensity image (*column of two small intensity images to compare before/after alignment*)
		fitSwitches				.. list of True/False for all parameters; True for variable ones (to be optimized), False for fixed	
		saveImagesCC			.. True/False whether to save color-coded images
		saveHistograms			.. True/False whether to save histograms
		savePlots				.. True/False whether to save plots with spectra
		"SaveNamePrefixSufix "	.. option: "Prefix " or "Sufix " - where to place the additional string
		"SaveNameImagesCC " 		.. option: string to be added to filename when exporting color-coded images
		"SaveNameHistograms "	.. option: string to be added to filename when exporting histograms	
		"SaveNameSpectra "		.. option: string to be added to filename when exporting plots with spectra
		"SaveNameAll "			.. option: string to be added to filename when exporting all graphics as one panel 
		"SaveImagesCCTogether "	.. option: True/False whether to save all color-coded images as one panel or separately
		"SaveImagesCCIntensity " .. option: True/False whether to include intensity image at the beginning of the panel with color-coded images
		"SaveImagesCCAlign "		.. option: True/False whether to include alignment comparison image at the end of the panel with color-coded images	
		"SaveAllImagesTogether "	.. option: True/False whether to save all graphics as one panel
		"SaveShowFileName "		.. option: True/False whether to include exportName at the top of the panel
		"SaveImageFiletype "       .. option: string for exported graphics filetype: ".png", ".jpg", ".gif", ".tif"
		"SaveImageResolution "     .. option: number for exported graphics resolution in dpi

output: none
*)
ClearAll[ExportImages]
Options[ExportImages] = {
	"SaveNamePrefixSufix " -> "Prefix ", 
	"SaveNameImagesCC " -> "ccImages ", 
	"SaveNameHistograms " -> "histograms ", 
	"SaveNameSpectra " -> "spectra ", 
	"SaveNameAll " -> "results ",
	"SaveImagesCCTogether " -> True, 
	"SaveImagesCCIntensity " -> True, 
	"SaveImagesRaw " -> False,
	"SaveImagesCCAlign " -> False,
	"SaveAllImagesTogether " -> False, 
	"SaveShowFileName " -> True,
	"SaveImageFiletype "->".png",
	"SaveImageResolution "->150};
ExportImages[exportFolder_,exportName_,crop_,imageIntensity_,imagesCC_,histograms_,plots_,imageIntensityCompareAlign_,fitSwitches_,
	{saveImagesCC_,saveHistograms_,savePlots_},opts:OptionsPattern[]] := 
	Block[{imagesCCout,histogramsOut,out,outRaw,outRawCC,outName,outNameRaw,outNameRawCC,parameterNames,imagesCCPop},
		
		SetDirectory[exportFolder];
		
        If[ OptionValue@"SaveAllImagesTogether " ,
        	
        	(*save all selected plots and images as one panel*)
        	
        	imagesCCout = ImagesCCTogether[imageIntensity,imagesCC,imageIntensityCompareAlign,
        		"SaveImagesCCIntensity " -> OptionValue@"SaveImagesCCIntensity ", 
        		"SaveImagesCCAlign " -> OptionValue@"SaveImagesCCAlign "];
        	histogramsOut = ImagesCCTogether[Null,histograms,Null,
        		"SaveImagesCCIntensity " -> OptionValue@"SaveImagesCCIntensity ", 
        		"SaveImagesCCAlign " -> OptionValue@"SaveImagesCCAlign "];
        	out = Pick[{imagesCCout,histogramsOut,plots},{saveImagesCC,saveHistograms,savePlots}];
        	(*if saving CC-images and histograms, align them into a grid*)
        	If[ saveImagesCC&&saveHistograms,
        		out[[1]] = Flatten[out[[;;2]],1];
        		out[[1]] = Grid[out[[1]],Spacings->{3,3},Alignment->Left,ItemSize->Full];
        		out = Drop[out,{2}];
        	];
        	out = ColumnOfImagesWithFilename[out,exportName,OptionValue@"SaveShowFileName "]; 
        	outName = ExportFileName[exportName,OptionValue@"SaveNameAll ",OptionValue@"SaveNamePrefixSufix ",OptionValue@"SaveImageFiletype "];
        	Export[outName,out,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];
        		
        	,(*else save separately*)
        	
	        If[ saveHistograms,
	        	out = ColumnOfImagesWithFilename[{Grid[histograms,Spacings->{3,3},ItemSize->Full]},exportName,OptionValue@"SaveShowFileName "];
	        	outName = ExportFileName[exportName,OptionValue@"SaveNameHistograms ",OptionValue@"SaveNamePrefixSufix ",OptionValue@"SaveImageFiletype "];
	        	Export[outName,out,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];
	        ];
	        
	        If[ savePlots,
	            out = ColumnOfImagesWithFilename[{plots},exportName,OptionValue@"SaveShowFileName "]; 
        		outName = ExportFileName[exportName,OptionValue@"SaveNameSpectra ",OptionValue@"SaveNamePrefixSufix ",OptionValue@"SaveImageFiletype "];
	        	Export[outName,out,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];
	        ];
	        
	        If[ saveImagesCC,
	        	If[OptionValue@"SaveImagesCCTogether ",
	        		(*all color-coded images together*)
	        		imagesCCout = ImagesCCTogether[imageIntensity,imagesCC,imageIntensityCompareAlign,
        						"SaveImagesCCIntensity " -> OptionValue@"SaveImagesCCIntensity ", 
        						"SaveImagesCCAlign " -> OptionValue@"SaveImagesCCAlign "];
	        		out = ColumnOfImagesWithFilename[{Grid[imagesCCout,Spacings->{3,3},ItemSize->Full]},exportName,OptionValue@"SaveShowFileName "];  
	        		outName = ExportFileName[exportName,OptionValue@"SaveNameImagesCC ",OptionValue@"SaveNamePrefixSufix ",OptionValue@"SaveImageFiletype "];
	        		Export[outName,out,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];
	        		,

	        		(*else save each color-coded image separately*)

					parameterNames = Pick[{"x ","y ","i ","i1 ","p1 ","lm1 ","w1 ","a1 ","b1 "(*,"i2 "*),"HemOxy","p2 ","lm2 ","w2 ","a2 ","b2 ","HemOxy","chi2"},
	        				If[fitSwitches[[5]],ReplacePart[fitSwitches,Partition[{4,10,11},1]->True],ReplacePart[fitSwitches,4->True]](*save also image for i1,2 and p2*)
	        				];
					Print[{"parameterNames",parameterNames}];
					imagesCCPop=Flatten[imagesCC];

	        		Do[
	        			imagesCCout = ImagesCCTogether[imageIntensity,{Cases[Flatten[imagesCC],Except[Null]][[{i}]]},imageIntensityCompareAlign,
        						"SaveImagesCCIntensity " -> OptionValue@"SaveImagesCCIntensity ", 
        						"SaveImagesCCAlign " -> OptionValue@"SaveImagesCCAlign "];
	        			out = ColumnOfImagesWithFilename[{Grid[imagesCCout,Spacings->{3,3},ItemSize->Full]},exportName,OptionValue@"SaveShowFileName "];  
	        			outName = ExportFileName[exportName,(OptionValue@"SaveNameImagesCC ")<>"-"<>parameterNames[[i]],OptionValue@"SaveNamePrefixSufix ",OptionValue@"SaveImageFiletype "];
						Export[outName,out,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];	        		
						
						outRaw=ImageCrop[imageIntensity,Length[ImageData[imageIntensity]]-20];
						outNameRaw=exportName<>"Intensity"<>OptionValue@"SaveImageFiletype ";	
						Export[outNameRaw,outRaw,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];

						outRawCC=ImageCrop[ImageResize[imagesCCPop[[i,1,1]],512],Length[ImageData[imageIntensity]]-20];
						outNameRawCC=ExportFileName[exportName<>"Raw",(OptionValue@"SaveNameImagesCC ")<>"-"<>parameterNames[[i]],OptionValue@"SaveNamePrefixSufix ",OptionValue@"SaveImageFiletype "];
						Export[outNameRawCC,outRawCC,ImageResolution->OptionValue@"SaveImageResolution "(*,ImageSize->All*)];	

	        		,{i,{1,3,5,4}} (*\[CHacek]e fitamo tudi \[SHacek]irino spektra moramo tu spremeniti tabelo!! \[CapitalCHacek]e fitamo w imamo 5 stevil {1,3,5,4,6}, \[CHacek]e ne fitamo w jih imamo 4!*)
	        		];
	        	];
	        ];
        ];
        
 	];

(*
ExportResults: ask for filename and export all the results as user defined

input:	folder						.. full path where to set the file-saving dialog
		name						.. experiment name to include in the filename
		crop						.. {{xmin,ymin},{xmax,ymax}} of the region analized [pix] 
		timeExposure				.. exposure time of each image in seconds
		fitSwitches							.. list of True/False for all parameters; True for variable ones (to be optimized), False for fixed	
		initialValues				.. list of initial paramter values - starting values for variable parameters or fixed for the rest; bleaching rates in [frames^-1]!
		parameters					.. list of all parameters+chi2 from all points; bleaching rates in [s^-1]
		spectraPoints                 ..spectral points for all image pixels
		histograms					.. list of histograms for chi2 and all variable parameters
		plots						.. plots of experimental and fitted spectra from 6 points with different intensities
		intensityData				.. data of all spectral images superimposed, normalized to the maximal value
		imageIntensityCompareAlign	.. aligned intensity image (*column of two small intensity images to compare before/after alignment*)
		imagesCC					.. 2D matrix of Row[{color-coded image, corresponding legend}] for all fitted parameters (2 rows for 2 spectral components)
		saveParameters				.. True/False whether to save the table of fitted parameters as a .txt file
		saveHistograms				.. True/False whether to save histograms
		savePlots					.. True/False whether to save plots with spectra
		saveImagesCC				.. True/False whether to save color-coded images
		adjustContrastMinMax     	.. list of two values between 0 and 1 to scale contrast for intensityImage: {minIntensity,maxInensity}
		"SaveNamePrefixSufix "		.. option: "Prefix " or "Sufix " - where to place the additional string
		"SaveNameParameters "		.. option: string to be added to filename to export parameters
		"SaveNameImagesCC " 			.. option: string to be added to filename when exporting color-coded images
		"SaveNameHistograms "		.. option: string to be added to filename when exporting histograms	
		"SaveNameSpectra "			.. option: string to be added to filename when exporting plots with spectra
		"SaveNameAll "				.. option: string to be added to filename when exporting all graphics as one panel 
		"SaveImagesCCTogether "		.. option: True/False whether to save all color-coded images as one panel or separately
		"SaveImagesCCIntensity " 	.. option: True/False whether to include intensity image at the beginning of the panel with color-coded images
		"SaveImagesCCIntensityCrop " .. option: True/False whether to crop intensity image to the same area as color-coded images
        "SaveImagesRaw "              .. option: True/False whether to export intensity image and analysed CC images separately
		"SaveImagesCCAlign "			.. option: True/False whether to include alignment comparison image at the end of the panel with color-coded images	
		"SaveAllImagesTogether "		.. option: True/False whether to save all graphics as one panel
		"SaveShowFileName "			.. option: True/False whether to include exportName at the top of the panel
		"SaveImageFiletype "       	.. option: string for exported graphics filetype: ".png", ".jpg", ".gif", ".tif"
		"SaveImageResolution "     	.. option: number for exported graphics resolution in dpi
		"BatchModeQ "                  .. option: True/False to control behavior in batch analysis mode

output: exportFolder				.. full path of the selected export folder ; returned to be used as input at next ExportResults
*)
ClearAll[ExportResults]
Options[ExportResults] = {"SaveNamePrefixSufix" -> "Prefix", "SaveNameParameters" -> "fitParameters", 
	"SaveNameImagesCC" -> "ccImages", "SaveNameHistograms" -> "histograms", "SaveNameSpectra" -> "spectra", "SaveNameAll" -> "results",
	"SaveImagesCCTogether" -> True, "SaveImagesCCIntensity" -> True, "SaveImagesCCIntensityCrop" -> False, "SaveImagesRaw" -> False, "SaveImagesCCAlign" -> False,
	"SaveAllImagesTogether" -> False, "SaveShowFileName" -> True,
	"SaveImageFiletype"->".png","SaveImageResolution"->150,
	"BatchModeQ"->False};
ExportResults[folder_,name_,crop_,timeExposure_,{fitSwitches_,initialValues_},
		{parameters_,spectraPoints_,specFit_,histograms_,plots_,intensityData_,imageAlign_,imagesCC_},
		{saveParameters_,saveHistograms_,savePlots_,saveImagesCC_},
		adjustContrastMinMax_,
		opts:OptionsPattern[]] :=
    Block[ {exportPath,exportFolder,exportName,imageIntensity},
        
        (*add crop info to the name to prevent overwritting*)        
        exportName = If[ crop=={{1,1},{512,512}},
                         name,
                         name<>"_crop-"<>StringJoin[Riffle[ToString/@Flatten[crop],"-"]]
                     ];
		
		(*if not in batch mode, ask for folder to export, name can be modified*)
        exportFolder = folder;
		exportPath = FileNameJoin[{exportFolder,exportName}];       
        If[!OptionValue@"BatchModeQ",
            exportPath = SystemDialogInput["FileSave",exportPath];
            If[ToString@exportPath=="$Canceled",Return[folder]];
            {exportFolder,exportName}=FileNameTake[exportPath,#]&/@{{1,-2},-1};
		];
		        
        (*export parameters*)
        If[saveParameters,
            ExportParameters[exportFolder,exportName,parameters,spectraPoints,specFit,initialValues,fitSwitches,timeExposure,
            	"SaveNamePrefixSufix" -> OptionValue@"SaveNamePrefixSufix", 
            	"SaveNameParameters" -> OptionValue@"SaveNameParameters"];
        ];
        
        (*export images*)
		If[saveImagesCC && OptionValue@"SaveImagesCCIntensity",
			If[OptionValue@"SaveImagesCCIntensityCrop",
				imageIntensity = IntensityImage[CropImages[{intensityData},crop][[1]],adjustContrastMinMax];
				,
				imageIntensity = IntensityImage[intensityData,adjustContrastMinMax];
				imageIntensity = IntensityImageCropRectangle[imageIntensity,crop];
			]
		];

        If[saveHistograms||savePlots||saveImagesCC,
        	ExportImages[exportFolder,exportName,crop,imageIntensity,imagesCC,histograms,plots,imageAlign,fitSwitches,{saveImagesCC,saveHistograms,savePlots},
        		"SaveNamePrefixSufix " -> OptionValue@"SaveNamePrefixSufix",
        		"SaveNameImagesCC " -> OptionValue@"SaveNameImagesCC",
        		"SaveNameHistograms " -> OptionValue@"SaveNameHistograms",
        		"SaveNameSpectra " -> OptionValue@"SaveNameSpectra",
        		"SaveNameAll " -> OptionValue@"SaveNameAll",
        		"SaveImagesCCTogether " -> OptionValue@"SaveImagesCCTogether",
        		"SaveImagesCCIntensity " -> OptionValue@"SaveImagesCCIntensity",
				"SaveImagesRaw " -> OptionValue@"SaveImagesRaw",
        		"SaveImagesCCAlign " -> OptionValue@"SaveImagesCCAlign",
        		"SaveAllImagesTogether " -> OptionValue@"SaveAllImagesTogether",
        		"SaveShowFileName " -> OptionValue@"SaveShowFileName",
				"SaveImageFiletype " -> OptionValue@"SaveImageFiletype",
				"SaveImageResolution " -> OptionValue@"SaveImageResolution"
        	];
        ];
        
        Return[exportFolder];
    ]










(* ::Section:: *)
(*(* FOR MANIPULATE *)*)


(*
AnalizeOnePoint: get SNR from green point, fit that spectrum, and plot the graph with green (point) and blue (background) spectra for manipulate    

input:	images					.. 3D matrix with 16-bit data from all images in the spectral series; sorted by wavelength, as saved
		cropPoint				.. {{xmin,ymin},{xmax,ymax}} of the region to calculate signal from (green point) [pix] 
		cropBackground			.. {{xmin,ymin},{xmax,ymax}} of the region to calculate background from (blue rectangle)[pix] 
		randomSampling			.. True for experiment with stochastic wavelength sampling, False for linear
		lambdaSample			.. list of wavelengths as sampled in the experiment
		timeExposure			.. exposure time of each image in seconds
		signalThreshold			.. intensity threshold to filter points with signal < signalThreshold * noise
		parameterSwitches		.. list of True/False for all parameters; True for variable ones (those to be optimized), False for fixed
		parameterFixedValues	.. list of values for all parameters; used as initial points for optimization of variable parameters, kept fixed for the rest; bleachRates in [frames^-1]
		parameterLimits			.. list of limits to keep parameters within during optimization: {{par1min,par1max},..} for all 16 parameters; bleachRates in [frames^-1]
		bleachRefMode			.. True for experiment with a reference wavelength for bleaching correction (used with linear sampling)
		bleachRefImages			.. 3D matrix with 16-bit data from images in the bleaching reference series; sorted by time taken
		bleachRefLambda			.. wavelength at which bleaching reference images were recorded
		bleachRefFrameNumbers	.. consecutive numbers of bleaching reference images in the whole series, e.g. {1,4,7,11,16,23} (not always equally spaced)  

output: noise					.. standard deviation of background signal
		parametersPoint			.. parameters+chi2 from fitted green point; bleaching rates in [s^-1]
		plotSpectrum			.. plot with experimental and fitted spectra of green point and background for manipulate 
		plotBleachRef			.. plot with bleaching reference time series from the green point; "" if reference data not measured
		parametersPointTable	.. table of fitted parameter values from the green poin, needed for display in manipulate 
		spectraPoint			.. table of numerical values for all plotted spectra, needed for saving
*)
ClearAll[AnalizeOnePoint]
AnalizeOnePoint[images_,cropPoint_, cropBackground_, randomSampling_, lambdaSample_, timeExposure_, signalThreshold_,
	{parameterSwitches_,parameterFixedValues_,parameterLimits_},{bleachRefMode_, bleachRefImages_, bleachRefLambda_, bleachRefFrameNumbers_}] :=
    
    Block[{snr, noise, spectrumBackground, lambdaOrder, lambdaShuffle, nAvg, 
    	spectrumBlock, spectrum, parametersPoint, parametersPointTable,
    	spectrumFit, spectrumFitComp1, spectrumFitComp2, spectrumFitNoBleach, spectrumFitNoBleachComp1, spectrumFitNoBleachComp2,
    	parameterNamesUnits, plotLabel, plotRange, plotSpectrum, 
    	bleachingSeries, bleachRate, BleachFit, plotBleachRef, spectraPoint},
    	 
    	(*signal from the selected (green) point, noise as StDev in background (blue rectangle)*)
    	{snr, noise} = SNRfromImages[images, cropPoint, cropBackground];
         
        (*spectrum in the background*)
        spectrumBackground = SpectrumAvg[images, cropBackground];
         

        (*preparation for fitting*)
		lambdaOrder = Ordering[lambdaSample];
		lambdaShuffle = Flatten[Position[Sort[lambdaSample],#]&/@lambdaSample];
		nAvg = Round[(Abs[cropPoint[[2,1]]-cropPoint[[1,1]]])/2];
		
			
		(* *** spectrum *** *)
		
		(*spectrumBlock contains all info needed for fitting: 
		{{list of parameters, symbols in places of variable ones},{spectrum},snr,fitSwitch,{bleachingSeries},{initialValues for variable parameters}}*)
		spectrumBlock = SpectraFromImages[CropImages[images,cropPoint],cropPoint,lambdaShuffle,nAvg,noise,0(*signalThreshold=0 -> fit always*),
		 				{bleachRefMode,If[bleachRefMode,CropImages[bleachRefImages,cropPoint],{}]},Count[parameterSwitches,True]];
		spectrum = spectrumBlock[[1,2,lambdaOrder]]; (*first (=only) element in the list, second element in the block (=spectrum)*)
		
(*Print[{"spectrum",spectrum}];*)


		(* *** fit *** *)
		
		parametersPoint = AnalizeSpectra[spectrumBlock,lambdaSample,lambdaOrder,noise,timeExposure,
						parameterSwitches,parameterFixedValues,parameterLimits,
						{bleachRefMode,bleachRefLambda,bleachRefFrameNumbers},
	 					"UseParallel" -> False, "GeneratePlots" -> False][[1,1]];(*first (=only) element in the list, parameters at position 1*)
(*Print["err"];*)
(*Print[{"parametersPoint",parametersPoint}];	*) 
		
        (*convert from s^-1 into frames^-1 for SpectraSumC*)
		parametersPoint[[{9,15}]]*=(timeExposure+.2);
		
		(*get fitted spectra - ordered by lambda*)         
        spectrumFit = SpectraSumC[parametersPoint,lambdaSample][[lambdaOrder]];
        spectrumFitComp1 = SpectraSumC[ReplacePart[parametersPoint,{3->parametersPoint[[4]],5->1,11->0}],lambdaSample][[lambdaOrder]];
        spectrumFitComp2 = SpectraSumC[ReplacePart[parametersPoint,{3->parametersPoint[[10]],5->0,11->1}],lambdaSample][[lambdaOrder]];
        spectrumFitNoBleach = SpectraSumC[ReplacePart[parametersPoint,{9->0,15->0}],lambdaSample][[lambdaOrder]];
        spectrumFitNoBleachComp1 = SpectraSumC[ReplacePart[parametersPoint,{3->parametersPoint[[4]],5->1,11->0,9->0,15->0}],lambdaSample][[lambdaOrder]];
        spectrumFitNoBleachComp2 = SpectraSumC[ReplacePart[parametersPoint,{3->parametersPoint[[10]],5->0,11->1,9->0,15->0}],lambdaSample][[lambdaOrder]];
         
        
        (* *** plots *** *)
         
        plotRange = {-1, 1.1 Mean[Sort[Flatten[{spectrum, spectrumBackground, spectrumFitNoBleach}]][[-5;;]]]};
        
        (*plot of reference bleaching series*)
        If[bleachRefMode,
		
			bleachingSeries = Transpose[{bleachRefFrameNumbers,spectrumBlock[[1,5]]}];
			bleachRate = parametersPoint[[9]];
			BleachFit[n_] := bleachingSeries[[1,2]] Exp[-bleachRate (n-1)];
						
	   		plotBleachRef = Show[
		    	ListPlot[bleachingSeries, PlotStyle -> Directive[Red,PointSize[.02]], PlotRange -> {{0,All},plotRange},
		    		Frame -> True, FrameLabel -> {"frame N","I [a.u.]"}, LabelStyle -> {FontSize -> 14}, 
		    		Epilog -> Inset[Row[{Style[Subscript["\[Lambda]","ref"],14], Style[" = " <> ToString[bleachRefLambda] <> " nm",14]}], ImageScaled[{.38,.35}]],
		    		ImageSize -> 300(*, AspectRatio -> 1*)],
		    	Plot[BleachFit[n],{n, 1, bleachRefFrameNumbers[[-1]]},PlotRange -> All, PlotStyle -> Red]
	    	,ImagePadding->{{60,5},{60,5}}];
	    	,
			plotBleachRef = "";
		];
         
        (*convert from frames^-1 back into s^-1*) 
        parametersPoint[[{9,15}]]/=(timeExposure+.2);
        
        (*list of parameters {par, unit} - used for plotLabel*) 
        parameterNamesUnits = StringCases[#, x__~~" ["~~y___~~"]" :> {x,(*" "<>*)y}][[1]]&/@$parameterSymbolsForm;
       
        
        (*table of parameters' values - Red if fitted, else Gray; if only 1 component used, replace 2nd by "-"*)        	        
		parametersPointTable = Grid[
			{
			(*x,y,I0,chi2*)			
			parameterNamesUnits[[{1,2,3,-1},1]], 
			parametersPoint[[{1,2,3,-1}]],
			{""},		
			(*1st spectral component*)	
			parameterNamesUnits[[4;;9,1]],
			MapThread[
				Style[#1,If[#2,Red,Gray]]&,
				{parametersPoint[[4;;9]],parameterSwitches[[4;;9]]}
			],
			{""},
			(*2nd spectral component*)	
			parameterNamesUnits[[10;;-2,1]],
			If[parameterSwitches[[5]]||(parameterFixedValues[[5]]<1),
				MapThread[
					Style[#1,If[#2,Red,Gray]]&,
					{parametersPoint[[10;;-2]],parameterSwitches[[10;;-2]]}
				],
				{"-","-","-","-","-","-","-"}
			]
			}
			];
		
        (*plot of spectra*) 
		 plotSpectrum = {ListLinePlot[{ spectrumBackground(*, spectrumFitNoBleach*),spectrumFitNoBleachComp1, spectrumFit}, 
           						PlotStyle -> { Blue, Directive[Green, Dashing[{.04,.02}]], Green}, 
           						PlotRange -> plotRange, DataRange -> {Min[lambdaSample], Max[lambdaSample]}
           						],
               		ListPlot[spectrum, PlotStyle -> Directive[Green,PointSize[.02]], PlotRange -> plotRange, DataRange -> {Min[lambdaSample], Max[lambdaSample]}]
         			};       
		
		If[parameterSwitches[[5]]||(parameterFixedValues[[5]]<1),
					PrependTo[plotSpectrum,
								ListLinePlot[{ spectrumFitNoBleachComp1, spectrumFitNoBleachComp2, spectrumFitComp1,spectrumFitComp2}, 
           						PlotStyle -> { Directive[Green, Dotted], Directive[Green, Dotted], Directive[Green,Dashed], Directive[Green,Dashed]}, 
           						PlotRange -> plotRange, DataRange -> {Min[lambdaSample], Max[lambdaSample]}
           						]];
		];
 
		plotSpectrum = Tooltip[
						Show[plotSpectrum,
							Frame -> True, FrameLabel -> {"\[Lambda] [nm]", "I [a.u.]"}, LabelStyle -> {FontSize -> 14},
              				(*PlotLabel -> plotLabel,*)
               				ImageSize -> 300,(* AspectRatio -> 1,*)ImagePadding->{{60,5},{60,5}},
               				Epilog -> {Gray,Line[{{Min[lambdaSample], noise signalThreshold}, {Max[lambdaSample],noise signalThreshold}}]}],
         				parametersPointTable];
		
		spectraPoint={Sort@lambdaSample, spectrum, spectrumFit, spectrumFitComp1, spectrumFitComp2, 
				spectrumFitNoBleach, spectrumFitNoBleachComp1, spectrumFitNoBleachComp2, spectrumBackground};
				
		Return[{noise, parametersPoint, plotSpectrum, plotBleachRef, parametersPointTable, spectraPoint}];	
    ];


(* NOT USED
(*
ManualSetValues: replace manualSet parameter values with fitted values from the green point, unless locked by manualSwitch

input:  manualSetValues .. list of numbers for all fitable parameters {manualSetPortion1,manualSetLambdaMax1,manualSetWidth1,manualSetAsymmetry1,manualSetBleach1,manualSetLambdaMax2,...}
        manualSwitches  .. list of True/False for all fitable parameters - where True, the corresponding manualSet value will not be replaced
        parametersPoint .. list of all parameter values from the fitted green point {x,y,int,p1,lm1,w1,a1,b1,p2,lm2,w2,a2,b2}

output: manualSetValuesNew .. same as manualSetValues, with some values taken from parametersPoint 
*)
ClearAll[ManualSetValues]
ManualSetValues[manualSetValues_,manualSwitches_,parametersPoint_]:=
	Block[{manualSetValuesNew,indexesToReset},
		manualSetValuesNew = manualSetValues;
		indexesToReset = Flatten[Position[manualSwitches,False]];
		manualSetValuesNew[[indexesToReset]]=parametersPoint[[{4,5,6,7,8,10,11,12,13}[[indexesToReset]]]];

		Return[manualSetValuesNew]
	];
*)




(*
ControlsFit: gathers manipulate controls for fitting into lists for all the parameters
			input bleaching rates in [s^-1], output in [frames^-1]

input:	fitNComponents			.. number of spectral components to be used for optimization: 1 or 2
		fitSwitchLambdaMax1		.. True/False whether to fit wavelength of max intensity for first component
		fitSwitch___			.. same for all other parameters
								.. 
		manualSetLambdaMax1 	.. value of lambdaMax that should be used as initial point (for variable parameters) or fixed value (for nonvariable)
		manualSet___			.. same for all other parameters (bleachRates in [s^-1])
								..
		limitLambdaMaxLow		.. lower limit for lambdaMax optimization
		limitLambdaMaxHigh		.. upper limit for lambdaMax optimization
		limit___ ___ ___Low		.. lower limit for optimization of all other parameters (bleachRates in [s^-1])
		limit___ ___ ___High		.. upper limit for optimization of all other parameters (bleachRates in [s^-1])
								..
		exposureTime			.. exposure time of each image in seconds
		initialPoint			.. parameters from green point spectrum optimization (bleachRates in [s^-1])

output:	fitSwitches				.. list of True/False for all parameters; True for variable ones (those to be optimized), False for fixed
		initialValues			.. list of values for all parameters; used as initial points for optimization of variable parameters, kept fixed for the rest
		parameterLimits			.. list of limits to keep parameters within during optimization: {{par1min,par1max},..} for all 15 parameters; 
									bleachRates in [frames^-1]
*)
ClearAll[ControlsFit]
ControlsFit[fitNComponents_,
	fitSwitchPortion1_, fitSwitchLambdaMax1_, fitSwitchWidth1_, fitSwitchAssymetry1_, fitSwitchBleach1_,
	fitSwitchLambdaMax2_, fitSwitchWidth2_, fitSwitchWidth2_, fitSwitchAssymetry2_, fitSwitchBleach2_, oxygenation_,
	manualSetPortion1_,
	manualSetLambdaMax1_, manualSetWidth1_,	manualSetAsymmetry1_, manualSetBleach1_,
	manualSetLambdaMax2_, manualSetWidth2_,	manualSetAsymmetry2_, manualSetBleach2_,manualSetOxygen_,
	limitPortionLow_, limitPortionHigh_, limitLambdaMaxLow_, limitLambdaMaxHigh_, limitWidthLow_,limitWidthHigh_,
	limitAsymmetryLow_,limitAsymmetryHigh_, limitBleachLow_, limitBleachHigh_, limitOxyLow_, limitOxyHigh_,
	exposureTime_] := 

	Block[{fitSwitches, initialValues, parameterLimits},
		
		fitSwitches = {(*x,y,intensity,intensity1 - never fit*)False, False, False, False,
	       	   (*portion1*)	fitSwitchPortion1 	&& fitNComponents == 2,
	       (*lm1,w1,a1,b1*) fitSwitchLambdaMax1, 
	       					fitSwitchWidth1, 
	       					fitSwitchAssymetry1, 
	       					fitSwitchBleach1,
	     (*portion2, int2*)	False, False,
	       (*lm2,w2,a2,b2*) fitSwitchLambdaMax2 && fitNComponents == 2, 
	       					fitSwitchWidth2 	&& fitNComponents == 2, 
	       					fitSwitchAssymetry2 && fitNComponents == 2,
	       					fitSwitchBleach2 	&& fitNComponents == 2,
							   oxygenation       && fitNComponents == 2,
	       		   (*chi2*) False};

	    initialValues = {(*x,y,intensity*)1, 1, 1.,
	    	(*int1,p1,lm1,w1,a1,b1*)1., manualSetPortion1,	  manualSetLambdaMax1, manualSetWidth1, manualSetAsymmetry1, manualSetBleach1,
	    	(*int2,p2,lm2,w2,a2,b2,HemOxy*)1., 1.-manualSetPortion1, manualSetLambdaMax2, manualSetWidth2, manualSetAsymmetry2, manualSetBleach2, manualSetOxygen};


		If[ fitNComponents == 1, 
			initialValues[[{5,11}]] = {1.,0.};
		];
		If[ fitNComponents == 2 && manualSetPortion1 == 1., 
			initialValues[[{5,11}]] = {.5,.5};
		];
		
		parameterLimits = {(*x,y,intensity,int1 - no meaning, just to keep lists of the same size*){0,1}, {0,1}, {0,1}, {0,1},
	       (*portion1*)	{limitPortionLow, 	limitPortionHigh},
	   (*lm1,w1,a1,b1*)	{limitLambdaMaxLow, limitLambdaMaxHigh}, 
	   					{limitWidthLow, 	limitWidthHigh}, 
	       				{limitAsymmetryLow, limitAsymmetryHigh},
	       				{limitBleachLow, 	limitBleachHigh},
	       	   (*int2*) {0,1},
	       (*portion2*) {1.-limitPortionHigh,1.-limitPortionLow},
	   (*lm2,w2,a2,b2,HemOxy*) {limitLambdaMaxLow, limitLambdaMaxHigh}, 
	   					{limitWidthLow, 	limitWidthHigh}, 
	       				{limitAsymmetryLow, limitAsymmetryHigh}, 
	       				{limitBleachLow, 	limitBleachHigh},
						   {limitOxyLow, 	limitOxyHigh}};
		

		(*convert time units from s^-1 into frames^-1*)
		initialValues[[{9,15}]] *= (exposureTime + .2);
	    parameterLimits[[{9,15}]] *= (exposureTime + .2);
		
		(*Print[{fitSwitches, initialValues, parameterLimits}];*)

		Return[{fitSwitches, initialValues, parameterLimits}];
	
	];

(*
ControlsColor: gathers manipulate controls for color coding one list for all the parameters

input: 	separate lower and upper limits for all parameters (bleaching rates in s^-1); intensity "limits" work differently -> see comments at ColorCodedImage
		parameters	.. list of all fitted parameters+chi2 from all points; bleaching rates in [s^-1]
		fitNComponents  .. number of spectral components to be fitted (1 or 2)

output: list of pairs {par1min,par1max} for all parameters (bleaching rates in s^-1)
*)
ClearAll[ControlsColor]
ControlsColor[(*{colorLimitIntensityLow_, colorLimitIntensityHigh_},*)
			{colorLimitIntensity1Low_, colorLimitIntensity1High_},
			{colorLimitIntensity2Low_, colorLimitIntensity2High_},
			{colorLimitPortionLow_, colorLimitPortionHigh_},
			{colorLimitLambdaMaxLow_, colorLimitLambdaMaxHigh_},
			{colorLimitWidthLow_, colorLimitWidthHigh_},
			{colorLimitAsymmetryLow_,colorLimitAsymmetryHigh_},
			{colorLimitBleachLow_, colorLimitBleachHigh_},
			{colorLimitOxyLow_,colorLimitOxyHigh_},
			parameters_,
			fitNComponents_
			]:=
	
	Block[{(*colorRangeIntensity,*) colorRangePortion, colorRangeLambda, colorRangeWidth, colorRangeAsymmetry, colorRangeBleach,colorRangeHemOxy, colorRangeIntensity1,colorRangeIntensity2, colorRanges},
		
		(*colorRangeIntensity = 	{colorLimitIntensityLow,	colorLimitIntensityHigh};*)
		colorRangePortion = 	{colorLimitPortionLow, 		colorLimitPortionHigh};
		colorRangeLambda = 		{colorLimitLambdaMaxLow, 	colorLimitLambdaMaxHigh};
		colorRangeWidth = 		{colorLimitWidthLow, 		colorLimitWidthHigh};
		colorRangeAsymmetry = 	{colorLimitAsymmetryLow,	colorLimitAsymmetryHigh};
		colorRangeBleach = 		{colorLimitBleachLow, 		colorLimitBleachHigh};
		colorRangeIntensity1 = 	{colorLimitIntensity1Low,	colorLimitIntensity1High};
		colorRangeIntensity2 = 	{colorLimitIntensity2Low,	colorLimitIntensity2High};
		colorRangeHemOxy    =      {colorLimitOxyLow,colorLimitOxyHigh};
		

		(*if partial intensities are set to default values, scale them according to the fitted values*)
		If[{colorLimitIntensity1Low, colorLimitIntensity1High} == {0., 1.} || ({colorLimitIntensity2Low, colorLimitIntensity2High} == {0., 1.}&&fitNComponents==2),
  			With[{parameters2show = Select[parameters, 0<#[[3]](*<2^16-1*)&]},
  				(*{colorRangeIntensity1, colorRangeIntensity2} = Round[Sort[parameters2show[[All,#]]][[Round[{.001, .999} Length[parameters2show]]+{1,0}]]]&/@{4, 10}*)
					colorRangeIntensity1 = {0,N[Max[Round[Sort[parameters2show[[All,#]]][[Round[.999 Length[parameters2show]]+{1,0}]]]&/@{4, 10}]/1.2]};
					colorRangeIntensity2 = N[colorRangeIntensity1/7];
  			]
  		];
  				
		colorRanges = {(*x,y*){0., 1.}, {0., 1.},(*intensity*){0,1}(*colorRangeIntensity*),
			(*int1*)colorRangeIntensity1,(*p1*)colorRangePortion,			(*lm1,w1,a1,b1*)colorRangeLambda, colorRangeWidth, colorRangeAsymmetry, colorRangeBleach,
   			(*int2*)colorRangeIntensity2,(*p2*)1.-Reverse@colorRangePortion,(*lm2,w2,a2,b2*)colorRangeLambda, colorRangeWidth, colorRangeAsymmetry, colorRangeBleach,colorRangeHemOxy
   			};

   		
   		(*if upper and lower limit the same, avoid division by zero*)	
   		If[Subtract@@(colorRanges[[#]])==0, colorRanges[[#]]+={0,1}]&/@Range[Length[colorRanges]];
   			
   		Return[colorRanges];
	];
  


  


(* ::Section:: *)
(*(* BOOKMARKS *)*)


(*bookmarks contain settings for all controls in Manipulate; accessible through (+) in top right corner of the panel, stored in Typeset`bookmarkList$$*)

(*
BookmarkContents: extract names and values of parameters from each bookmark in the list

input: 	bookmarkList		.. list of bookmarks in FullForm, converted to string: ToString@FullForm@Typeset`bookmarkList$$

output:	bookmarkNames		.. list of bookmark names
		bookmarkSettings	.. list of lists {"parameter1Name = parameter1Value", ...} for all bookmarks
*)
ClearAll[BookmarkContents]
BookmarkContents[bookmarkList_] :=
    Block[ { bookmarks, bookmarkNames, bookmarkSettings},
     	(*split the whole string into separate bookmarks; 
     	they come in the form RuleDelayed["bookmarkName", List[Set[parameter=value], ... ,Set[...]]];
     	x takes the name, y a string of Set rules*)
     	bookmarks = StringCases[bookmarkList,
  						Shortest[__ ~~ "RuleDelayed[\"" ~~ x__ ~~ "\", List[" ~~ y__ ~~ "]]]" ~~ __] :> {x, y <> "]"}
  						(*the last "]", cut by StringCases, needs to be added back manually*)
  					];
		 If[ Length[bookmarks] == 0,
        	Return[],
        	{bookmarkNames, bookmarkSettings} = Transpose[bookmarks];
        ];
        
		(*delete any additional \" in the names*)
		bookmarkNames = StringReplace[bookmarkNames, "\\\"" -> ""];
		
		(*get parameter names and their values; 
		they appear as Set[FE`parameterName$$someNumber = parameterValue];
        return in the form "parameterName = parameterValue"*)
		bookmarkSettings = 
		 	StringCases[#, 
		    	 Shortest["Set[FE`" ~~ x__ ~~ "$$" ~~ ___ ~~ ", " ~~ y__ ~~ "]"] :> {x <> " = " <> y}][[All, 1]] 
		     & /@ bookmarkSettings;
		
		(*2D controls come as List[x,y; replace with {x,y}*)
		bookmarkSettings = 
		 	StringReplace[#, {"List[" ~~ x__ :> "{" <> x <> "}"}] & /@ 
		  	bookmarkSettings;
		
		Return[{bookmarkNames, bookmarkSettings}];
     	
    ];

(*
BookmarksExport: write all created bookmarks into separate .dat files

input:	bookmarkList	.. list of bookmarks in FullForm, converted to string: ToString@FullForm@Typeset`bookmarkList$$

output: .dat files with all Manipulate settings, one for each bookmark (selected by user through dialog)
*)
ClearAll[BookmarksExport]
BookmarksExport[bookmarkList_] :=
    Block[ {bookmarks, bookmarkNames, bookmarkParameters, bookmarkSaveSwitches, userInput, switches, names, exportPath, exportFolder, exportFileName},
        (*get a list of parameter names and a list of values*)
        bookmarks = BookmarkContents[bookmarkList];
            
        (*if there were no bookmarks in the list, warn the user and exit*)
        If[ Length[bookmarks] == 0,
            DialogInput[
              DialogNotebook[{TextCell["No bookmarks to export - create one first!\n 
              \n- Click \[CirclePlus] in to top-right corner
              \n- Select \"Add To Bookmarks\"
              \n- Enter a name of the bookmark
              \n- Confirm by clicking \"add\""], 
                DefaultButton[]}]];
            Return[];
        ];
        
        (*split the list into names and parametr values*)
        {bookmarkNames, bookmarkParameters} = bookmarks;
        
                
        (*check which bookmarks already appear in the "bookmarks" folder; set saveSwitch True for those that do not and select only those*)
        SetDirectory[FileNameJoin[{NotebookDirectory[], "bookmarks"}]];
        bookmarkSaveSwitches = FileNames["FMSfit-bookmark_" ~~ # ~~ ".dat"] == {} & /@bookmarkNames;
        bookmarkNames = Pick[bookmarkNames, bookmarkSaveSwitches];
        bookmarkParameters = Pick[bookmarkParameters, bookmarkSaveSwitches];
        bookmarkSaveSwitches = Pick[bookmarkSaveSwitches, bookmarkSaveSwitches];
        
        (*allow user to check which bookmarks to export and to change their names*)
        userInput = DialogInput[{switches = bookmarkSaveSwitches, names = bookmarkNames},
 			Column[
 				{"Select bookmarks to save: ",
   				Grid[
   					Table[With[ {i = i},
              			{Checkbox[Dynamic@switches[[i]]], "  FMSfit-bookmark_", InputField[Dynamic@names[[i]], String, FieldSize -> 10], ".dat"}
              			],{i, Length[switches]}]
    			, Spacings -> {0,1}],   
   				
   				Row[{
   					Button["Save", DialogReturn[{switches, names}], ImageSize -> {100,25}],
   					CancelButton[ ImageSize -> {100,25}]
   				},"  "]
   
   			}, Spacings -> 2]
   		];
   		If[
   			ToString[userInput] == "$Canceled", 
   			Return[];
   			,
   			{bookmarkSaveSwitches, bookmarkNames} = userInput;
   			If[Count[bookmarkSaveSwitches,True] == 0, Return[];];
   		];
   		
   		(*select bookmarks that were chosen to be exported*)
   		bookmarkNames = Pick[bookmarkNames, bookmarkSaveSwitches];
   		bookmarkParameters = Pick[bookmarkParameters, bookmarkSaveSwitches];
   			
   		(*export*)
   		MapThread[
   			Export["FMSfit-bookmark_" <> #1 <> ".dat", Column[#2]]&,
   			{bookmarkNames, bookmarkParameters}
   		];
   		       
    ];

(*
BookmarkImport: import a bookmark from a file

input:	fileName	.. full path to the bookmark file; not necessary

output:	bookmark	.. a string "bookmarkName":>{par1name = par1value, ...}; in Manipulate used as AppendTo[Typeset`bookmarkList$$, ImportString[bookmark, "Package"]]
*)
ClearAll[BookmarkImport]
BookmarkImport[fileName_: "Null"] :=
    Block[ {bookmarkFileName, bookmarkParameters, bookmarkName, bookmark},
        (*if no fileName is given as argument, ask for it*)
        bookmarkFileName = If[ fileName == "Null",
                               SystemDialogInput["FileOpen", 
                                FileNameJoin[{NotebookDirectory[], "bookmarks", "FMSfit-bookmark_sample.dat"}]],
                               fileName
                           ];
        If[ ToString@bookmarkFileName == "$Canceled",
            Return[]
        ];
        (*read the file by lines*)
        bookmarkParameters = ReadList[bookmarkFileName, String];
        (*in each line, replace parameterName by $CellContext`parameterName$$*)
        bookmarkParameters = 
	         StringReplace[#, 
	            x__ ~~ " = " ~~ y__ :> "$CellContext`" <> x <> "$$ = " <> y 
	            ] & /@ bookmarkParameters;
        (*convert the list of parameters into a string*)
        bookmarkParameters = ToString[bookmarkParameters];
        (*extract the name of the bookmark from filename - whatever is between a _ and .dat*)
        bookmarkFileName = FileNameTake[bookmarkFileName, -1];
        bookmarkName = StringCases[bookmarkFileName,"FMSfit-bookmark_" ~~ x__ ~~ ".dat" -> x];
        bookmarkName = If[bookmarkName =={},(*if the filename is of a different form, just keep the filename*)
        					StringDrop[bookmarkFileName,-4],
        					bookmarkName[[1]]
        				];
	         
        (*final bookmark form: "bookmarkName":>{...}*)
        bookmark = 
         "\"" <> bookmarkName <> "\"" <> ":>" <> bookmarkParameters;
        Return[bookmark]
    ];
