function [ima_fil,final_estimate]=denoise_clusters(ima_patchs_vect,func_denoising_patches,func_thresholding,func_recontruction,IDX,Patch_width,nb_axis,nb_clusters,m,n,newton,parallel)
% denoise_clusters computes the denoising operation cluster by cluster
%
%   INPUT:
%     ima_patchs_vect        : collection of patches (matrix)
%     func_denoising_patches : choice of the denoising perform
%                                (gausiann or poisson are handle)
%     Patch_width            : width of the square patches to consider
%     func_thresholding      : only no thresholding is yet supported
%  	  func_recontruction     : build the representaion of the patches
%                              thanks to the model used (gausiann or
%                              poisson are handle)
%     IDX                    : indexes of the cluster
%  	  nb_clusters		     : number of clusters choosen
%  	  nb_axis		         : number of axis/components choosen
%     Patch_width            : width of the square patches to consider
%     m,n                    : size of the original image
%   OUTPUT:
%     ima_fil                : final denoised image
%     final_estimate         : collection of denoised patches
%
%  [ima_fil,final_estimate]=denoise_clusters(ima_patchs_vect,...
%                           func_denoising_patches,func_thresholding,...
%                           func_recontruction,IDX,Patch_width,nb_axis,...
%                           nb_clusters,m,n)
%   Perform the denoising of the collection of patches and the reprojection
%   to get an estimation of the original image (ima_fil). The function can
%   handle both poisson and gaussian (after Anscombe transform) case.
%
%   Copyright (C) 2012 NL-PCA project
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)

%---------------------------------------------------------------------
%
%   This file is part of NL-PCA.
%
%   NL-PCA is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as
%   published by the Free Software Foundation, either version 3 of
%   the License, or (at your option) any later version.
%
%   NL-PCA is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public
%   License along with NL-PCA.  If not, see
%   <http://www.gnu.org/licenses/>.
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)

final_estimate=zeros(size(ima_patchs_vect));
size_cluster = zeros(nb_clusters,1);
indexes=cell(nb_clusters,1);

switch parallel
	case 1 % Use parallelization
		parfor i = 1:nb_clusters
			indexes{i} = find(IDX == i); % finds the location in IDX for a specific cluster number
			size_cluster(i) = size((indexes{i}),1); % number of patches in each cluster
		end
	case 0 % No parallelization
		for i = 1:nb_clusters
			indexes{i} = find(IDX == i); % finds the location in IDX for a specific cluster number
			size_cluster(i) = size((indexes{i}),1); % number of patches in each cluster
		end
end

