TrueEBSD distortion correction on a WC-Co composite edit page

Authors: Vivian Tong; Stefan Olovsjö, Seco Tools AB, R&D Materials and Technology, 737 82 Fagersta, Sweden. Contact: vivian.tong@extern.tu-freiberg.de

An EBSD map and an SEM image of the same area never quite line up. The beam drifts during the scan, the camera moves between acquisitions, the specimen is tilted. Overlay the two and features are out by tens of pixels, which is enough to ruin any measurement that needs both.

TrueEBSD corrects that, so every pixel of the map and every pixel of the images refer to the same point on the sample. The method is described in Tong et al., arXiv 2605.00703.

This script needs the TrueEBSD toolbox, which is not part of MTEX and is distributed separately under Apache-2.0: https://github.com/vtvivian/mtex-trueebsd 

addpath(genpath('<path to mtex-trueEbsd>'))

Use the mtex7-compat branch with MTEX 7. It also needs MATLAB R2024a or newer and the Image Processing, Curve Fitting, and Statistics and Machine Learning toolboxes. Expect minutes, not seconds.

The data

An EBSD map of a WC-Co composite taken at 20 kV, and four SEM images of the same area, all in ebsd.opt.trueEbsdImgs:

  1. fsdB3 — colour image from the three FSD detectors below the EBSD camera, camera retracted 20 mm from the mapping position;
  1. fsdT3 — the same beam scan as fsdB3, from the FSD detectors above the camera;
  1. fsdT1 — top FSD detectors, camera at the mapping position;
  1. fsdT10 — as fsdT1, but with the beam at 10 kV.

pixSzImg is their pixel size in microns, the same for all four. The map itself is represented by its band contrast, ebsd.bc.

The convention set here only decides which way up the figures come out. It has no effect on the correction.

plottingConvention.default('y↓→x')

mtexdata trueEbsdWCCo

display(ebsd.opt.trueEbsdImgs)
ans = plottingConvention (y↓→x)
 
ebsd = EBSDsquare (y↓→x)
 
 Phase  Orientations     Mineral         Color  Symmetry  Crystal reference frame
     0  6657 (0.85%)  notIndexed          none                                   
     1  714592 (91%)         W C  LightSkyBlue      -6m2        X||a*, Y||b, Z||c
     2    62791 (8%)      Co-fcc  DarkSeaGreen      m-3m                         
     3   2392 (0.3%)      Co-hcp     Goldenrod     6/mmm        X||a*, Y||b, Z||c
 
 Properties: bc, bs, bands, MAD, quality, oldId
 Scan unit : um
 X x Y x Z : [0 → 41] x [0 → 30] x [0 → 0]
 Normal vector: (0,0,1)
 Square grid  :768 x 1024
 
  struct with fields:

       fsdB3: [1536×2048×3 double]
       fsdT3: [1536×2048 double]
       fsdT1: [1536×2048 double]
      fsdT10: [1536×2048 double]
    pixSzImg: 0.0199

Build the sequence

TrueEBSD does not jump straight from the EBSD map to the reference image. It steps through the images one pair at a time, correcting one kind of distortion at each step. That is the point of having four images rather than one: each pair differs by something simple enough to model, where the map and the final image differ by everything at once.

So the order matters. It runs from the most distorted to the ground truth, and each entry names the distortion between itself and the next one:

  • 'shift-drift' — beam drift during the scan, plus a rigid offset
  • 'shift' — a rigid offset, from the camera moving
  • 'tilt' — the specimen tilt, seen differently at a different kV
  • 'true' — nothing to correct; the pair came from one beam scan

The reference image is the one you trust: here fsdT10, taken at the mapping position. Everything is brought onto it.

'name' is what each image will be called once it is attached to the map at the end, so you read the result as ebsd.fsdT1.

'highContrast' says the image has usable edges to match on. 'edgePadWidth' widens those edges where boundaries are blurred.

