% Matlab script for computing the Aurofacial Asymmetry as predicted by the Axial Twist hypothesis. % The first draft was created by Niko Troje (2016). % Changes were made by Marc de Lussanet (2018) % The script needs the symmetric face and the individual faces from the Faces database, and the Facerenderer functions that come with it. Access can be obtained by contacting http://faces.kyb.tuebingen.mpg.de addpath ./Faces ./Faces/MAT %% flags: CalculateNewFaces = 0; % set to 1 to compute the new rendered faces CalculateAverages = 0; % set to 1 to compute the averages of all faces % The reference head seems to be perfectly symmetric, but rotated a bit. % These data structures have 6 columns. The first three encode the x- % (left/right), y- (height) and z- (depth) coordinates of the 75337 mesh vertices. % The next three columns are RGB texture values. XX=1;YY=2;ZZ=3;RR=4;GG=5;BB=6;XYZ=XX:ZZ; %% Load reference head (both sides are the average left) Tmp = load('Tex200L'); % der referenz-kopf aus zwei linke haelfte head200L = Tmp.head - [-104.1387 10455 17270.5822 0 0 0]; % subtract offset with respect to ear Tragi symhd200L= rotatehead(head200L,'y',-0.34585); % rotate to symmetry (degrees) %% Sort the vertices such that corresponding ones on the left and right side are kept together. % in the range of the eyes, many vertices have almost the same ZZ value, so sort by abs(XX) in this range % 23985 43665 (kritische range voor de YY waardes) Sort = [symhd200L(:,YY),symhd200L(:,ZZ),abs(symhd200L(:,XX))]; CritIdx = find((Sort(:,1)>=23985) & (Sort(:,1)<=43665)); Sort(CritIdx,:) = [Sort(CritIdx,1),Sort(CritIdx,3),Sort(CritIdx,2)]; [~,Idx] = sortrows(Sort); sortedsym = symhd200L(Idx,XYZ); [~,rIdx] = sort(Idx); % YY : 0 615 1230 1845 2460 ... 11070/615 % now, the matched X-coordinates should be following each other d = zeros(1,length(sortedsym)-1)==0; for i=1:length(sortedsym)-1 d(i) = sign(sortedsym(i,XX))==sign(sortedsym(i+1,XX)); end f = find(d); % in most of these cases, the order of the left and right matched vertex is % reversed, so it is not found. Search them by comparing the distance dist = zeros(1,length(sortedsym)); for i=1:length(sortedsym)-1 dist(i) =(sortedsym(i,XX)+sortedsym(i+1,XX)).^2 + ... (sortedsym(i,YY)-sortedsym(i+1,YY)).^2 + ... (sortedsym(i,ZZ)-sortedsym(i+1,ZZ)).^2; end dist(end) = dist(end-1); dist = sqrt(dist); dist(2:end-1) = min(dist(2:end-1),dist(3:end)); g=find(dist>100); Idx2 = Idx; Idx2(g)=[]; [~,rIdx2] = sort(Idx2); % most of these are still matched, but I don't not how to find the pairs % diff(g) % figure; plot(1:length(dist),dist) % figure(); plot(abs(sortedsym(g(18:37),1)),sortedsym(g(18:37),3),'o'); fprintf('Es gibt %d nicht-gematchte vertices\n',length(g)); %% % We can assign correspondence to all the vertices % The vector contains for each vertex the index of the corresponding % vertex on the other side. Note that this correspondence holds for all % faces in the data base. correspondence = zeros(size(symhd200L,1),1); correspondence(Idx2(1:2:end)) = Idx2(2:2:end); correspondence(Idx2(2:2:end)) = Idx2(1:2:end); %% ================================================================ % For Marc's question, we may not even need interhemispheric correspondence. % We just look at differences between the perfectly symmetric head and the average head. %% rotate and translate such that the ear tragus are symmetric to the origin roty = -0.8 - 0.015; % de tweede zet de tragus symmetrisch rotz = 0.5 - 0.120; % de tweede zet de tragus symmetrisch Offset = [-311 10294 16661]; %% calculate/load average and symmetric faces for rendering if CalculateNewFaces load facenames %#ok<*UNRCH> hd200sym =zeros(75337,6); hd200mean =zeros(75337,6); for i=1:200 headi = load(names(i,:)); h = headi.head; h = rotatehead(h,'y',roty); h = rotatehead(h,'z',rotz); h(:,XYZ) = h(:,XYZ) - Offset; % h_flipped is going to be the mirror flipped version of h h_flipped = h; h_flipped(:,XX) = -h_flipped(:,XX); h_flipped(find(correspondence),:) = h(correspondence(find(correspondence)),:); h_flipped(:,XX) = -h_flipped(:,XX); % hsym is the mean of h and h_flipped, that is, the perfectly symmetric version. hsym = 0.5*(h+h_flipped); hd200sym = hd200sym + hsym /200; hd200mean = hd200mean + h /200; end save(fullfile('Faces','MAT','hd200sym')); save(fullfile('Faces','MAT','hd200mean')); else load hd200sym load hd200mean end % check: % figure; renderhead(rotatehead(hd200mean,'y',150)); title('the average of all faces'); % figure; renderhead(hd200mean); title('the average of all faces'); % figure; renderhead(rotatehead(hd200sym,'y',150)); title('the average of all faces'); % figure; renderhead(hd200sym); title('the average of all faces'); %% calculate angles of the matched pairs if CalculateAverages load facenames mhsym =zeros(75337,3); mhead =zeros(75337,3); angles =zeros(75337,200); radius =zeros(75337,200); for i=1:200 headi = load(names(i,:)); h = headi.head; h = rotatehead(h,'y',roty); h = rotatehead(h,'z',rotz); h = h(:,XYZ); % get rid of the colored texture h = h - Offset; % h_flipped is the mirror flipped version of h h_flipped =zeros(75337,3); h_flipped(find(correspondence),XYZ) = h(correspondence(find(correspondence)),XYZ); h_flipped(:,XX) = -h_flipped(:,XX); % hsym is the mean of h and h_flipped, that is, the perfectly symmetric version of the current face. h(correspondence==0,:) = 0; hsym = 0.5*(h+h_flipped); mhsym = mhsym + hsym/200; mhead = mhead + h/200; % angles is the angle of each voxel with the corresponding one of its symmetric version angles(:,i) = atan2(h(:,XX),h(:,ZZ)) - atan2(hsym(:,XX),hsym(:,ZZ)); %angles(:,i) = angles(:,i) - mean(angles(:,i)); radius(:,i) = sqrt(sum(h.^2,2)); end angles(angles> 10*pi/180) = nan; angles(angles<-10*pi/180) = nan; radius(~radius) = nan; figure; plot(angles); figure; plot(radius); save mhead; save mhsym; save angles; save radius; else load mhead; load mhsym; load angles; load radius; end %% Indices of the tragi % TrIx are the indices of the blobs to mark the tragi on the figures TrIx = find(mhsym(:,YY)>-1500 & mhsym(:,YY)<1500 & mhsym(:,ZZ)>-1500 & mhsym(:,ZZ)<1500); TrIx(mhsym(TrIx,1)==0) = []; % omit the non-matched ones in the center Tragus = mhsym(TrIx,:); % TrL and TrL are the single vortex indices of the left and right tragus [~, TrR] = min(mhsym(TrIx( 1:10),YY).^2 + mhsym(TrIx( 1:10),ZZ).^2); TrR = TrIx(TrR+10); [~, TrL] = min(mhsym(TrIx(11:20),YY).^2 + mhsym(TrIx(11:20),ZZ).^2); TrL = TrIx(TrL); % correct the rotational offsets of each head with respect to the tragus angles = angles(:,:)-angles(TrL,:); %% Asymmetry of the radius radFlip(find(correspondence),:) = radius(correspondence(find(correspondence)),:); %#ok<*FNDSB> radDif = radius - radFlip; % figure; plot(radDif); %% Landmark coordinates of Klingenberg et al, taken from the figure of the mhsym % 7 nasion 0.00 3.39 9.51 * 10^4 % 14 pronasale 0.00 0.00 1.19 % 8 subnasale 0.00 -1.28 1.03 % 13 labiale superior 0.00 -2.87 1.03 % 17 menton 0.00 -7.28 9.33 % 9 endocanthion -4.45 3.09 7.19 4.45 3.09 7.19 % 11 exocanthion -1.64 2.95 7.82 1.76 3.02 7.72 % 15 frontotemporale -4.79 5.85 7.46 4.81 5.66 7.52 NasionSymXYZ = [ 0.00 3.39 9.51] .* 10^4; PronasSymXYZ = [ 0.00 0.10 11.9] .* 10^4; SubnasSymXYZ = [ 0.00 -1.28 10.3] .* 10^4; LabsupSymXYZ = [ 0.00 -2.87 10.3] .* 10^4; MentonSymXYZ = [ 0.00 -7.28 9.33] .* 10^4; EnCanLSymXYZ = [-4.45 3.00 7.19] .* 10^4; EnCanRSymXYZ = [ 4.45 3.00 7.19] .* 10^4; ExCanLSymXYZ = [-1.64 3.00 7.75] .* 10^4; ExCanRSymXYZ = [ 1.76 3.00 7.75] .* 10^4; FrTemLSymXYZ = [-4.80 5.55 7.50] .* 10^4; FrTemRSymXYZ = [ 4.80 5.75 7.50] .* 10^4; % search the closest voxel [~, NasionV] = min((mhsym(:,XX)).^2 + (mhsym(:,YY) - NasionSymXYZ(YY)).^2); [~, PronasV] = min((mhsym(:,XX)).^2 + (mhsym(:,YY) - PronasSymXYZ(YY)).^2); [~, SubnasV] = min((mhsym(:,XX)).^2 + (mhsym(:,YY) - SubnasSymXYZ(YY)).^2); [~, LabsupV] = min((mhsym(:,XX)).^2 + (mhsym(:,YY) - LabsupSymXYZ(YY)).^2); [~, MentonV] = min((mhsym(:,XX)).^2 + (mhsym(:,YY) - MentonSymXYZ(YY)).^2); [~, EnCanLV] = min((mhsym(:,XX) - EnCanLSymXYZ(XX)).^2 + (mhsym(:,YY) - EnCanLSymXYZ(YY)).^2); [~, EnCanRV] = min((mhsym(:,XX) - EnCanRSymXYZ(XX)).^2 + (mhsym(:,YY) - EnCanRSymXYZ(YY)).^2); [~, ExCanLV] = min((mhsym(:,XX) - ExCanLSymXYZ(XX)).^2 + (mhsym(:,YY) - ExCanLSymXYZ(YY)).^2); [~, ExCanRV] = min((mhsym(:,XX) - ExCanRSymXYZ(XX)).^2 + (mhsym(:,YY) - ExCanRSymXYZ(YY)).^2); [~, FrTemLV] = min((mhsym(:,XX) - FrTemLSymXYZ(XX)).^2 + (mhsym(:,YY) - FrTemLSymXYZ(YY)).^2); [~, FrTemRV] = min((mhsym(:,XX) - FrTemRSymXYZ(XX)).^2 + (mhsym(:,YY) - FrTemRSymXYZ(YY)).^2); Nasion = angles(NasionV,:); Pronas = angles(PronasV,:); Subnas = angles(SubnasV,:); Labsup = angles(LabsupV,:); Menton = angles(MentonV,:); EnCant = (angles(EnCanRV,:) + angles(EnCanLV,:))./2; ExCant = (angles(ExCanRV,:) + angles(ExCanLV,:))./2; FrTemp = (angles(FrTemRV,:) + angles(FrTemLV,:))./2; MeanLandmark = mean(reshape([Nasion Pronas Subnas Labsup Menton EnCant ExCant FrTemp],[200,8]),2); Landmark = {'Nasion','Pronas','Subnas','Labsup','Menton','EnCant','ExCant','FrTemp','Mean'}; MeanLM = [ ... round(mean(Nasion)*180/pi,2); round(mean(Pronas)*180/pi,2); ... round(mean(Subnas)*180/pi,2); round(mean(Labsup)*180/pi,2); ... round(mean(Menton)*180/pi,2); round(mean(EnCant)*180/pi,2); ... round(mean(ExCant)*180/pi,2); round(mean(FrTemp)*180/pi,2); ... round(mean(MeanLandmark)*180/pi,2)]; SD_LM = [ ... round(std(Nasion)*180/pi,2); round(std(Pronas)*180/pi,2); ... round(std(Subnas)*180/pi,2); round(std(Labsup)*180/pi,2); ... round(std(Menton)*180/pi,2); round(std(EnCant)*180/pi,2); ... round(std(ExCant)*180/pi,2); round(std(FrTemp)*180/pi,2); ... round(std(MeanLandmark)*180/pi,2)]; Landmarks = table(MeanLM,SD_LM,'RowNames',Landmark'); writetable(Landmarks,'Landmarks.txt','WriteRowNames',true,'Delimiter','tab'); %% mark the rendered faces hd200sym(TrIx,RR) = 0; hd200sym(TrIx,GG) = 1; hd200sym(TrIx,BB) = 1; hd200mean(TrIx,RR)= 0; hd200mean(TrIx,GG) = 1; hd200mean(TrIx,BB) = 1; %% average distorted head exaggerate = 5; Diff = hd200mean-hd200sym; Diff(:,4:6) = 0; Diff(abs(Diff)>5000) = 0; %remove non-matched Diff = hd200sym + exaggerate*Diff; Diff(TrIx,RR)= 0; Diff(TrIx,GG) = 1; Diff(TrIx,BB) = 1; % mark tragi figure; renderhead(Diff); title('the exaggerated asymmetry'); figure; renderhead(rotatehead(hd200sym,'y',90)); title('the average symmetric of all faces'); figure; renderhead(hd200sym); title('the average symmetric of all faces');% %figure; renderhead(hd200mean); title('the average of all faces'); %figure; renderhead(rotatehead(hd200mean,'y',90)); title('the average of all faces'); %% Statistics % The student-t (or Mahalanobis, or d-prime) is best: mangle = nanmean(angles, 2)*180/pi; sangle = nanstd( angles,0,2)*180/pi; mradii = nanmean(radDif, 2); sradii = nanstd( radDif,0,2); % remove extreme vals in non-matched coordinates sangle(abs(mangle)>1.5) = nan; mangle(abs(mangle)>1.5) = nan; mangle( sangle<0.01) = nan; sangle( sangle<0.01) = nan; sangle(abs(mradii)>2500) = nan; mradii(abs(mradii)>2500) = nan; tAng = sqrt(200) * mangle./sangle; tRad = sqrt(200) * mradii./sradii; sangle(abs(tRad)>25) = nan; mangle(abs(tRad)>25) = nan; tRad(abs(tRad)>25) = nan; % figure; plot(mradii) % figure; plot(1./sradii) % figure; plot(tRad) %% t-tests (Bonferroni corrected) % ttest can cope with nans Bonf = 1/length(angles); Alpha = 0.01; % =0.01/75337 [Hleft, ~,~] = ttest(angles',0,'alpha',Alpha*Bonf,'tail','left'); [Hright,~,~] = ttest(angles',0,'alpha',Alpha*Bonf,'tail','right'); Hleft( isnan(Hleft)) = 0; Hright(isnan(Hright))= 0; Faceregion = find(logical(Hright)); FAngle = nanmean(mangle(Faceregion)); EAngle = nanmean(mangle(logical(Hleft))); %% this gives pos and neg t values a different color: % hd = Diff; % hd(:,4) =-min(tAng, 0)/16; % hd(:,5) = max(tAng, 0)/16; % hd(:,6) = hd(:,5); % figure; renderhead(hd); % title('t-statistic of the difference head -angle'); %% the same for the radius: % hd = Diff; % hd(:,4) =-min(tRad, 0)/16; % min(tRad) = -20.6 % hd(:,5) = max(tRad, 0)/16; % max(tRad) = 20.6 % hd(:,6) = hd(:,5); % H=figure; renderhead(hd); % title('t-statistic of the difference head -radius'); %% Angle % hd = Diff; % hd(:,4:6) = repmat(mean(hd(:,4:6),2),[1 3]); % hd(logical(Hleft),RR) = mangle(logical(Hleft)) / min(mangle); % hd(logical(Hleft),GG) = 0; % hd(logical(Hleft),BB) = 0; % hd(logical(Hright),RR) = 0; % hd(logical(Hright),GG) = mangle(logical(Hright)) / max(mangle); % hd(logical(Hright),BB) = mangle(logical(Hright)) / max(mangle); % figure; renderhead(hd); % title(sprintf('angle of the difference head (av %.2f deg in facial region)',FAngle)); %% Radius % hd = Diff; % hd(:,4:6) = repmat(mean(hd(:,4:6),2),[1 3]); % hd(mradii<0,RR) = mradii(mradii<0) / min(mradii); % hd(mradii>0,GG) = mradii(mradii>0) / max(mradii); % hd(mradii>0,BB) = mradii(mradii>0) / max(mradii); % figure; renderhead(hd); % title('radius of the difference head'); %figure; contourc(hd200mean(:,XX),hd200mean(:,YY),mangle) % YY : 0 615 1230 1845 2460 ... 11070/615 figure; hold on; for i=1:3 xpo=mhead(find(symhd200L(:,YY)==(i-2)*615*65), XX); ang=mangle(symhd200L(:,YY)==(i-2)*615*65); %ang=mradii(symhd200L(:,YY)==(i-2)*615*65); % = 39975 plot(xpo,ang); end %40000/615 % Fig-significant : hd = Diff; hd(:,4:6) = hd200mean(:,4:6); hd(logical(Hleft), RR) = 0.3*hd(logical(Hleft), RR); hd(logical(Hright),GG) = 0.3*hd(logical(Hright),GG); figure; renderhead(hd); title(sprintf('significant rotation at alpha=%0.2f (Bonferroni corrected); non corrected=%3.1e)',Alpha,Alpha*Bonf)); %% plot the landmarks and age from Klingenberg et al 2010, along with the adult data Age = [4.400 5.200 13.20 13.70 18]; PronasYr = [1.627 1.523 0.379 0.603 mean(Pronas)*180/pi]; LabialYr = [1.704 1.606 0.637 0.603 mean(Labsup)*180/pi]; NasionYr = [1.978 1.614 0.576 0.864 mean(Nasion)*180/pi]; SubnasYr = [1.847 1.863 0.634 0.753 mean(Subnas)*180/pi]; FrTempYr = [1.536 1.853 0.879 1.290 mean(FrTemp)*180/pi]; MentonYr = [2.013 2.437 1.305 0.925 mean(Menton)*180/pi]; ExCantYr = [2.357 3.084 0.977 0.993 mean(ExCant)*180/pi]; EnCantYr = [2.736 3.288 1.087 1.201 mean(EnCant)*180/pi]; AvLandYr = [1.975 2.158 0.809 0.904 mean(MeanLandmark)*180/pi]; FontSzDef = get(gca,'FontSize'); H=figure; hold on; set(gca,'FontSize',14); plot(Age, AvLandYr,'LineWidth',4); plot(Age, EnCantYr, '->'); plot(Age, ExCantYr, '-<'); plot(Age, MentonYr, '-*'); plot(Age, FrTempYr, '-+'); plot(Age, LabialYr, '-d'); plot(Age, SubnasYr, '-p'); plot(Age, NasionYr, '-x'); plot(Age, PronasYr, '-o'); legend('average','endocanthion','exocanthion','menton','frontotemporale','labiale superior', ... 'subnasale','nasion','pronasale'); xlabel('age (years)');ylabel('asymmetry (deg)') saveas(H,'FigAgeRelation.eps') saveas(H,'FigAgeRelation.pdf') set(gca,'FontSize',FontSzDef); % return to default fontsize FacePictures = tdfread('FacePicsResult.txt'); FacePicsAngle= FacePictures.Angle; H=figure; subplot(1,2,1); histfit(EnCant*180/pi); title('3D analysis'); xlabel('asymmetry (deg)');ylabel('N'); subplot(1,2,2); histfit(FacePicsAngle,15);title('Picture analysis');xlabel('asymmetry (deg)');ylabel('N'); saveas(H,'FigHistogram.pdf');% saveas(H,'FigHistogram.eps') %% regression analysis H=figure; [R1,P1] = corrcoef(ExCant,EnCant);Lab1 = sprintf('correlation: r=%.2f',R1(1,2)); [R2,P2] = corrcoef(Nasion,Menton);Lab2 = sprintf('correlation: r=%.2f',R2(1,2)); [R3,P3] = corrcoef(Subnas,Pronas);Lab3 = sprintf('correlation: r=%.2f',R3(1,2)); [R4,P4] = corrcoef(ExCant,Pronas);Lab4 = sprintf('correlation: r=%.2f',R4(1,2)); % [R4,P4] = corrcoef(ExCant,FrTemp);Lab4 = sprintf('correlation: r=%.2f',R4(1,2)); subplot(2,2,4); plot(ExCant,EnCant,'.'); daspect([1 1 1]);xlabel('Exocanthion');ylabel('Endocanthion');pbaspect([5 4 5]);title(Lab1); subplot(2,2,2); plot(Nasion,Menton,'.'); daspect([1 1 1]);xlabel('Nasion');ylabel('Menton');pbaspect([5 4 5]);title(Lab2); subplot(2,2,1); plot(Subnas,Pronas,'.'); daspect([1 1 1]);xlabel('Subnasion');ylabel('Pronasion');pbaspect([5 4 5]);title(Lab3); subplot(2,2,3); plot(ExCant,Pronas,'.'); daspect([1 1 1]);xlabel('Exocanthion');ylabel('Pronasion');pbaspect([5 4 5]);title(Lab4); % subplot(2,2,4); plot(ExCant,FrTemp,'.'); daspect([1 1 1]);xlabel('Exocanthion');ylabel('Frontotemporale');pbaspect([5 4 5]);title(Lab4); saveas(H,'FigCorrel.pdf'); %saveas(H,'FigCorrel.eps') % Nasion % Pronas % Subnas % Labsup % Menton % EnCant % ExCant % FrTemp