switch parallel
	case 1 % Use parallelization
		
		temp_clusters = cell(1,nb_clusters);
		
		for k = 1:(nb_clusters)
			% this finds all the patches in a specific cluster and inserts
			% it into 'temp_clusters'
			temp_clusters(1,k) = {ima_patchs_vect(indexes{k},:)};
		end
		
		temp_final_est = cell(1,nb_clusters);
		ima_ppca_cluster = cell(1,nb_clusters);
		ima_patchs_fil_cluster = cell(1,nb_clusters);
		
		parfor p = 1:(nb_clusters)
			
			%temp_clusters = ima_patchs_vect(indexes{k},:); % takes the patches of the kth cluster of ima_patchs_vect
			
			[utmp1,stmp1,vtmp1] = svd_new((temp_clusters{1,p}),nb_axis); % est 1
			est1 = utmp1*stmp1*vtmp1';
			est1 = est1.*(est1>0); % sets negative elements of est1 to 0
			
			[utmp2,stmp2,vtmp2] = svd_new(log((temp_clusters{1,p})+1),nb_axis); % est 2
			est2 = exp(utmp2*stmp2*vtmp2');
			
			[utmp3,stmp3,vtmp3] = svd_new(log(est1+1),nb_axis); % est 3
			est3 = exp(utmp3*stmp3*vtmp3');
			
			negLogLike1 = negLogLike((temp_clusters{1,p}),est1);
			negLogLike2 = negLogLike((temp_clusters{1,p}),est2);
			negLogLike3 = negLogLike((temp_clusters{1,p}),est3);
			[~,bestEst] = min([negLogLike1,negLogLike2,negLogLike3]); % compares est1, est2, est3
			
			switch bestEst % fills final_estimate with the best estimate for those values
				case 1
					temp_final_est(1,p) = {est1};
				case 2
					temp_final_est(1,p) = {est2};
				case 3
					temp_final_est(1,p) = {est3};
			end
			
			if newton % SPIRALTAP % Determines if SPIRALTAP is used
				
				if negLogLike2 < negLogLike3
					utmp = utmp2;
					stmp = stmp2;
					vtmp = vtmp2;
				else
					utmp = utmp3;
					stmp = stmp3;
					vtmp = vtmp3;
				end
				
				startv = vtmp'; % finds starting matrices
				startu = utmp*stmp;
				
				ima_ppca_cluster{1,p} = func_denoising_patches({temp_clusters{1,p};startu;startv}); % est 4
				ima_patchs_fil_cluster{1,p} = func_recontruction(ima_ppca_cluster{1,p});
				
				if negLogLike(temp_clusters{1,p},ima_patchs_fil_cluster{1,p}) < negLogLike(temp_clusters{1,p},temp_final_est{1,p})
					% compares the best of est(1,2,3) with est(4)
					temp_final_est(1,p) = {ima_patchs_fil_cluster{1,p}};
					disp(strcat('Cluster  ',num2str(p),': Newton iteration did improve subspace estimate'))
				else
					disp(strcat('Cluster  ',num2str(p),': Newton iteration did not improve subspace estimate'))
				end
				
				ima_ppca_cluster{1,p} = []; % clears the cell at cluster k
				ima_patchs_fil_cluster{1,p} = []; % clears the cell at cluster k
			end
		end
		
		for q = 1:nb_clusters % fills 'final_estimate' with the denoised clusters
			final_estimate(indexes{q},:) = temp_final_est{1,q};
		end
		
		ima_fil = reprojection_UWA(reshape(final_estimate,[(m-Patch_width+1),(n-Patch_width+1),Patch_width^2]));
		
	case 0 % No parallelization
		
		ima_ppca_cluster = cell(1,nb_clusters);
		ima_patchs_fil_cluster = cell(1,nb_clusters);
		
		for k = 1:(nb_clusters)
			
			temp_clusters = ima_patchs_vect(indexes{k},:); % takes the patches of the kth cluster of ima_patchs_vect
			
			[utmp1,stmp1,vtmp1] = svd_new(temp_clusters,nb_axis); % est 1
			est1 = utmp1*stmp1*vtmp1';
			est1 = est1.*(est1>0); % sets negative elements of est1 to 0
			
			[utmp2,stmp2,vtmp2] = svd_new(log(temp_clusters+1),nb_axis); % est 2
			est2 = exp(utmp2*stmp2*vtmp2');
			
			[utmp3,stmp3,vtmp3] = svd_new(log(est1+1),nb_axis); % est 3
			est3 = exp(utmp3*stmp3*vtmp3');
			
			negLogLike1 = negLogLike(temp_clusters,est1);
			negLogLike2 = negLogLike(temp_clusters,est2);
			negLogLike3 = negLogLike(temp_clusters,est3);
			[~,bestEst] = min([negLogLike1,negLogLike2,negLogLike3]); % compares est1, est2, est3
			
			switch bestEst % fills final_estimate with the best estimate for those values
				case 1
					final_estimate(indexes{k},:) = est1;
				case 2
					final_estimate(indexes{k},:) = est2;
				case 3
					final_estimate(indexes{k},:) = est3;
			end
			
			if newton % Determines if SPIRALTAP is used
				
				if negLogLike2 < negLogLike3
					utmp = utmp2;
					stmp = stmp2;
					vtmp = vtmp2;
				else
					utmp = utmp3;
					stmp = stmp3;
					vtmp = vtmp3;
				end
				
				startv = vtmp'; % finds starting matrices
				startu = utmp*stmp;
				
				ima_ppca_cluster{k} = func_denoising_patches({ima_patchs_vect(indexes{k},:);startu;startv}); % est 4
				ima_patchs_fil_cluster{k} = func_recontruction(ima_ppca_cluster{k});
				
				if negLogLike(temp_clusters,ima_patchs_fil_cluster{k}) < negLogLike(temp_clusters,final_estimate(indexes{k},:))
					% compares the best of est(1,2,3) with est(4)
					final_estimate(indexes{k},:) = ima_patchs_fil_cluster{k};
					disp(strcat('Cluster ',num2str(k),': Newton iteration did improve subspace estimate'))
				else
					disp(strcat('Cluster ',num2str(k),': Newton iteration did not improve subspace estimate'))
				end
				
				ima_ppca_cluster{k} = []; % clears the cell at cluster k
				ima_patchs_fil_cluster{k} = []; % clears the cell at cluster k
			end
		end
		ima_fil = reprojection_UWA(reshape(final_estimate,[(m-Patch_width+1),(n-Patch_width+1),Patch_width^2]));
end