% a light box filter first — cross-correlation dislikes noise
img = ebsd.opt.trueEbsdImgs;
img.fsdB3  = rescale(imboxfilt(img.fsdB3,3));
img.fsdT3  = rescale(imboxfilt(img.fsdT3,3));
img.fsdT1  = rescale(imboxfilt(img.fsdT1,3));
img.fsdT10 = rescale(imboxfilt(img.fsdT10,3));

imgList = createArray(5,1,'distortedImg');
imgList(1) = distortedImg('bc','shift-drift', ebsd, 'name','bcImg', ...
  'highContrast',1, 'edgePadWidth',3);
imgList(2) = distortedImg(img.fsdB3, 'true',  'dxy', img.pixSzImg, 'name','fsdB3',  'highContrast',1, 'edgePadWidth',5);
imgList(3) = distortedImg(img.fsdT3, 'shift', 'dxy', img.pixSzImg, 'name','fsdT3',  'highContrast',1, 'edgePadWidth',5);
imgList(4) = distortedImg(img.fsdT1, 'tilt',  'dxy', img.pixSzImg, 'name','fsdT1',  'highContrast',1, 'edgePadWidth',5);
imgList(5) = distortedImg(img.fsdT10,'true',  'dxy', img.pixSzImg, 'name','fsdT10', 'highContrast',1, 'edgePadWidth',3)
imgList = distortedImg (y↓→x)
 
              image    pixel     frame   distortion  contrast                 EBSD
 1       768 x 1024  0.04 µm      y↓→x  shift-drift      high  W C, Co-fcc, Co-hcp
 2  1536 x 2048 x 3  0.02 µm  row↓→col         true      high                    -
 3      1536 x 2048  0.02 µm  row↓→col        shift      high                    -
 4      1536 x 2048  0.02 µm  row↓→col         tilt      high                    -
 5      1536 x 2048  0.02 µm  row↓→col         true      high                    -

Plot the sequence to check the images really do cover the same area. Note how different the contrasts look — that is why matching is done on edges rather than on grey values.

plot(imgList,'TrueEBSD starting image sequence')

Say how the map and the images are related

Comparing two pictures pixel by pixel only means something if you know which way round each one is stored. TrueEBSD asks rather than guesses, because guessing wrong does not fail loudly — it returns a finished-looking result that is quietly nonsense.

Usually there is nothing to do: a map and images collected in one session are stored the same way and TrueEBSD confirms it. This file is the exception. Its EBSD map was written with the columns in the opposite order from the images, so pixel (1,1) of the map and pixel (1,1) of an image are at opposite ends of the sample.

You will not see this in the plot above. Every plot draws each entry the way that entry says it should be drawn, so on screen they already agree — the difference is in how the numbers are stored, and that is what matters when they are compared. transformReferenceFrame puts them all in one frame, working the relation out from what the map and the images each already say about themselves.

Where that is not enough — an image from a separate session, at an unknown rotation — you state the relation instead. See the frames page.

imgList = imgList.transformReferenceFrame(imgList(1).frame)
imgList = distortedImg (y↓→x)
 
              image    pixel  frame   distortion  contrast                 EBSD
 1       768 x 1024  0.04 µm   y↓→x  shift-drift      high  W C, Co-fcc, Co-hcp
 2  1536 x 2048 x 3  0.02 µm   y↓→x         true      high                    -
 3      1536 x 2048  0.02 µm   y↓→x        shift      high                    -
 4      1536 x 2048  0.02 µm   y↓→x         tilt      high                    -
 5      1536 x 2048  0.02 µm   y↓→x         true      high                    -

Now the job can be built. It holds the whole workflow, and each step below adds to it. It also re-checks that the sequence is in one frame, and stops with trueEbsd:frameMismatch if not.

job = trueEbsd(imgList)
job = trueEbsd (as imported)
 
              image   distortion  shift  residual
 1       768 x 1024  shift-drift      -         -
 2  1536 x 2048 x 3         true      -         -
 3      1536 x 2048        shift      -         -
 4      1536 x 2048         tilt      -         -
 5      1536 x 2048         true      -         -

