function [faces surftri surfNormali T2F F2T]=surftriMod(p,t)
% Form all faces, non-duplicates are surface triangles
Nt=size(t,1);
faces=[t(:,[1,2,3]);
       t(:,[1,2,4]);
       t(:,[1,3,4]);
       t(:,[2,3,4])];
node4=[t(:,4);t(:,3);t(:,2);t(:,1)];
faces=sort(faces,2);
[foo,ix,jx]=unique(faces,'rows');
vec=histc(jx,1:max(jx));
qx=find(vec==1);
surftri=faces(ix(qx),:);
node4=node4(ix(qx));
faces=foo;
T2F=[jx(1:Nt) jx((Nt+1):(2*Nt)) jx((2*Nt+1):(3*Nt)) jx((3*Nt+1):(4*Nt))];
F2T=-1*ones(size(faces,1),2);

% NON RIESCO A VETTORIZZARE QUESTA PARTE DI CODICE!!!! ODIO MATLAB!!!!!
% SONO INCAZZATISSIMO.... ADESSO ME NE VADO SENNO' SPACCO TUTTO GGRRRRR

for i=1:Nt
    App=union(F2T(T2F(i,1),:),[i]);
    F2T(T2F(i,1),:)=App((end-1):end);
    App=union(F2T(T2F(i,2),:),[i]);
    F2T(T2F(i,2),:)=App((end-1):end);
    App=union(F2T(T2F(i,3),:),[i]);
    F2T(T2F(i,3),:)=App((end-1):end);
    App=union(F2T(T2F(i,4),:),[i]);
    F2T(T2F(i,4),:)=App((end-1):end);
end

% App=[(1:Nt)'; (1:Nt)'; (1:Nt)'; (1:Nt)']*ones(1,size(faces,1));
% [Row Col Ind]=find(jx*ones(1,size(faces,1))-((1:size(faces,1))'*ones(1,length(jx)))'==0);
% F2T(1,:)=union(App(Row(find(Col==1)),1));

% Orientation and normals
v1=p(surftri(:,2),:)-p(surftri(:,1),:);
v2=p(surftri(:,3),:)-p(surftri(:,1),:);
v3=p(node4,:)-p(surftri(:,1),:);
surfNormali=cross(v1,v2,2);
ix=find(dot(surfNormali,v3,2)>0);
surfNormali=surfNormali./(sqrt(sum(surfNormali.^2,2))*[1 1 1]);
surftri(ix,[2,3])=surftri(ix,[3,2]);
surfNormali(ix,:)=-surfNormali(ix,:);