% Read ROI and tmap
roi = readmr('ROI.bxh');
tmap = readmr('Tmap.bxh');
if ~isequal(size(roi.data), size(tmap.data))
  error('roi and tmap must be the same size');
end

% Find the voxels in the ROI 
roiIndex = 1;  % ROI from the IRIS mseg to process
voxels = find(roi.data == roiIndex);

% Find the max and the location of the max in the ROI
[val,loc] = max(tmap.data(voxels));
[x,y,z] = ind2sub(size(roi.data), voxels(loc));
disp(sprintf('Max of %g found at [%d,%d,%d]', val, x, y, z));