Skip to content
Snippets Groups Projects
Commit 3ef96488 authored by Laura Nervo's avatar Laura Nervo Committed by Nicola Vigano
Browse files

Color maps : fixing bugs for c-axis and IPF maps.


Signed-off-by: default avatarLaura Nervo <laura.nervo@esrf.fr>

git-svn-id: https://svn.code.sf.net/p/dct/code/trunk@629 4c865b51-4357-4376-afb4-474e03ccb993
parent 43d9cc96
No related branches found
No related tags found
No related merge requests found
function c_axismap = gtCaxisCmap(phaseid)
function c_axismap = gtCaxisCmap(phaseid, bodge)
% GTCAXISCMAP Make a colourmap according to c-axis orientation
% c_axismap = gtCaxisCmap(phaseid)
% c_axismap = gtCaxisCmap(phaseid, bodge)
% --------------------------------
% Loads parameters.mat and 4_grains/phase_##/r_vectors.mat
% Saves c_axis.map and c_axismap.mat
% Saves c_axis.mat, map.mat and c_axismap.mat
%
% INPUT:
% phaseid = phase number <int> {1}
%
% OUTPUT:
% c_axismap = c_axis colourmpa calculated from r_vectors list as produced by gtReadGrains
% c_axismap = c_axis colourmap calculated from r_vectors list as produced by gtReadGrains
%
% r_vectors -> c-axis -> phi/psi -> hsl -> rgb colormap
% use psi to give brightness; phi to give colour
......@@ -20,6 +20,9 @@ function c_axismap = gtCaxisCmap(phaseid)
if ~exist('phaseid','var') || isempty(phaseid)
phaseid = 1;
end
if ~exist('bodge','var') || isempty(bodge)
bodge = 0;
end
parameters = [];
r_vectors = [];
......@@ -54,6 +57,10 @@ all_hkls = all_hkils;
for i=1:size(r_vectors,1)
g = Rod2g(r_vectors(i,2:4));
disp([num2str(bodge) ' degrees bodge here!'])
g=[cosd(bodge) sind(bodge) 0; -sind(bodge) cosd(bodge) 0; 0 0 1]*g;
all_normals = (g * normalised_hkls')';
c_axis(i,1)=r_vectors(i,1);
if all_normals(2,3)>=0
......@@ -108,4 +115,7 @@ end
%save the caxis colourmap
save('c_axismap.mat','c_axismap');
map = gtIVPhexCmap(r_vectors,bodge);
sa('map.mat','map');
end % end of function
\ No newline at end of file
function cmap = gtIVPhexCmap(r_vectors)
% map = gtIVPhexCmap(r_vectors)
function cmap = gtIVPhexCmap(r_vectors, bodge)
% map = gtIVPhexCmap(r_vectors, bodge)
if ~exist('bodge','var') || isempty(bodge)
bodge = 0; % degrees
end
cmap = [];
RD=[0 0 1];
......@@ -7,12 +12,12 @@ RD=[0 0 1];
% six fold rotation around z
for i=1:6
a=(i-1)*60;
sym(i).g3 = [...
symm(i).g3 = [...
cosd(a) sind(a) 0;...
-sind(a) cosd(a) 0;...
0 0 1];
% include the mirror
sym(i+6).g3=[1 0 0; 0 -1 0; 0 0 1]*sym(i).g3;
symm(i+6).g3=[1 0 0; 0 -1 0; 0 0 1]*symm(i).g3;
end
map = zeros(length(r_vectors+1),3);
......@@ -20,18 +25,22 @@ map = zeros(length(r_vectors+1),3);
for i=1:length(r_vectors)
g=Rod2g(r_vectors(i, 2:4));
disp([num2str(bodge) ' degrees bodge here!'])
g=[cosd(bodge) sind(bodge) 0; -sind(bodge) cosd(bodge) 0; 0 0 1]*g;
% RD in cryst axes
RDc = g*RD';
% all equivilents by sym
RDc_sym=[];
for jj=1:12
RDc_sym(jj, :)=(sym(jj).g3*RDc)';
RDc_sym(jj, :)=(symm(jj).g3*RDc)';
end
RDc_sym=[RDc_sym; -RDc_sym];
RDc_psi=acosd(RDc_sym(:,3));
RDc_phi=rad2deg(atan2(RDc_sym(:,2), RDc_sym(:,1)));
% pick equiv that is in out triangle
ndx=find(RDc_psi<90 & RDc_phi>=60 & RDc_phi<90);
......
function gtMakePoleFigure2(poles, weights)
% as for gtMakePoleFigure... however, given poles, gtMakePoleFigure uses a
% search radius to find all poles close to each data point on the pole
% figure, thus doing some inherent smoothing. This is also slow for large
% numbers of poles. If the data going in does not require smoothing
% (remains to be seen!) maybe can construct the pole figue data in a faster
% way.
%at the mo this is quick (~10x cf gtMakePoleFigure) but not clever enough
%because the sampling areas get too small near the pole, and thus noise
% here overpowers the signal. Needs something cleverer... Triangle mesh
% from isosurface on a sphere?
%step size in pole figure
inc=0.02;
%change this if you want a 360deg or 90deg view
phimax=pi/2;
%phimax=2*pi;
psimax=pi/2;
phi_range=0:inc:(phimax-inc);
psi_range=0:inc:(psimax-inc);
%before building figure, take only those poles that fall in the range 0-phimax,
%0-psimax (allowing for seach radius)
if phimax==pi/2
%for a 90 degree section of the pole figure
poles=poles(find(poles(:,1)>0 & poles(:,2)>0 & poles(:,3)>0),:);
elseif phimax==2*pi
%for 360 pole figure - reduce only using psi
poles=poles(find(poles(:,3)>cos(psimax)-sin(searchR)),:);
else
disp('dont know how take a subset of the poles input')
end
%preallocate density
density=zeros(length(phi_range), length(psi_range));
%calc phi and psi for all poles
psi=acos(poles(:,3));
phi=atan2(poles(:,2), poles(:,1));
for i=1:length(phi_range)
i
for j=1:length(psi_range)
dum=find(phi>=phi_range(i) & phi<phi_range(i)+inc & psi>=psi_range(j) & psi<psi_range(j)+inc);
density(i, j)=sum(weights(dum));
end
end
%normalise density values
%needs to be wrt random population, and account for the search area
%increment size as function of psi:
inc_area=inc*[-cos(psi_range+inc)+cos(psi_range)];
inc_area=repmat(inc_area, length(phi_range), 1);
density2=density./inc_area;
%total density/total area
total_density=sum(density(:));
total_area=sum(inc_area(:));
density2=10*density2/max(density2(:));
keyboard
phi=0:inc:phimax; %around z-axis
psi=0:inc:psimax; %away from z axis
%radius distorted for stereographic proj
r=tan((psi_range+(inc/2))/2);
[phi,r]=meshgrid((phi_range+(inc/2)),r);
[x,y]=pol2cart(phi,r);
figure
surf(x,y,density2')
hold on
shading('interp')
axis equal
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment