Visualizing B1 field using Matlab

% clear all;
% 
% path = 'D:\work\matlab_work\B1-field\Run0001\MultiPoint_Planar_Sensor_0_';
% %path = 'D:\work\matlab_work\B1-field\25\MultiPoint_XY_0_'
% 
% infoFile    = strcat(path,'info.bin');
% geomFile     = strcat(path,'geom.bin');
% 
% dataFileExr = strcat(path,'steady_Bxr_total_freq0.bin');
% dataFileExi = strcat(path,'steady_Bxi_total_freq0.bin');
% 
% dataFileEyr = strcat(path,'steady_Byr_total_freq0.bin');
% dataFileEyi = strcat(path,'steady_Byi_total_freq0.bin');
% 
% dataFileEzr = strcat(path,'steady_Bzr_total_freq0.bin');
% dataFileEzi = strcat(path,'steady_Bzi_total_freq0.bin');
% 

%-------------------------- begin '..._info.bin'
fid = fopen(infoFile);

Rmpt = fread(fid, 4, 'char*1');
version = fread( fid, 1, '*uint8', 0, 'a' );

bitmask = fread( fid, 1, '*uint32', 0, 'a' );   

hasTimeDomainScatteredE = bitget( bitmask, 32 );          % 0
hasTimeDomainTotalE = bitget( bitmask, 31 );              % 0 
hasTimeDomainScatteredH = bitget( bitmask, 30 );          % 0
hasTimeDomainTotalH = bitget( bitmask, 29 );              % 0
hasTimeDomainScatteredB = bitget( bitmask, 28 );          % 0
hasTimeDomainTotalB = bitget( bitmask, 27 );              % 0
hasTimeDomainJ = bitget( bitmask, 26 );                   % 0
hasDiscreteFrequencyTotalE = bitget( bitmask, 25 );       % 1 if steady state E is requested
hasDiscreteFrequencyTotalH = bitget( bitmask, 24 );       % 0
hasDiscreteFrequencyJ = bitget( bitmask, 23 );            % 0
hasDiscreteFrequencyTotalB = bitget( bitmask, 22 );       % 0

geoSize = fread( fid, 1, '*uint32', 0, 'a' );

fclose(fid);

%-------------------------- end '..._info.bin'



%-------------------------- begin '..._geom.bin'
fid = fopen( gemFile );

xPos = zeros(1,geoSize);
yPos = zeros(1,geoSize);
zPos = zeros(1,geoSize);

for kk = 1:geoSize; 
    xPos(kk) = fread( fid, 1, '*uint32' );
    yPos(kk) = fread( fid, 1, '*uint32' );
    zPos(kk) = fread( fid, 1, '*uint32' );
end

fclose(fid);

max_xPos = max(xPos);
max_yPos = max(yPos);
max_zPos = max(zPos);

min_xPos = min(xPos);
min_yPos = min(yPos);
min_zPos = min(zPos);

xdim = max_xPos - min_xPos +1;
ydim = max_yPos - min_yPos +1;
zdim = max_zPos - min_zPos +1;

%-------------------------- end '..._geom.bin'


%-------------------------- get data values
fid = fopen( dataFileBxr );
dataExr = fread( fid, geoSize, 'single' );
fclose(fid);

fid = fopen( dataFileBxi );
dataExi = fread( fid, geoSize, 'single' );
fclose(fid);

fid = fopen( dataFileByr );
dataEyr = fread( fid, geoSize, 'single' );
fclose(fid);

fid = fopen( dataFileByi );
dataEyi = fread( fid, geoSize, 'single' );
fclose(fid);

fid = fopen( dataFileBzr );
dataEzr = fread( fid, geoSize, 'single' );
fclose(fid);

fid = fopen( dataFileBzi );
dataEzi = fread( fid, geoSize, 'single' );
fclose(fid);


%-------------------------- 

dataEx = complex(dataExr, dataExi);
dataEy = complex(dataEyr, dataEyi);
dataEz = complex(dataEzr, dataEzi);

