function []=CalcTalStat(DegOfFree,TdFilePath,ValFilePath,MaxRange);

%	CalcTalStat.m

%   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.  

  global BATValVec BAChar BACoorChar BAStrucChar RangeVec
%  BATValVec=[]; BACoorChar=[]; BAStrucChar=[];RangeVec=[];
%
%=================================================================
if nargin<4; MaxRange=3; end
if nargin<3; ValFilePath=[]; end
if nargin<2; TdFilePath=[]; end
if nargin<1; DegOfFree=21; end


MinNVox=2;

if isempty(DegOfFree)
    MaxRange=MaxRange
    DegOfFree=DegOfFree
    pause(2)
end
if isempty(TdFilePath)
    DefFolder=['D:\ValAro\00_Stats_NoScaling\RandStatsSmooth08\Male\u010k2\Filtered\018_AroExp1HiGTLoAll']
    Mask=['TD*.*'];
    TdFilePath=spm_get(1,Mask,'Choose a Talairach Daemon file',DefFolder);
end
if isempty(ValFilePath)
    [TdFile,TdPath]=SepFilePath(TdFilePath)
    ValFile=TdFile(3:length(TdFile)-4);
    ValFile=[ValFile,'Val.txt'];
    ValFilePath=[TdPath,ValFile];
    if ~exist(ValFilePath)
        Mask=['*Val.txt'];
        ValFilePath=spm_get(1,Mask,'Choose a t-value text file',DefFolder)
        if isempty(ValFilePath); return; end
    end
end
fidTD=fopen(TdFilePath,'r');
fidVal=fopen(ValFilePath,'r');
LineInd=0;
TDFileStatus=1;
TValVec=[];
while 1
    LineInd=LineInd+1;
	TdLine = fgetl(fidTD);
    if isempty(findstr(TdLine,'Record'))
        ValLine= fgetl(fidVal);
	    if ~isstr(TdLine)
		    break
	    end
	    if ~isempty(findstr(TdLine,'Range='))
		    TDFileStatus=2;
	    end
        TValVec(LineInd)=str2num(ValLine);
    end
end
fclose(fidVal);
fseek(fidTD,0,'bof');
TdLine = fgetl(fidTD);%First line
Ind=0;
BAInd=0;
CerebellumInd=0;
ThalamusInd=0;
HypothalamusInd=0;
AmygdalaInd=0;
CaudateInd=0;
HippocampusInd=0;
PutamenInd=0;
PallidusInd=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');
            if Ind==0
                fprintf('No BA voxel',Ind);
                return;
            else
                fprintf('%g total voxel ...  ',Ind);
                break;
            end
        end
        if ~isempty(findstr(TdLine,'Brodmann')) | ...
           ~isempty(findstr(TdLine,'Cerebellum')) | ...
           ~isempty(findstr(TdLine,'Thalamus')) | ...
           ~isempty(findstr(TdLine,'Hypothalamus')) | ...
           ~isempty(findstr(TdLine,'Amygdala')) | ...
           ~isempty(findstr(TdLine,'Caudate')) | ...
           ~isempty(findstr(TdLine,'Hippocampus')) | ...
           ~isempty(findstr(TdLine,'Putamen')) | ...
           ~isempty(findstr(TdLine,'Pallidus'))
            Ind=Ind+1;
            
            Tmp=findstr(TdLine,',');
            TdCoorTmp=TdLine(Tmp(1)+2:Tmp(4)-1);
            if ~isempty(findstr(TdLine,'Brodmann'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+15:length(TdLine)-1)];
                BAInd=BAInd+1;
                if round(BAInd./25)==BAInd./25
                    fprintf('Found %g BA voxel ...  \n',BAInd);
                end
            elseif ~isempty(findstr(TdLine,'Cerebellum'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(9)-1)];
                CerebellumInd=CerebellumInd+1;
                if round(CerebellumInd./25)==CerebellumInd./25
                    fprintf('Found %g Cerebellum voxel ...  \n',CerebellumInd);
                end
            elseif ~isempty(findstr(TdLine,'Thalamus'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                ThalamusInd=ThalamusInd+1;
                fprintf('Found %g Thalamus voxel ...  \n',ThalamusInd);
            elseif ~isempty(findstr(TdLine,'Amygdala'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                AmygdalaInd=AmygdalaInd+1;
                fprintf('Found %g Amygdala voxel ...  \n',AmygdalaInd);
            elseif ~isempty(findstr(TdLine,'Caudate'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                CaudateInd=CaudateInd+1;
                fprintf('Found %g Caudate voxel ...  \n',CaudateInd);
            elseif ~isempty(findstr(TdLine,'Hippocampus'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                HippocampusInd=HippocampusInd+1;
                fprintf('Found %g Hippocampus voxel ...  \n',HippocampusInd);
            elseif ~isempty(findstr(TdLine,'Putamen'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                PutamenInd=PutamenInd+1;
                fprintf('Found %g Putamen voxel ...  \n',PutamenInd);
            elseif ~isempty(findstr(TdLine,'Pallidus'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                PallidusInd=PallidusInd+1;
                fprintf('Found %g Pallidus voxel ...  \n',PallidusInd);
            elseif ~isempty(findstr(TdLine,'Hypothalamus'))
                TdLineTmp=[TdLine(Tmp(4)+2:Tmp(7)+4),',',TdLine(Tmp(8)+1:length(TdLine)-1)];
                HypothalamusInd=HypothalamusInd+1;
                fprintf('Found %g Hypothalamus voxel ...  \n',HypothalamusInd);
            end
			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
    if Ind==0
        BATValVec=[];
    else
	    BATValVec=TValVec(BAValIndVec);
    end
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');
            if Ind==0
			    fprintf('No BA voxel',Ind);
                return;
            else
                 fprintf('%g BA voxel ...  ',Ind);
                 break;
             end
		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
    if Ind==0
        BATValVec=[];
        BATValVec=[];
    else
        BATValVec=TValVec(BAValIndVec);
        %[MaxBATValVec,MaxBATValVecInd]=max(BATValVec)
        %BAChar(MaxBATValVecInd,:)
        %BACoorChar(MaxBATValVecInd,:)
        %pause
        MaxRangeIndVec=find(RangeVec<=MaxRange);
        BAChar=BAChar(MaxRangeIndVec,:);
        BACoorChar=BACoorChar(MaxRangeIndVec,:);
        BATValVec=BATValVec(MaxRangeIndVec);
    end
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);
   MaxBACoorChar=BACoorCharTmp(VoxInStrucIndVec(MaxVoxTInd),:);
   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(round(MaxVoxTVal.*100)./100),',',num2str(round(MeanTValInStrucVec(i).*100)./100)];
	TmpLine3=[',',num2str(round(MaxVoxZVal.*100)./100),',',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,:);
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(12)-1)];
		else
			Line=[Line(1:5),Line(Tmp(2)+1:Tmp(3)-1),Line(Tmp(4):Tmp(12)-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(12)-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(12)-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

