function[]=CalcTalStat(TdFilePath,ValFilePath,MaxRange,DegOfFree);
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
%==================================================================
%
%	CalcTalStat.m
%
%	Markus Junghoefer	[2002]
%
%
%	Function definition
%
 global BATValVec BAChar BACoorChar BAStrucChar RangeVec
 BATValVec=[]; BACoorChar=[]; BAStrucChar=[];RangeVec=[];
%
%=================================================================
if nargin<4; DegOfFree=17; end
if nargin<3; MaxRange=3; end
if nargin<2; ValFilePath=[]; end
if nargin<1; TdFilePath=[]; end

DefFilePath=SetDefPath(1,'TD*.txt');
[TdFile,TdPath,TdFilePath]=ReadFilePath(TdFilePath,DefFilePath,'Choose a TD file:');
if TdFile==0; return; end
SetDefPath(2,TdPath);
DefFilePath=SetDefPath(1,'*Val.txt');
[ValFile,ValPath,ValFilePath]=ReadFilePath(ValFilePath,DefFilePath,'Choose a Val file:');
if ValFile==0; return; end
SetDefPath(2,ValPath);

fidTD=fopen(TdFilePath,'r');
fidVal=fopen(ValFilePath,'r');
LineInd=0;
TDFileStatus=1;
TValVec=[];
while 1
    LineInd=LineInd+1;
	TdLine = fgetl(fidTD);
    ValLine= fgetl(fidVal);
	if ~isstr(TdLine)
		break
	end
	if ~isempty(findstr(TdLine,'Range='))
		TDFileStatus=2;
	end
    TValVec(LineInd)=str2num(ValLine);
end
fclose(fidVal);
fseek(fidTD,0,'bof');
Ind=0;
LineInd=0;
if TDFileStatus==1
	fprintf(1,'Talairach Daemon file status: Simple structure\n\n')
	while 1
		TdLine = fgetl(fidTD);
		LineInd=LineInd+1;
		if ~isstr(TdLine)
			fprintf('\n\n');
			fprintf('Found %g BA voxel ...  ',Ind);
			break
		end
		if ~isempty(findstr(TdLine,'Brodmann'))
			Ind=Ind+1;
			if round(Ind./25)==Ind./25
				fprintf('\n');
				fprintf('Found %g BA voxel ...  ',Ind);
			end
			Tmp=findstr(TdLine,',');
			TdCoorTmp=TdLine(Tmp(1)+2:Tmp(4)-1);
			TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+15:length(TdLine))];
			LTdLineVec(Ind)=length(TdLineTmp);
			if Ind==1
				BAChar=char(TdLineTmp);
				BACoorChar=char(TdCoorTmp);
			else
				BAChar=char(BAChar,TdLineTmp);
				BACoorChar=char(BACoorChar,TdCoorTmp);
			end
			BAValIndVec(Ind)=LineInd;
		end
	end
	BATValVec=TValVec(BAValIndVec);
elseif TDFileStatus==2
	fprintf(1,'Talairach Daemon file status: Nearest gray matter\n\n')
	while 1
		TdLine = fgetl(fidTD);
		LineInd=LineInd+1;
		if ~isstr(TdLine)
			fprintf('\n\n');
			fprintf('Found %g BA voxel ...  ',Ind);
			break
		end
		if isempty(findstr(TdLine,'No GM'))
			Ind=Ind+1;
			if round(Ind./25)==Ind./25
				fprintf('\n');
				fprintf('Found %g BA voxel ...  ',Ind);
			end
			Tmp=findstr(TdLine,',');
			TdCoorTmp=TdLine(Tmp(1)+2:Tmp(4)-1);
			TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+15:Tmp(9)-1)];
			RangeVec(Ind)=str2num(TdLine(Tmp(9)+7:length(TdLine)));
			LTdLineVec(Ind)=length(TdLineTmp);
			if Ind==1
				BAChar=char(TdLineTmp);
				BACoorChar=char(TdCoorTmp);
			else
				BAChar=char(BAChar,TdLineTmp);
				BACoorChar=char(BACoorChar,TdCoorTmp);
			end
			BAValIndVec(Ind)=LineInd;
		end
	end
	BATValVec=TValVec(BAValIndVec);
	MaxRangeIndVec=find(RangeVec<=MaxRange);
	BAChar=BAChar(MaxRangeIndVec,:);
	BACoorChar=BACoorChar(MaxRangeIndVec,:);
	BATValVec=BATValVec(MaxRangeIndVec);