data_X = reshape(dataEx, xdim, ydim, zdim);
data_Y = reshape(dataEy, xdim, ydim, zdim);
data_Z = reshape(dataEz, xdim, ydim, zdim);

 %data_X = reshape(dataEx, zdim, ydim, xdim);
 %data_Y = reshape(dataEy, zdim, ydim, xdim);
 %data_Z = reshape(dataEz, zdim, ydim, xdim);

clc; 
close all;
clear all;

% Select the files
[files, path] = uigetfile('*info.*');
infoFile    = strcat(path,files);
cd(path);

[files, path] = uigetfile('*geom.*');
gemFile    = strcat(path,files);

[files, path] = uigetfile('*_B*.*', 'MultiSelect', 'on');

dataFileBxi = strcat(path,files{:,1});
dataFileBxr = strcat(path,files{:,2});

dataFileByi = strcat(path,files{:,3});
dataFileByr = strcat(path,files{:,4});

dataFileBzi = strcat(path,files{:,5});
dataFileBzr = strcat(path,files{:,6});

cd('C:\Users\UOUWIn\Desktop\Matlab code');

% Read Data

[dataBx, dataBy, dataBz]= ReadData(infoFile, gemFile, dataFileBxr, dataFileBxi, dataFileByr, dataFileByi,dataFileBzr, dataFileBzi);


%%
x = abs(dataBx);
y = abs(dataBy);
z = abs(dataBz);

img = sqrt(x.*x + y.*y+z.*z);
img_size = size(img);
clims = [min(min(img)) max(max(img))*1.5E-1];



%% Display Results

h = figure('Name',path);

%%
subplot(3,2,[1 3]);
h_img = imagesc(img, clims);
axis square; colormap jet;


%% ROI : Avg, std %
img_center = img_size /2;
e_position = [img_center 0 0] + [ -40 -40 40 40];

%e = imellipse(gca,[40 40 114 114]);
%e = imellipse(gca,[40 40 40 40]);
e = imellipse(gca, e_position);
%e = imellipse(gca, [20 20 40 40]);


BW = createMask(e, h_img);
Num = sum(sum(BW));
Avg = sum(sum(img.*BW))./Num;
tmp = ((img.*BW)-(Avg.*BW)).*((img.*BW)-(Avg.*BW));
std = sqrt(sum(sum(tmp))/(Num-1));

title(sprintf('ROI = %d\n Avg= %5.4f[uT], std= %5.4f[uT]', Num, ...
    Avg*1E+6, std*1E+6 ));


%%
subplot(3,2,[2 4]); contour(imrotate(img',90), 100); axis square

% ROI : Max, Min
max_img = max(max(img.*BW));
tmp =img+(1-BW);
min_img = min(min(tmp));

title(sprintf('ROI = %d\n Max = %5.4f[uT], Min= %5.4f[uT]', Num, ...
    max_img*1E+6, min_img*1E+6 ));

%%
subplot(3,2, [5 6]);

plot(img(ceil(img_size(1,1)/2),:), '-r');
hold on; grid on;
plot(img(:,ceil(img_size(1,1)/2)), '-b');
legend('x-dir', 'y-dir','Location','NorthEastOutside');


%%

info_A = sprintf('%s\n',path);
info_B = sprintf('Avg[uT]= \t%f\n',Avg*1E+6);
info_C = sprintf('Std[uT]= \t%f\n',std*1E+6);
info_D = sprintf('Max[uT]= \t%f\n',max_img*1E+6);
info_E = sprintf('Min[uT]= \t%f\n',min_img*1E+6);
info_F = sprintf('Std/Avg= \t%f\n',std/Avg);


disp(info_A);
disp(info_B);
disp(info_C);
disp(info_D);
disp(info_E);
disp(info_F);

% % open the file with write info of Exp.
% fid = fopen('~/Expinfo.txt', 'a+');
% fprintf(fid,'\t%f \t%f \t%f \t%f \t%f\n', ...
%     Avg*1E+6, std*1E+6, ...
%     max_img*1E+6, min_img*1E+6, std/Avg);
% fclose(fid);



%%
cd(path);
saveas(h,'Result.jpg');
cd('.\');