Put everything on one pixel grid

The map and the images cover the same area at different pixel sizes. pixelSizeMatch resamples them all onto the finest one, so pixel (i,j) means roughly the same place in each.

Images are interpolated linearly. EBSD data is not: orientations and phase labels have no meaningful average, so the nearest measured point is used.

Nothing has been corrected yet — this is only bookkeeping.

job.pixelSizeMatch
using default pixel size of 0.019875 um, minimum from imgList
 
ans = trueEbsd (pixel size matched)
 
              image   distortion  shift  residual
 1      1536 x 2048  shift-drift      -         -
 2  1536 x 2048 x 3         true      -         -
 3      1536 x 2048        shift      -         -
 4      1536 x 2048         tilt      -         -
 5      1536 x 2048         true      -         -
 
 common grid: 1536 x 2048 at 0.02 µm

The resampled sequence is in job.resizedList.

job.resizedList
ans = distortedImg (y↓→x)
 
              image    pixel  frame   distortion  contrast                 EBSD
 1      1536 x 2048  0.02 µm   y↓→x  shift-drift      high  W C, Co-fcc, Co-hcp
 2  1536 x 2048 x 3  0.02 µm   y↓→x         true      high                    -
 3      1536 x 2048  0.02 µm   y↓→x        shift      high                    -
 4      1536 x 2048  0.02 µm   y↓→x         tilt      high                    -
 5      1536 x 2048  0.02 µm   y↓→x         true      high                    -

[Optional] Adjust the matching windows

Distortions are measured by cutting both images of a pair into small boxes and cross-correlating each box with its partner. That gives a local shift at each box, and those shifts are fitted to the distortion model.

Two settings control it, both in pixels, so they have to be written after pixelSizeMatch and before calcDistortion:

  • ROISize — the box width, a power of two. Rule of thumb: at least four times the shift you expect to measure. Too small and the boxes have nothing in common to lock onto; too large and real distortion varies within a box.
  • NumROI — how many boxes across and down. Roughly one per grain works.

The defaults suit a typical polycrystal map, so most people never touch this. Here the first box size is set deliberately far too small, to show what happens when it is wrong.

customSetXCF1.ROISize = 2^round(log2(32));   % deliberately too small
customSetXCF1.NumROI = struct;
customSetXCF1.NumROI.x = 40;                 % about one box per grain across
customSetXCF1.NumROI.y = round(customSetXCF1.NumROI.x * ...
  size(job.resizedList(1).img,1)/size(job.resizedList(1).img,2)); % keep the aspect ratio
customSetXCF1.xcfImg = 'edge';               % match on edges, or 'img' for grey values

customSetXCF2 = customSetXCF1;
customSetXCF2.ROISize = 2^round(log2(128));

% a whole settings struct at once
job.resizedList(1).setXCF(2) = customSetXCF1;
job.resizedList(3).setXCF(1) = customSetXCF2;

% or one setting at a time
job.resizedList(1).setXCF(1).ROISize = 2^round(log2(64));
job.resizedList(3).setXCF(1).xcfImg = 'img';
job.resizedList(4).setXCF(1).xcfImg = 'img';
job.resizedList(5).setXCF(1).xcfImg = 'img';

Measure the distortion

These are the pictures that will actually be matched: edge transforms where xcfImg is 'edge', grey values where it is 'img'. A band contrast map and a backscatter image have nothing in common as grey values, but their grain boundaries fall in the same places — which is why edges are the default.

plot(job.resizedList,'TrueEBSD image sequence for cross-correlation','xcf')

'fitErr' re-measures the shifts after each correction and reports what is left over. That residual is how you tell whether it worked: around a pixel or less is good.

Because the first box size was set too small, the EBSD map step leaves a residual above two pixels. TrueEBSD notices, doubles the box size and tries again, repeating until the residual comes down or the box outgrows the image. Watch the printed table.