end
fclose(fidTD);

BAZValVec=spm_t2z(BATValVec,DegOfFree,10^-16);

BATValTmp=BATValVec;
BAZValTmp=BAZValVec;

i=0;
BACharTmp=BAChar;
BACoorCharTmp=BACoorChar;
BAVec=[];
NRightStruc=0;
NLeftStruc=0;
RightBAStrucChar=[];
LeftBAStrucChar=[];
while ~isempty(BACharTmp)
	i=i+1;
	Line=deblank(BACharTmp(1,:));
	VoxInStrucIndVec=strmatch(Line,BACharTmp);
	NVoxInStruc=length(VoxInStrucIndVec);
	NVoxInStrucVec(i)=NVoxInStruc;
	VoxTValVec=BATValTmp(VoxInStrucIndVec);
	VoxZValVec=BAZValTmp(VoxInStrucIndVec);
	[MaxVoxTVal,MaxVoxTInd]=max(VoxTValVec);
	[MaxVoxZVal,MaxVoxZInd]=max(VoxZValVec);
	SumVoxTValVec=sum(VoxTValVec);
	SumVoxZValVec=sum(VoxZValVec);
	MeanTValInStrucVec(i)=mean(VoxTValVec);
	MeanZValInStrucVec(i)=mean(VoxZValVec);
	CoorMat=zeros(NVoxInStruc,3);
	for m=1:NVoxInStruc
		CoorLine=deblank(BACoorCharTmp(m,:));
		Tmp=findstr(CoorLine,',');
		Coor1=CoorLine(1:Tmp(1)-1);
		Coor2=CoorLine(Tmp(1)+2:Tmp(2)-1);
		Coor3=CoorLine(Tmp(2)+2:length(CoorLine));
		CoorMat(m,:)=[str2num(Coor1),str2num(Coor2),str2num(Coor3)];
	end
	if NVoxInStruc>1
		MeanCoorVec=round(mean(CoorMat));
	else
		MeanCoorVec=round(CoorMat);
	end
	MeanTWeightCoorVec=MeanCoorVec;
	MeanZWeightCoorVec=MeanCoorVec;
	CoorMatTTmp=CoorMat;
	CoorMatZTmp=CoorMat;
	for j=1:3
		CoorMatTTmp(:,j)=CoorMatTTmp(:,j).*VoxTValVec';
		CoorMatZTmp(:,j)=CoorMatZTmp(:,j).*VoxZValVec';
		MeanTWeightCoorVec(j)=round(sum(CoorMatTTmp(:,j))./SumVoxTValVec);
		MeanZWeightCoorVec(j)=round(sum(CoorMatZTmp(:,j))./SumVoxZValVec);
	end
	TmpLine1=[',',num2str(NVoxInStrucVec(i))];
	TmpLine2=[',',num2str(MaxVoxTVal),',',num2str(round(MeanTValInStrucVec(i).*100)./100)];
	TmpLine3=[',',num2str(MaxVoxZVal),',',num2str(round(MeanZValInStrucVec(i).*100)./100)];
	TmpLine4=[',',num2str(round(CoorMat(MaxVoxTInd,:)))];
	TmpLine5=[',',num2str(MeanCoorVec)];
	TmpLine6=[',',num2str(MeanTWeightCoorVec)];
	TmpLine7=[',',num2str(MeanZWeightCoorVec)];
	
	if ~isempty(findstr(Line,'Right')) 
		NRightStruc=NRightStruc+1;
		NRightVoxInStrucVec(NRightStruc)=NVoxInStruc;
		RightBAStrucChar=char(RightBAStrucChar,[Line,TmpLine1,TmpLine2,TmpLine3,TmpLine4,TmpLine5,TmpLine6,TmpLine7]);
	else
		NLeftStruc=NLeftStruc+1;
		NLeftVoxInStrucVec(NLeftStruc)=NVoxInStruc;
		LeftBAStrucChar=char(LeftBAStrucChar,[Line,TmpLine1,TmpLine2,TmpLine3,TmpLine4,TmpLine5,TmpLine6,TmpLine7]);
	end
	NRestLines=size(BACharTmp,1);
	l=0;
	RestIndVec=[];
	for k=1:NRestLines
		if all(VoxInStrucIndVec-k)
			l=l+1;
			RestIndVec(l)=k;
		end
	end
	BATValTmp=BATValTmp(RestIndVec);
    BAZValTmp=BAZValTmp(RestIndVec);
	BACharTmp=BACharTmp(RestIndVec,:);
	BACoorCharTmp=BACoorCharTmp(RestIndVec,:);
