Calculate block width error

[wangkang88max]

@AzdiarGazder
Hello, When I use ORtools to calculate block width, using my ebsd data to run the code（example 7） will result in an error. Have you encountered this problem before?
The version I am using is MTEX-5.10.2.
Here is the code I ran

% *********************************************************************
%             ORTools - Example 7
% *********************************************************************
% Reconstructing the prior austenite microstructure from lath martensite
% with the variant graph approach and analyzing variant pairing and block
% width
% The MTEX description for reconstruction is given here:
% https://mtex-toolbox.github.io/MaParentGrainReconstruction.html
% https://www.youtube.com/watch?v=E9IFyFUQNl4&t=166s
% *********************************************************************
% Dr. Azdiar Gazder, 2022, azdiaratuowdotedudotau
% Dr. Frank Niessen, 2022, contactatfniessendotcom
% (Remove "dot" and "at" to make this email address valid)
% *********************************************************************
home; close all; clear variables;
currentFolder;
set(0,'DefaultFigureWindowStyle','normal');
screenPrint('StartUp','ORTools - Example 7 - Variant Graph Approach');

%% Initialize MTEX
% Startup and set some settings
startup_mtex;
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxPlane');
setMTEXpref('FontSize',14);   
setInterp2Tex;

% Default directories - Do not modify
Ini.dataPath = [strrep(pwd,'\','/'),'/data/'];
Ini.cifPath = [Ini.dataPath,'input/cif/'];
Ini.ebsdPath = [Ini.dataPath,'input/ebsd/'];
Ini.texturePath = [Ini.dataPath,'output/texture/'];
Ini.imagePath = [Ini.dataPath,'output/images/'];

%% Import Script for EBSD Data
%
% This script was automatically created by the import wizard. You should
% run the whoole script or parts of it in order to import your data. There
% is no problem in making any changes to this script.

%% Specify Crystal and Specimen Symmetries

% crystal symmetry
CS = {... 
  'notIndexed',...
  crystalSymmetry('m-3m', [3.7 3.7 3.7], 'mineral', 'Iron fcc', 'color', [0.53 0.81 0.98]),...
  crystalSymmetry('m-3m', [2.9 2.9 2.9], 'mineral', 'Iron bcc (old)', 'color', [0.56 0.74 0.56])};

% plotting convention
setMTEXpref('xAxisDirection','east');
setMTEXpref('zAxisDirection','intoPlane');

%% Specify File Names

% path to files
pname = 'D:\kw\0.2 Ni焊接\1';

% which files to be imported
fname = [pname '\项目 1 样品 1 1 分布图数据 8 jiangzao.cpr'];

%% Import the Data

% create an EBSD variable containing the data
ebsd = EBSD.load(fname,CS,'interface','crc',...
  'convertEuler2SpatialReferenceFrame');

%% Compute, filter and smooth grains
screenPrint('SegmentStart','Computing, filtering and smoothing grains');
% Grains are calculated with a 3� threshold
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',3*degree);
% EBSD data in small grains are removed
ebsd(grains(grains.grainSize < 3)) = [];
% Recalculate the grains from the remaining data ...
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'),'angle',3*degree);
% ... and smooth the grain boundaries
grains = smooth(grains,5);

%% Rename and recolor phases 
screenPrint('SegmentStart','Renaming and recoloring phases');
phaseNames = {'Gamma','AlphaP'};
% Rename 'Iron bcc (old)'to 'AlphaP' and 'Iron fcc' to 'Gamma'
ebsd = renamePhases(ebsd,phaseNames);
% Choose your favourite colors
ebsd = recolorPhases(ebsd);

%% Define and refine parent-to-child orientation relationship
screenPrint('SegmentStart','Define and refine parent-to-child OR');
% Define 'Gamma" as the parent and 'AlphaP' as the child phase
job = setParentGrainReconstructor(ebsd,grains,Ini.cifPath);
% Give an initial guess for the OR: Kurdjumov-Sachs ...
KS = orientation.KurdjumovSachs(job.csParent,job.csChild);
job.p2c = KS;
% ... and refine it based on the fit with boundary misorientations
job.calcParent2Child;
% Let us check the disorientation and compare it with K-S and N-W
% (The disorientation is the misfit between the grain misorientations
% and the misorientation of the OR)
NW = orientation.NishiyamaWassermann(job.csParent,job.csChild);
plotHist_OR_misfit(job,[KS,NW],'legend',{'K-S OR','N-W OR'});
% Display information about the OR
ORinfo(job.p2c);
%    - There are 24 martensitic variants
%    - And a ~2.1� disorientation exists from the Nishiyama-Wassermann OR

%% Plotting (with ORTools functions)
screenPrint('SegmentStart','Plotting some ORTools maps');
% Use some of the ORTools functions to visualize the determined OR
% and its relation to the microstructure

% Phase map
plotMap_phases(job,'linewidth',2);
%    - There is no retained austenite (gamma)

% Parent and child IPF maps
plotMap_IPF_p2c(job,vector3d.Z,'linewidth',2);
%    - It is clear that martensite has formed from multiple prior 
%    - austenite grains and that some surface scratches led to bad indexing

% Child-child grain boundary misorientation map
plotMap_gB_c2c(job,'linewidth',2);
%    - Misorientation angles of ~15-50� are not present within prior
%    austenite grains and thus delinitate prior austenite grain
%    boundaries

% Plot a map of the OR boundary disorientation, or misfit
plotMap_gB_misfit(job,'linewidth',2, 'maxColor',5);
%    - By setting a threshold at 5 degrees, the prior austenite grain
%    boundaries are identified by their large misfit with the OR

% Plot a map of the OR boundary probability 
plotMap_gB_prob(job,'linewidth',2);
%   - the same can be visualized by calculating the probability that a
%     boundary belongs to the OR

%% Reconstruct parent microstructure
%   - Reconstruct the microstructure with the variant graph based approach
job.calcVariantGraph('threshold',4*degree,'tolerance',3.5*degree);
%job.calcVariantGraph('threshold',2.5*degree,'tolerance',2.5*degree,'mergeSimilar')
job.clusterVariantGraph;
% ... plot the votes (high values show high certainty)
figure; plot(job.grains,job.votes.prob(:,1))
mtexColorbar
% ... and calculate the parent orientations
job.calcParentFromVote('minProb',0.5)
% Plot the reconstructed parent microstructure
plotMap_IPF_p2c(job,vector3d.Z,'linewidth',1,'parent');


%% Remove badly reconstructed clusters
% In order to reconstruct the remaining parent grains, we can calculate the
% votes for surrounding parent grains by the already reconstructed parent
% grains

% compute the votes
job.calcGBVotes('p2c','reconsiderAll')
% assign parent orientations according to the votes
job.calcParentFromVote
% plot the result
plot(job.parentGrains,job.parentGrains.meanOrientation)

%% Clean reconstructed grains
% Now clean the grains by: 
% - merging grains with similar orientation
job.mergeSimilar('threshold',7.5*degree);
% - merging small inclusions
job.mergeInclusions('maxSize',150);
% This is the cleaned reconstructed parent microstructure
plotMap_IPF_p2c(job,vector3d.Z,'linewidth',2,'parent');

%% Variant analysis
% Now that both the alpha and the associated prior gamma orientations are
% available, variant analysis can be conducted.
% Plot the possible variants in a pole figure
plotPDF_variants(job);
% Then calculate the variant IDs of all alpha grains ...
job.calcVariants;
% ... and plot them
plotMap_variants(job,'grains','linewidth',2); %Grain data
plotMap_variants(job,'linewidth',2); %EBSD data
plotMap_variants(job,'grains','bc','linewidth',2); %Plot as overlay on BC data


% We see that each child grain contains several variants on EBSD level
% For analyzing variant pairing, we need the variants on the EBSD level
% to be reconstructed as grains
%% Variant pairing (block boundary) analysis
[~,maxGrainId] = max(job.grains.area);
[newGrains,newEBSD] = computeVariantGrains(job);
% Compare variant indexing for old and new grains
figure; % Old child grains
plot(grains(job.csChild),grains(job.csChild).meanOrientation);
figure; % New child grains
plot(newGrains(job.csChild),newGrains(job.csChild).meanOrientation);
% We see that all the variant detail of the variantIds in the EBSD map are
% absent at the grain level. This allows us to analyze the boundaries
% between variants.
variantBoundaries_map = plotMap_KSvariantPairs(job,'linewidth',1.5);
% We can also analyze and plot the same for individual prior austenite grains.
% We also return the grain object as "PAG" with the variant data stored in
% PAG.prop
[variantBoundaries_PAG, PAG] = plotMap_KSvariantPairs(job,'parentGrainId',maxGrainId,'linewidth',2);
% In theory, one could use the reindexed grains to redo the parent grain
% reconstruction based on these grains. This does however not lead to a
% significantly better reconstruction in the present dataset.

%% Calculate martensite block widths
% ... for a PAG by specifying a gamma grain id
plotMap_blockWidths(job,'parentGrainId',maxGrainId,'linewidth',1.5);
% ... for the entire map (more statistics, but can be slow for large maps)
plotMap_blockWidths(job,'linewidth',1);

%% Display variant information for a PAG by specifying a gamma grain id
plotStack(job,'parentGrainId',maxGrainId,'linewidth',1.5);
plotStack(job,'grains','noFrame','parentGrainId',maxGrainId,'linewidth',1.5);

%% Save images
saveImage(Ini.imagePath);

%% Display variant information for a PAG by interactively clicking on a gamma grain
grainClick(job,'noScalebar','noFrame'); %Based on EBSD data

%% Do the same at the grain level 
grainClick(job,'grains','noFrame'); %Based on grain data

%% Display variant pairing (block boundary) analysis for a PAG by interactively clicking on a gamma grain
grainClick(job,'variantPairs','noFrame','linewidth',2);

%% Calculate martensite block widths for a PAG by interactively clicking on a gamma grain
grainClick(job,'blockWidth','noFrame','linewidth',2);

%% Refine gamma twins for a PAG by interactively clicking on a gamma grain
grainClick(job,'parentTwins');`

The following is the error situation

Incorrect use of subscript (line 68)
Out of range; index must be a positive integer or logical.

Error grain2d/subsref (line 61)
    ind = subsind(grains,s(1).subs);

Error plotMap_KSvariantPairs (line 95)
        vP_pId = variantGrains(temp_variantPairs_boundary.grainId).prop.packetId;

Error ORTools_exemplar07 (line 198)
variantBoundaries_map = plotMap_KSvariantPairs(job,'linewidth',1.5);`

I would greatly appreciate it if you could help me.

[frankNiessen]

Hi @wangkang88max ,

it is unfortunately not possible to reproduce this error without having access to your data.

Best wishes
Frank

[wangkang88max]

Hi @frankNiessen
The following is my file and code，

项目 1 02ni 720热处理 10 分布图数据 9 - 副本.zip

Code.zip

Best wishes
Wangkang

[maohonghua]

Brother, have you solved your problem? Can you add your email to discuss it?

[AzdiarGazder]

Hi @wangkang88max,

Thank you for reaching out.
My apologies for the late reply - I am currently away on leave and do not have access to my computer.

I shall have a more definitive answer to give you once I'm back and have had a chance to run your script using your data to diagnose the fault.
I'd be grateful for a bit more of your patience in the meantime.

Warm regards,
Azdi

[AzdiarGazder]

Hi @wangkang88max,

Before we begin, I advise you to please download MTEX v5.11.2 and the latest developer version of ORTools.
Both these versions have been extensively tested to ensure they work together without error.

I've done the following:

1. I downloaded your map and re-named the files as "testMap.cpr" and "testMap.crc".

2. I opened this map in OI Channel-5 Tango, and found that you had already processed your map. (Good job! :))
   However, many of your grains are smaller than 10 pixels - so I filtered them out and re-filled the holes they left behind using these settings:
   Under Edit -> Detect Grains, the critical misorientation = 2 deg and the "Allow boundary completion" = 2 deg.
   Under Edit -> Noise Reduction, critical angle deg = 2 deg.
   Following that, I saved the filtered map with the file names "testMap_C.cpr" and "testMap_C.crc".

3. I copied your map files to the folder: ~\ORTools\data\input\ebsd

4. I downloaded your script and found that it was essentially ORTools - Example 7.
   Instead of using your script, I chose to modify ORTools - Example 7 to input your map.
   To do this step, replace this section in ORTools - Example 7:

%% Load data
% Load an MTEX dataset into 'ebsd'
mtexDataset = 'martensite';
screenPrint('SegmentStart',sprintf('Loading MTEX example data ''%s''',mtexDataset));
ebsd = mtexdata(mtexDataset);

with this command:

%% Load data
% Load your map data into 'ebsd'
ebsd = readCPR(Ini.ebsdPath);

When you run this command, a GUI-based window opens and you can double-click your mouse to select the map file you want to load - in this case "testMap_C.cpr".
ORTools takes care of everything - you do not need to manually input any crystal information.

5. After that, to line 186 in ORTools - Example 7, I added the command:

return

This is because all the commands below this line should only ever be run one at a time - never together in a sequence.

6. The command on line 172 in ORTools - Example 7:

variantBoundaries_map = plotMap_KSvariantPairs(job,'linewidth',1.5);

... takes some time to execute on my laptop (about 7 minutes) as your map has 1,149,785 pixels.
However, I can confirm that there is absolutely no error in its execution.
I'm sorry but I simply cannot replicate your issue.

**Lastly, and importantly, if an explanation from the maintainers resolves your issue, please press the "Mark as answer" button on the appropriate response. Your assistance in this regard enables us to close discussions as complete. **

Hope this helps.

Warm regards,
Azdi