==== Automatic EBSD map stiching ==== ---- So let's say you have two maps A and B and you want to stich them automatically to get the big map AB without any seams: {{:recherche:methodes:ebsd:stich.png?400|}} The following matlab function merges two mtex ebsd objects using [[https://fr.mathworks.com/help/images/ref/imregcorr.html|the phase correlation method]] as explained in [[https://seafile.emse.fr/d/501a9886d50b4466bbdf/files/?p=%2FSimonBreumier_Stiching.pptx|this presentation]]: function [ebsd_s,t_x, t_y] = stich_map(ebsd1, ebsd2, fact, axis, phase) % Debruitage des carto ebsd1_g = fill(ebsd1(phase)); ebsd2_g = fill(ebsd2(phase)); % Conversion de la carto en image (nuance de gris IPFZ) im1 = ebsd1_g.gridify.orientations; im2 = ebsd2_g.gridify.orientations; ipfKey = ipfHSVKey(ebsd1(phase)); [im1y, im1x] = size(im1); [im2y, im2x] = size(im2); im1 = im2uint8(rgb2gray(reshape(ipfKey.orientation2color(im1), [im1y, im1x, 3]))); im2 = im2uint8(rgb2gray(reshape(ipfKey.orientation2color(im2), [im2y, im2x, 3]))); [im1x, im1y] = invert_xy(im1x, im1y); [im2x, im2y] = invert_xy(im2x, im2y); % Calcul de la conversion pixel => coordonnees carto p2m_x1 = (max(ebsd1.x) - min(ebsd1.x))/im1y; p2m_y1 = (max(ebsd1.y) - min(ebsd1.y))/im1x; p2m_x2 = (max(ebsd2.x) - min(ebsd2.x))/im2y; p2m_y2 = (max(ebsd2.y) - min(ebsd2.y))/im2x; % Choix des zones de recherche if axis == 'y' taille_recouv = round(im2x/fact); x1d=round(im1x-2.5*taille_recouv); x1f=im1x; x2d=1; x2f=taille_recouv; im1p = im1(x1d:x1f,:); im2p = im2(x2d:x2f,:); shift11_x = 0; shift11_y = im1x-2.5*taille_recouv; elseif axis == 'x' taille_recouv = round(im2y/fact); x1d=round(im1y-2.5*taille_recouv); x1f=im1y; x2d=1; x2f=taille_recouv; im1p = im1(:,x1d:x1f); im2p = im2(:,x2d:x2f); shift11_x = im1y-2.5*taille_recouv; shift11_y = 0; else fprintf("Wrong axis specification ('x' or 'y' only)"); end % Calcul de la transformation a appliquer tformEstimate21 = imregcorr(im2p,im1p,'transformtype','translation'); t_x = (tformEstimate21.T(3,1) + shift11_x)*p2m_x1; t_y = (tformEstimate21.T(3,2)+ shift11_y)*p2m_y1; %Fusion des objets EBSD ebsd_t = ebsd2; ebsd_t.x = ebsd_t.x + t_x; ebsd_t.y = ebsd_t.y + t_y; ebsd_t = ebsd_t(ebsd_t.y > max(ebsd1.y) | ebsd_t.x > max(ebsd1.x)); ebsd_s = [ebsd1 ebsd_t]; end The function takes the following entry parameters: * **ebsd1, ebsd2**: the two ebsd object to stich (2 must be at the right of 1 for the x axis, or at the bottom of 1 for the y axis) * **fact**: the length ratio of the zone to use in image 2 for the phase correlation method * **axis**: chose either to stich map 1 on top of map 2 or map1 at the left of map 2 (could be generalized I guess..) * **phase**: the phase to use to generate the orientation field used for stiching Here's an example using the function to automatically stich several maps along the y and x axis successively. The scripts pauses after each column and after stiching two columns to ask if the stiching is correct. Is not, the **fact** parameter can be changed until it works. clear close all clc %% Choix des paramètres % Nom des images carto1 = 'carto/1/Echantillon L14 1800C 100s grande zone_0_0.crc'; preset = 'carto/1/Echantillon L14 1800C 100s grande zone_'; ix_max = 7; iy_max = 7; verif=0; phase = 'Tungstène'; %% CS = {... 'notIndexed',... crystalSymmetry('m-3m', [3.2 3.2 3.2], 'mineral', phase, 'color', [0.53 0.81 0.98])}; % plotting convention setMTEXpref('xAxisDirection','east'); setMTEXpref('zAxisDirection','intoPlane'); f = waitbar(0,'1','Name','Stiching maps...',... 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)'); setappdata(f,'canceling',0); factx = 7.; facty = 7.; diffx = 0; diffy = 0; for j=0:iy_max done = 0; while ~done ebsd1 = EBSD.load(carto1,CS,'interface','crc','convertEuler2SpatialReferenceFrame'); ebsd_s = ebsd1; for i=0:ix_max waitbar((j*ix_max+i)/((ix_max+1)*(iy_max+1)),f,"Stitching maps..."); if getappdata(f,'canceling') break end if ~(i==0) carto_actu = [preset num2str(j) '_' num2str(i) '.crc']; ebsd2 = EBSD.load(carto_actu,CS,'interface','crc','convertEuler2SpatialReferenceFrame'); [ebsd_s, tx, ty] = stich_map(ebsd_s, ebsd2, factx, 'y', phase); end end plot(ebsd_s(phase), ebsd_s(phase).orientations) done = input("Is the column correct ? (0=no - 1=yes)"); if ~done factx = input("Specify a new covering zone parameter: "); end end if j==0 ebsd_final = ebsd_s; else done = 0; while ~done carto_j_0 = [preset num2str(j-1) '_0.crc']; carto_jp1_0 = [preset num2str(j) '_0.crc']; ebsd_j_0 = EBSD.load(carto_j_0,CS,'interface','crc','convertEuler2SpatialReferenceFrame'); ebsd_jp1_0 = EBSD.load(carto_jp1_0,CS,'interface','crc','convertEuler2SpatialReferenceFrame'); [tmp, tx, ty] = stich_map(ebsd_j_0, ebsd_jp1_0, facty, 'x', phase); ebsd_s_tmp = ebsd_s; ebsd_s_tmp.x = ebsd_s_tmp.x + tx + diffx; ebsd_s_tmp.y = ebsd_s_tmp.y + ty + diffy; ebsd_f_t = [ebsd_final ebsd_s_tmp]; plot(ebsd_f_t(phase), ebsd_f_t(phase).orientations) done = input("Is the new global map correct ? (0=no - 1=yes)"); if ~done facty = input("Specify a new covering zone parameter: "); end end diffx = diffx + tx; diffy = diffy + ty; ebsd_final = ebsd_f_t; % ebsd_final.x = ebsd_final.x - min(ebsd_final.x); % ebsd_final.y = ebsd_final.y - min(ebsd_final.y); end carto1 = [preset num2str(j+1) '_0.crc']; end delete(f) ebsd_final.export("Carto_global.ctf")