end	
if NRightStruc>0
	RightBAStrucChar=RightBAStrucChar(2:size(RightBAStrucChar,1),:);
end
if NLeftStruc>0
	LeftBAStrucChar=LeftBAStrucChar(2:size(LeftBAStrucChar,1),:);
end
[NonUse,SortNRightVoxInStrucVec]=sort(-NRightVoxInStrucVec);
RightBAStrucChar=RightBAStrucChar(SortNRightVoxInStrucVec,:);
[NonUse,SortNLeftVoxInStrucVec]=sort(-NLeftVoxInStrucVec);
LeftBAStrucChar=LeftBAStrucChar(SortNLeftVoxInStrucVec,:);
MinNVox=5;
NBilatStruc=0;
NJustRightStruc=0;
NJustLeftStruc=0;
BilatBAStrucChar=[];
JustRightBAStrucChar=[];
JustLeftBAStrucChar=[];
for i=1:NRightStruc	
	RightLine=RightBAStrucChar(i,:);
	Tmp=findstr(RightLine,',');
	RightLineTmp=deblank(RightLine(Tmp(1)+1:Tmp(5)-1));
	BilatStatus=0;
	j=0;
	while j<NLeftStruc
		j=j+1;
		LeftLine=LeftBAStrucChar(j,:);
		if ~isempty(findstr(LeftLine,RightLineTmp))
			BilatStatus=1;
			NBilatStruc=NBilatStruc+1;
			LeftBilatVec(NBilatStruc)=j;
			BilatBAStrucChar=char(BilatBAStrucChar,RightLine);
			BilatBAStrucChar=char(BilatBAStrucChar,LeftLine);
			break
		end
	end
	if ~BilatStatus
		NJustRightStruc=NJustRightStruc+1;
		JustRightBAStrucChar=char(JustRightBAStrucChar,RightLine);
	end
end
for j=1:NLeftStruc
	if all(LeftBilatVec-j)
		NJustLeftStruc=NJustLeftStruc+1;
		JustLeftBAStrucChar=char(JustLeftBAStrucChar,LeftBAStrucChar(j,:));
	end
end
ShortBilatBAStrucChar=[];
if NBilatStruc>0
	BilatBAStrucChar=BilatBAStrucChar(2:size(BilatBAStrucChar,1),:);
	for i=1:size(BilatBAStrucChar,1)
		Line=deblank(BilatBAStrucChar(i,:));
		Tmp=findstr(Line,',');
		NBilatVox=Line(Tmp(5)+1:Tmp(6)-1);
		NBilatVoxVec(i)=str2num(NBilatVox);
		if ~isempty(findstr(Line,'Right'))
			Line=[Line(1:5),' ',Line(Tmp(2)+1:Tmp(3)-1),Line(Tmp(4):Tmp(10)-1)];
		else
			Line=[Line(1:5),Line(Tmp(2)+1:Tmp(3)-1),Line(Tmp(4):Tmp(10)-1)];
		end
		ShortBilatBAStrucChar=char(ShortBilatBAStrucChar,Line);
	end
	ShortBilatBAStrucChar=ShortBilatBAStrucChar(2:size(ShortBilatBAStrucChar,1),:);
	NMinVoxIndVec=find(NBilatVoxVec>=MinNVox);
	MinBilatBAStrucChar=BilatBAStrucChar(NMinVoxIndVec,:);
	MinShortBilatBAStrucChar=ShortBilatBAStrucChar(NMinVoxIndVec,:);