Steps marked 'true' are skipped — nothing separates that pair, so their shift is taken as zero whatever the residual says.

job.calcDistortion('fitErr')
fitting the distortion across 5 maps, 4 hops

  hop   distortion    model               ROI      shift   residual
  2→1   shift-drift   poly11            64 px   11.53 px
                      linearinterp      32 px    5.60 px
                      after fitting                         5.00 px
                      residual over 2 px, retrying with a 128 px ROI
  2→1   shift-drift   poly11           128 px   14.09 px
                      linearinterp      64 px    6.00 px
                      after fitting                         3.96 px
                      residual over 2 px, retrying with a 256 px ROI
  2→1   shift-drift   poly11           256 px   13.64 px
                      linearinterp     128 px    5.50 px
                      after fitting                         1.92 px
  3→2   true          not correlated        -    0.00 px
                      not fitted                            2.05 px
  4→3   shift         poly11           128 px   12.63 px
                      after fitting                         0.95 px
  5→4   tilt          projective       512 px   14.49 px
                      poly11           512 px    0.41 px
                      poly22           512 px    0.33 px
                      after fitting                         0.29 px

 
ans = trueEbsd (shifts calculated)
 
              image   distortion    shift  residual
 1      1536 x 2048  shift-drift   6.2 px    1.9 px
 2  1536 x 2048 x 3         true     0 px      2 px
 3      1536 x 2048        shift    13 px   0.95 px
 4      1536 x 2048         tilt  0.42 px   0.29 px
 5      1536 x 2048         true        -         -
 
 common grid: 1536 x 2048 at 0.02 µm

Correct it

Each map is now moved by the shifts of every step between it and the reference. The reference itself does not move.

Resampling is nearest-neighbour throughout, so no orientation and no phase label is ever invented by averaging two real measurements.

job.undistort

plot(job.undistortedList,'TrueEBSD image sequence after alignment')
ans = trueEbsd (undistorted)
 
              image   distortion    shift  residual
 1      1536 x 2048  shift-drift   6.2 px    1.9 px
 2  1536 x 2048 x 3         true     0 px      2 px
 3      1536 x 2048        shift    13 px   0.95 px
 4      1536 x 2048         tilt  0.42 px   0.29 px
 5      1536 x 2048         true        -         -
 
 common grid: 1536 x 2048 at 0.02 µm

Use the result

Every image is now attached to the EBSD map as a per-pixel property, under the 'name' given earlier. So ebsd.fsdT1 is just another map property, and plot(ebsd,ebsd.fsdT1) works like any other plot — no conversion, and it stays with the map through cropping, gridding and indexing.

Plotting them back onto the map is also the quickest check that nothing came out the wrong way round.

fsdB3 is a colour image and keeps all three channels. Plotting onto a map needs one value per pixel, so it is averaged to grey here.

ebsdOut = job.undistortedList(1).ebsd;

figure
nextAxis
plot(ebsdOut('W C'), ebsdOut('W C').orientations, 'coordinates','on')
title('Undistorted MTEX EBSD map (WC IPF out of screen)','Color','k')

for n = 1:numel(job.undistortedList)

  im = ebsdOut.(job.undistortedList(n).name);
  if size(im,3) > 1, im = mean(im,3); end

  nextAxis
  plot(ebsdOut, im, 'coordinates','on')
  mtexColorMap gray
  title(['Undistorted ' job.undistortedList(n).name],'Color','k')
end

Finish

The map and the images now overlay pixel for pixel, and ebsdOut is an ordinary MTEX EBSD map that happens to carry four SEM images as properties. Anything you would normally do with a map works from here — including turning a thresholded image into a phase, which is what the copper page does with its voids.

For this sample the next step is measuring WC grain contiguity: which boundaries are WC against WC, and which are WC against binder. That needs the orientation data and the image contrast to agree pixel for pixel, which is what the last few minutes bought.