Image Analysis: Matlab Script

The following is an expert from Dewar - Characterization and Evaluation of Aged 20Cr32Ni1Nb Stainless Steels, and references a stainless steel used in this thesis. The content of this article can also be applied to a general case, and should not be limited to to the presented example.

1. Thresholding scripts

1.1. Setting up the Photoshop connections to Matlab

In every distribution for Photoshop there is a plugin for connecting it to Matlab. This is useful if a file has multiple layers that need to be sequentially called to Matlab. Instruction for how to set up this connections can be found in the Adobe/Photoshop/MATLAB directory typically found in Program Files.

1.2. Image analysis toolbox for Matlab

Image processing with Matlab, and the following scripts require the Image Processing Toolbox that can be purchased from Matlab from their website,

1.3. Thresholding Scripts

1.3.1. Photoshop Script

The following section will go through the photoshop script used to threshold EPMA EPMA element maps. The script EPMAColorThreshold can be downloaded at in the Photoshop Scripts folder. This script is specifically developed for the 2032Nb alloy, and uses the nickel, niobium, silicon, and chromium maps. The script is written where the first function SelectAll Phases is the main function, and all subsequent functions are called from it.

Lines 32-55 will initiate a cropping function to crop out the false color map, and get rid of the color scale and the measurement scale. It does this with the nickel map, as nickel is found at a relatively high concentration throughout the entire microstructure. The function CropArea selects the black area outside of the map in the nickel layer, and then inverts the selection, so that the map area is highlighted. This function returns the width and height of this bounding box. The If statement on line 35 is needed if the map is a merged map that has multiple images stacked on the y-axis.