end
ShortJustRightBAStrucChar=[];
if NJustRightStruc>0
	JustRightBAStrucChar=JustRightBAStrucChar(2:size(JustRightBAStrucChar,1),:);
	for i=1:size(JustRightBAStrucChar,1)
		Line=deblank(JustRightBAStrucChar(i,:));
		Tmp=findstr(Line,',');
		NJustRightVox=Line(Tmp(5)+1:Tmp(6)-1);
		NJustRightVoxVec(i)=str2num(NJustRightVox);
		Line=[Line(Tmp(2)+1:Tmp(3)-1),Line(Tmp(4):Tmp(10)-1)];
		ShortJustRightBAStrucChar=char(ShortJustRightBAStrucChar,Line);
	end
	ShortJustRightBAStrucChar=ShortJustRightBAStrucChar(2:size(ShortJustRightBAStrucChar,1),:);
	NMinVoxIndVec=find(NJustRightVoxVec>=MinNVox);
	MinJustRightBAStrucChar=JustRightBAStrucChar(NMinVoxIndVec,:);
	MinShortJustRightBAStrucChar=ShortJustRightBAStrucChar(NMinVoxIndVec,:);
end
ShortJustLeftBAStrucChar=[];
if NJustLeftStruc>0
	JustLeftBAStrucChar=JustLeftBAStrucChar(2:size(JustLeftBAStrucChar,1),:);
	for i=1:size(JustLeftBAStrucChar,1)
		Line=deblank(JustLeftBAStrucChar(i,:));
		Tmp=findstr(Line,',');
		NJustLeftVox=Line(Tmp(5)+1:Tmp(6)-1);
		NJustLeftVoxVec(i)=str2num(NJustLeftVox);
		Line=[Line(Tmp(2)+1:Tmp(3)-1),Line(Tmp(4):Tmp(10)-1)];
		ShortJustLeftBAStrucChar=char(ShortJustLeftBAStrucChar,Line);
	end
	ShortJustLeftBAStrucChar=ShortJustLeftBAStrucChar(2:size(ShortJustLeftBAStrucChar,1),:);
	NMinVoxIndVec=find(NJustLeftVoxVec>=MinNVox);
	MinJustLeftBAStrucChar=JustLeftBAStrucChar(NMinVoxIndVec,:);
	MinShortJustLeftBAStrucChar=ShortJustLeftBAStrucChar(NMinVoxIndVec,:);
end
SortBAStrucChar=char('Bilateral',BilatBAStrucChar,'Right',JustRightBAStrucChar,'Left',JustLeftBAStrucChar)
MinSortBAStrucChar=char('Bilateral',MinBilatBAStrucChar,'Right',MinJustRightBAStrucChar,'Left',MinJustLeftBAStrucChar)
ShortSortBAStrucChar=char('Bilateral',MinShortBilatBAStrucChar,'Right',MinShortJustRightBAStrucChar,'Left',MinShortJustLeftBAStrucChar)
%==============================================================================
OutFilePath=[TdFilePath(1:length(TdFilePath)-4),'Sort.txt'];
fidOut=fopen(OutFilePath,'w');
fprintf(fidOut,'Right/Left,Lobe,Gyrus,Gray/White,BA,NVox,Max. t,Mean t,Max. Z,Mean Z,Max. xyz,Mean xyz,t-weighted xyz,Z-weighted xyz\n');
for i=1:size(SortBAStrucChar,1)
	Line=deblank(SortBAStrucChar(i,:));
	fprintf(fidOut,Line);
	fprintf(fidOut,'\n');
end
fclose(fidOut);

%==============================================================================
OutFilePath=[TdFilePath(1:length(TdFilePath)-4),'Min',int2str(MinNVox),'Sort.txt'];
fidOut=fopen(OutFilePath,'w');
fprintf(fidOut,'Right/Left,Lobe,Gyrus,Gray/White,BA,NVox,Max. t,Mean t,Max. Z,Mean Z,Max. xyz,Mean xyz,t-weighted xyz,Z-weighted xyz\n');
for i=1:size(MinSortBAStrucChar,1)
	Line=deblank(MinSortBAStrucChar(i,:));
	fprintf(fidOut,Line);
	fprintf(fidOut,'\n');
end
fclose(fidOut);
%==============================================================================
OutFilePath=[TdFilePath(1:length(TdFilePath)-4),'ShortMin',int2str(MinNVox),'Sort.txt'];
fidOut=fopen(OutFilePath,'w');
fprintf(fidOut,'Gyrus,BA,NVox,Max. t,Mean t,Max. Z,Mean Z,Max. xyz,Mean xyz\n');
for i=1:size(ShortSortBAStrucChar,1)
	Line=deblank(ShortSortBAStrucChar(i,:));
	fprintf(fidOut,Line);
	fprintf(fidOut,'\n');
end
fclose(fidOut);
return