var refLay = docRef.layers["Ni"];
crop = cropArea (refLay);
if (crop.y > 512){

Lines 57-70 selects the niobium map and chooses and red colors that show up on the map. If there are no red pixels on the map then any yellow pixels are selected. A NbC layer is made, and all of the selected pixels are copied over to this layer and changed to white.

var actLay = docRef.layers["Nb"];
var layerRef = docRef.layers[1];

docRef.activeLayer = actLay;

multiColorChooser (10, c.colorRed.lab.l , c.colorRed.lab.a , c.colorRed.lab.b)
//IF NbC does not have red Choose a yellow color to define layer
var isRed = true;
	multiColorChooser (10, c.colorYellow.lab.l , c.colorYellow.lab.a , c.colorYellow.lab.b);
	isRed = false;

var nbcLayer = docRef.layers[0]; = "NbC";

docRef.activeLayer = nbcLayer;


Any other color pixels that have been identified to be associated with the NbC phase are then selected and copied over to the NbC layer, after this is complete all the other pixels in the layer are filled black, and the layer is deselected.

colSpecPick (c.colorDarkRed, c.colorDarkOrg, nbcLayer);
colSpecPick (c.colorOrg, undefined, nbcLayer);
//Add Background layer to NbC layer
var backLayer = docRef.layers[0]; = "Background";
backLayer.move (nbcLayer, ElementPlacement.PLACEAFTER);
docRef.activeLayer = backLayer;
docRef.activeLayer = actLay;

docRef.activeLayer = nbcLayer;

For all of the other phases the lines 99-130 are used to copy relative pixel color intensities over the new phase layers. Lines 99-130 need to formatted for each phase by changing the layer selected in line 104, defining the new layer to store the binary selection (lines 108-109), and the colors selected (lines 113-117). The function colSpecPick finds the colors specified and copies them to the layer specified in the last argument. All of the colors in the EPMA spectra are defined in the spectraColors function and can be found there for naming declarations used for the first two arguments of colSpecPick. The rest of the script makes new layers for a silicon rich phase, and a chromium rich phase.

 //Deselect NbC contents hide NbC layer
//make active layer the Nb layer
docRef.activeLayer = actLay;
// Make GNbRich Layer
var GNbLayer = docRef.layers[0] = "NbRich" 

 //-------------Nb for G-phase --------------------------
colSpecPick (c.colorMidYellow, c.colorYellow, GNbLayer);
colSpecPick (c.colorBritGreen, c.colorMidGreen,GNbLayer);
colSpecPick (c.colorDarkGreen, c.colorTeal, GNbLayer);
colSpecPick (c.colorSkyBlue, c.colorPurple, GNbLayer);
colSpecPick (c.colorMidPurple, undefined, GNbLayer);
//Add Background layer to carbide layer
var backLayer = docRef.layers[0]; = "Background";
backLayer.move (GNbLayer, ElementPlacement.PLACEAFTER);
docRef.activeLayer = backLayer;
docRef.activeLayer = actLay;

docRef.activeLayer = GNbLayer;
1.3.2. Matlab Script

Color thresholding operations can also be performed in Matlab, where this might be more efficient than calling Photoshop to run its thresholding script, and then sending the data (aka. layers) into Matlab for further processing. The function discussed in this section is called EPMAsegmentation, and can be found at under the Matlab/Threshold Scripts directory. Color choices are more limited compared to the Photoshop script where each color designation covers a few color intensities in the EPMA maps.

Calling the thresholding function requires defining 5 arguments; the first argument is the image, the second argument is the color designation (eg. ’r’ for red, ’b’ for blue), and the third argument is a boolean for displaying the image after the function. The forth argument is a boolean for a black background otherwise a grayscale background of the other pixels, and the fifth argument is for changing the selected pixels to white, or otherwise not changing their color. For thresholding multiple colors the following format must used:

nb = EPMAsegmentation(p1, 'w',0,1,1)+...
EPMAsegmentation(p1, 'r',0,1,1)+...
EPMAsegmentation(p1, 'o',0,1,1)+...
EPMAsegmentation(p1, 'y',0,1,1)+...
EPMAsegmentation(p1, 'g',0,1,1)+...
EPMAsegmentation(p1, 't',0,1,1)+...
EPMAsegmentation(p1, 'b',0,1,1);

1.4. Matlab phase fraction script for element mapping

Batch Matlab scripts were developed to calculate area fractions for phases characterized in the microstructure of sequentially aged 2032Nb alloys. The element maps processed by EPMA were compiled into .psd files, and cropped using scripts in the Maps2PSD, and Photoshop Scripts folders at The Matlab script discussed in this section is only applicable to the fully aged component of the in-service 2032Nb alloy discussed in his document, but can be modified to work for other systems, and aging times following the same steps.
 The Matlab script can be downloaded in the Matlab/EPMA Script directory at The Matlab script contains one function for batch processing .psd element maps provided two arguments: the top directory level of the .psd files, and how many subfolders deep the .psd files are located. Lines 5-20 initiate a looping process to go through each .psd file in the directory. Valid .psd files are found by breaking down the path string until the characters ‘cc’ or ‘sc’ are found designating a folder with .psd files from the same samples. The folder title will be used to label the data and should be descriptive in providing what sample it is and how long it was aged for.

function [data] = EPMAInter_GPhase(toplevel, depth)
%UNTITLED3 Summary of this function goes here
%   Detailed explanation goes here

[files,total] = file_list(toplevel,'*.psd', depth);
for i=1:total
 %Load .psd in photoshop
 [pathstr, name, ext] = fileparts(files{i});
 remain = pathstr;
 while true
	 [str, remain] = strtok(remain, '\\');
	 if isempty(str),  break;  end
	 if (strncmpi(str, 'cc', 2)== 1)
	 if (strncmpi(str, 'sc', 2)== 1)

Line 25 and 26 run the color thresholding script written for Photoshop designated by the path where the script is found on the hard drive. Line 29 sets the active layer in Photoshop to the layer with the name ‘NbC’, and Line 30 sends this layer to Matlab and converts the image into a binary format. Lines 31-32 import the SiRich binary image.

pstext = ['var run_EPMA = File ( "D:\\University\\EPMA\\Data\\2011-11Nov-30\\EPMAColorThreshold.jsx");'...
psresult = psjavascriptu(pstext);


Line 35 calculates the binary area of the ‘NbC’ map. Lines 36-40 separate the interdendritic niobium rich phases from the intradendritic phases, where line 36 defines a 2 × 2 structuring element and dilates the ‘NbC’ map with it in line 37. Line 38 closes the dilated images by the same element to fill any holes in the connected components. An area threshold in Line 39 is then set to make a mask of the interdendritic area of the microstructure, which is then multiplied with the original ‘NbC’ map in line 40 resulting in only the interdendritic precipitates remaining. Lines 41-45 convert the image back to a proper format which is then passed back to Photoshop as a new layer with a name specified in line 44.

Anb = bwarea(nb);
se = strel(ones(2,2));
I2 = imdilate(nb,se);
I2 = imclose(I2,se);
inter = bwareaopen(I2,25);
inter = immultiply(nb,inter);
p = bw2rgb(inter);
p = im2uint8(p);
pstext = 'app.activeDocument.layers[0].name="NbC_Gphase_Inter";';

Lines 47-52 multiply the interdendritic area with the silicon map which is then sent back to Photoshop. The area fraction for G-phase is then calculated from the silicon map, and the NbC fraction is found by subtracting the interdendritic silicon map from the interdendritic niobium map. The M23C6 area fraction is simply calculated from the chromium binary map produced from the Photoshop script.

p = bw2rgb(B);
p = im2uint8(p);
pstext = 'app.activeDocument.layers[0].name="G-Phase";';

The area fractions are then sent to populate and array which is then sent to a database configured in Matlab. The Photoshop file is saved and closed, where a new .psd file in opened in the loop.

	%populate two area arrays, name array, and file array
	data(i,1) = {str};
	data(i,2) = {name};        
	data(i,3) = {interNbC};
	data(i,4) = {M23C6Tot};
	data(i,5) = {GPhaseInter};

	pstextsave = [';'];
	psresult = psjavascriptu(pstextsave);            
%export to database
conn = database('EPMADB', '', '');
colnames = {'Sample_Name', 'Sample_Number', 'NbC_Inter', 'M23C6_Tot', 'GPhase_Inter'};
get(conn, 'AutoCommit');
fastinsert(conn, 'EPMA_Area_Table', colnames, data);

1.5. Python phase fraction script for backscattered images

For the MetalTek samples area fraction and precipitate size calculations were done using the OpenCV computer vision library with Python. Instructions on how to install OpenCV and connect its library to Python are provided at The current section will go through the script process_36D located in the Python/MetalTek Scripts directory found at Many functions defined in this script will not be explained past a general understanding of what they do, however, more insight into function from the OpenCV library can be found at The script is organized with the main function at he bottom of the document with all of its internal functions ordered above it. The script is specifically set to handle BSE micrographs from the Tescan Vega-3 SEM, and would need to be modified to work on other SEM micrographs.

After the image is loaded on line 306, a function was written to extract the scale bar that is located in the same position on every micrograph. This functions used Hough Lines to find the length of the scale bar, and OCR to extract the text to find out the real length of the scale bar. A ratio is then defined to convert all pixel measurements into their actual size. The excess parts of the image are then cropped out leaving the micrograph, where its overall area is calculated to be used for area fraction calculations later.

im ="overall2.tif")
scale = scaleExtract(im)
scaleRatio = scale[0]
scaleAmt = scale[2]
img = cropImg(im)
pil = cv2pil(img)
imdem = pil.size
area = imdem[0]*imdem[1]*(scaleRatio**2)

Lines 319-326 will threshold the greyscale image between the values specified in line 319, and then find the contours of each connected component, where the area fraction and the size of the NbC precipitates can then be determined. The contrast threshold is predetermined by the user and must be done for each phase and each micrograph, as there is no way in fully automating this process. An area threshold is done on lines 321 and 326 for any components with an area greater than 2px. This is done to eliminate any noise that has been picked up after thresholding. ’contourWrapFill’ differs from ’contourWrap’ in that the contour boundaries are filled disregarding any holes in the components. Two dilation masks are produced in lines 323 and 324 which will be used by the subsequent phases to determine the interdendritic regions of the microstructure.

nbC = binArray(pil, 255, 124) # ****** Change Binary Threshold
nbC= np2cv(nbC)
[nbC, dump] = contourWrapFill(nbC, 2, high, color= False)
mask = dilateImg(nbC, 3, 3, 2)
cv.SaveImage('mask.jpg', mask)
mask2 = dilateImg(nbC, 4, 4, 4)
cv.SaveImage('mask2.jpg', mask2)    
[nbPhase, dump] = contourWrap(nbC, 2, high, color= False) # '-1' means no area threshold ***** Inter Array

The mask in line 323 is used for the TiC phase to separate it out from the chromium carbide phases. TiC are generally seen to be situated inside of NbC precipitates as multiphase aggregated carbides with a TiC composition at their center where niobium is progressively exchanged with titanium as you reach the rim of the precipitates. The smaller mask defined in line 323 will separate the TiC from the chrome carbides that are typically seen adjacent to the NbC precipitates. Line 333 subtracts the masked TiC image from the NbC image to get rid of any overlapping regions. An area threshold is then conducted between 0-200 pixels.

TiCBefore = cv.CreateImage(cv.GetSize(img), 8, 1)
cv.Mul(carbides, mask, TiCBefore) # get inter M23C6 from G-phase inter mask
TiC = cv.CreateImage(cv.GetSize(img), 8, 1)
cv.Sub(TiCBefore, nbC, TiC) # subtract any G-phase overlap with M23C6
[TiCphase, dump] = contourWrapFill(TiC, 0, 200, color= False) # area <= 10 removed
#Subtract TiC From NbC

The resulting TiC image is then subtracted from the original carbide image leaving behind the chromium carbides in line 348. The intradendritic carbides are then separated from the interdendritic carbides, and labeled as M7C3 and M23C6 respectively. Both carbides undergo a size threshold from 0-200 pixels. The size and area fraction data is then output as both a csv file, and a histogram plot like that shown in Figure 1. An example .zip file of the resulting output after image analysis can be downloaded at

M7 = cv.CreateImage(cv.GetSize(img), 8, 1)
cv.Sub(carbides, TiCphase, M7)
[M7C3, M7C3Arr] = contourWrapFill(M7, 0,  200, color= True)
cv.SaveImage('M7C3.jpg', M7C3)
M23 = cv.CreateImage(cv.GetSize(img), 8, 1)
cv.Mul(M7, mask2, M23)
[M23C6, M23C6Arr] = contourWrapFill(M23, 0,  200, color= True)
cv.SaveImage('M23C6.jpg', M23C6)  

Fig. 1: Python histogram output of NbC precipitate size in design treatment c (commerical 2032Nb alloy, 1” wall thickness, and homogenized) aged after 2 months.

2. References

[1]    Jasiński, W.. Solutionizing of IN 519 superalloy after long-term exposure. Adv Manuf Sci Technol 2008;32(1):69–80.

[2]    Balk, T., Nychka, J.. Microscopy lessons. Tech. Rep.; University of Kentucky; 2006.

[3]    Shinozaki, K., Kuroki, H., Nakao, Y., Nishimoto, K., Inui, M., Takahashi, M.. Deterioration of weldability of long term aged HP heat-resistant cast steel containing Nb, Mo and W. Quarterly Journal of the Japan Welding Society 1998;16(2):223–232.

[4]    Berghof-Hasselbcher, E., Gawenda, P., Schorr, M., Schtze, M., Hoffman, J.. Atlas of Microstructures. Materials Technology Institute; 2008.

[5]    Werner, W.. Interaction of electron beams with matter. 2012.

[6]    Landon, P., Caligiuri, R., Duletsky, P.. The influence of the M23(C,N)6 compound on the mechanical properties of type 422 stainless steel. Metall Mater Trans A 1983;14(7):1395–1408.

[7]    Piekarski, B.. Effect of Nb and Ti additions on microstructure, and identification of precipitates in stabilized Ni-Cr cast austenitic steels. Mater Charact 2001;47(3-4):181–186.

[8]    Zhang, T., Wang, Q., Song, X., Fu, H.. Effect of electromagnetic centrifugal casting on solidification microstructure of cast high speed steel roll. Materialwissenschaft und Werkstofftechnik 2011;42(6):557–561.

[9]    Wolczynski, W., Krajewski, W., Ebner, R., Kloch, J.. The use of equilibrium phase diagram for the calculation of non-equilibrium precipitates in dendritic solidification. theory. Calphad 2001;25(3):401–408.

[10]    Shi, S., Lippold, J.. Microstructure evolution during service exposure of two cast, heat-resisting stainless steels – HP-Nb modified and 20-32Nb. Mater Charact 2008;59(8):1029–1040.

[11]    Gonzalez, R.C., Woods, R.E.. Digital Image Processing. Addison-Wesley Longman Publishing Co., Inc.; 1992.

[12]    Suzuki, S., be, K.. Topological structural analysis of digitized binary images by border following. Computer Vision, Graphics, and Image Processing 1985;30(1):32–46.


3. Glossary *

A specific level of a factor in a factorial design.. 3

4. Acronyms *

Auger Electron Microscope. 4, 9, 10
Backscattered Electron. 7
Energy Dispersive X-ray Spectroscopy. 6, 8–10
Electron Probe Microanalysis. 4, 10
Secondary Electron. 6, 7
Scanning Electron Microscope. 4, 6, 9
Wavelength Dispersive Spectroscopy. 9, 10
X-ray Diffraction. 4, 10


© 2012 CCWJ. All Rights Reserved.