% Author: Owen Mays
% This script generates a number of files for importing breast phantom models into CST with thermal properties and dispersive dielectric properties assigned.
% Many of them (the .vox, .lat and Material.txt files) are interlocking, they all refer to each other.
% Usage: edit this file to select the frequency range, output directory, and phantom class you wish to use.
% In CST: go to Modeling -> Import/Export -> Import -> 3D Models -> Voxel Data and select the .vox file.
% In the import window SET THE MESH PRIORITY TO -1. Otherwise the voxel data will supercede any other material in the domain (including PEC). Also, drag the box in the preview pane on the right to select the whole rectangle.
% Next, go to Home-> Macros-> Import VBA Macro and select the .mcs file created by this script.
% Finally, click Macros -> [macro name]. This will run the macro reassigning all the materials to their correct dispersive dielectric properties.
% The macro also sets the frequency range for the simulation to match the dielectric properties data.
% The imported voxel data will have a Working Coordinate System aligned with 0,0,0 from the matlab data for ease of cross-referencing.
% To convert from matlab i,j,k, 1-indexed, on a 0.5 mm grid, x(mm) = (i-1)/2, etc.

cc=299792458;
mu0=4*pi*1e-7;
eps0=1/(cc^2*mu0);

freqs = linspace(.5,5,451); %0.5-5GHz
w=2*pi*freqs*1e9;
tau=15e-12; % fixed tau for all tissues

% Min and max values for fat and fibro,
% from Shea et al., "Microwave imaging of realistic numerical breast phantoms," 2010. (Single-pole debye model with tau=15ps everywhere.)
prmDelMax=38.1;
prmDelFat=0.141;

prmInfMax = 29.1;
prmInfFat=2.28;

cndStaMax = 1.38;
cndStaFat = 0.0023;

dx=1; %[mm] Specify desired sampling density. Phantom is uniform 0.5 mm. This code will downsample if the specified dx is coarser
dy=dx;
dz=dx;



% Phantom folder - set this to point to your folder with the phantoms.
phantomDir = 'Matlab formatted Phantoms/'
% Output folder
outputDir='./'
% load prmDel

Class=3

switch Class
  case 1
    load(strcat(phantomDir,'ClassIPhantom.mat'),'prmDel','prmInf','cndSta')
    load(strcat(phantomDir,'region_mask_classI.mat'))
    PhantomSpecString='Class1PhantomSpec.lat';
    VoxFileString = 'Class1PhantomVoxFile.vox';
    MaterialFileString='Class1MaterialFile.txt';
    vbaScriptFileString='Class1MaterialRedefinition.mcs';

  case 2
    load(strcat(phantomDir,'ClassIIPhantom.mat'),'prmDel','prmInf','cndSta')
    load(strcat(phantomDir,'region_mask_classII.mat'))
    PhantomSpecString='Class2PhantomSpec.lat';
    VoxFileString = 'Class2PhantomVoxFile.vox';
    MaterialFileString='Class2MaterialFile.txt';
    vbaScriptFileString='Class2MaterialRedefinition.mcs';

  case 3
    load(strcat(phantomDir,'ClassIIIPhantom.mat'),'prmDel','prmInf','cndSta')
    load(strcat(phantomDir,'region_mask_classIII.mat'))
    PhantomSpecString='Class3PhantomSpec.lat';
    VoxFileString = 'Class3PhantomVoxFile.vox';
    MaterialFileString='Class3MaterialFile.txt';
    vbaScriptFileString='Class3MaterialRedefinition.mcs';

  case 4
    load(strcat(phantomDir,'ClassIVPhantom.mat'),'prmDel','prmInf','cndSta')
    load(strcat(phantomDir,'region_mask_classIV.mat'))
    PhantomSpecString='Class4PhantomSpec.lat';
    VoxFileString = 'Class4PhantomVoxFile.vox';
    MaterialFileString='Class4MaterialFile.txt';
    vbaScriptFileString='Class4MaterialRedefinition.mcs';


  otherwise
    error('Invalid class selection')
end

% Flip things so the orientation corresponds to the view given by, for example,
% imagesc(prmDelReference(100,:,:).'), axis xy
% These flips are applying the axis xy to the data before feeding it to CST.
%



prmDelReference=flipdim(prmDel,3);

prmDel = flipdim(prmDel,2);
prmInf = flipdim(prmInf,2);
cndSta = flipdim(cndSta,2);
regions = flipdim(regions,2);


% Native resolution of the matlab phantom is 0.5 mm. Downsampling to desired dx here:

ds=round(dx/0.5);
prmDel=prmDel(1:ds:end,1:ds:end,1:ds:end);
prmInf=prmInf(1:ds:end,1:ds:end,1:ds:end);
cndSta=cndSta(1:ds:end,1:ds:end,1:ds:end);
regions=regions(1:ds:end,1:ds:end,1:ds:end);

levelMask = zeros(size(prmDel));

materialFileID = fopen(strcat(outputDir,MaterialFileString),'wt');
phantomSpecFileID = fopen(strcat(outputDir,PhantomSpecString),'w');
voxFileID = fopen(strcat(outputDir,VoxFileString),'wt');
vbaScriptFileID = fopen(strcat(outputDir,vbaScriptFileString),'wt');

% Write header information
% Note, as of 2015, the header needed to say "Version 2009" or else an error occurred.
fprintf(materialFileID,'### Version 2009 ###\n');
fprintf(materialFileID,'// Tissue \t Num \t Eps \t Mue \t Kappa \t Rho \t K \t HeatCap \t BloodFlow \t MetabolicRate Red \t Green \t Blue\n');
fprintf(materialFileID,'// \t \t \t \t [s/m] \t [kg/m^3] \t [W/mK] \t [kJ/K/kg] \t [W/K/m^3] \t [W/m^3] \t 0-255 (optional) \n');
fprintf(materialFileID,'// \n');

%Put air and skin in the file
fprintf(materialFileID,'AirBG\t0\t1.000000\t1.0\t0.000000\t1.3\t0.025\t1.005\t0\t0\n');
fprintf(materialFileID,'Air\t1\t1.000000\t1.0\t0.000000\t1.3\t0.025\t1.005\t0\t0\n');
fprintf(materialFileID,'Skin\t2\t35.0\t1.0\t0.000000\t1085\t0.397\t3.765\t5929\t1620\t252\t219\t10\n');



regionMaskBackground = 0;
regionMaskInterior = 1;
regionMaskSkin = 2;

lastReservedLevelNum = 2; % reserving 0 for background, 1 for air, and 2 for skin

% Thresholds for deltaEps, percentage fibroglandular

thresholdPercentages = linspace(0,100,100)/100;
thresholdCenters = (1-thresholdPercentages)*prmDelFat+(thresholdPercentages)*prmDelMax;

% Thermal properties taken from Lazebnik et al. "3D computational study of non-invasive patient-specific microwave hyperthermia treatment of breast cancer," PMB 2010. DOI 10.1088/0031-9155/55/13/003
eps=1.0;
rhoAdipose=1069;
rhoFibro = 1050;
KAdipose=0.306;
KFibro = 0.5;
heatCapAdipose=2.279;
heatCapFibro = 3.6;
bloodFlowAdipose=2229;
bloodFlowFibro = 2700;
metabolicRateAdipose=350;
metabolicRateFibro = 690;


fprintf(vbaScriptFileID,''' %s\n\nSub Main ()\n\n',vbaScriptFileString);
% Write frequencies into VBA script.
fprintf(vbaScriptFileID,'Solver.FrequencyRange "%2.2f","%2.2f"\n',min(freqs),max(freqs));
for regionNum = 1:length(thresholdCenters)
    tissName=sprintf('%1.2fFibro',thresholdPercentages(regionNum));
    tissNum=regionNum+lastReservedLevelNum;
    if(regionNum==1) % first region includes 0 to the border between regions 1 and 2
        lowerBound=0;
    else
        lowerBound=(thresholdCenters(regionNum-1)+thresholdCenters(regionNum))/2;
    end
    if(regionNum==length(thresholdCenters)) % last region spans border between n-1 and n up to max value of prmDel
        upperBound = prmDelMax;
    else
        upperBound = (thresholdCenters(regionNum)+thresholdCenters(regionNum+1))/2;
    end

    levelMask((prmDel>=lowerBound)&(prmDel<=upperBound)) = tissNum;

    redVal=round(255*(thresholdPercentages(regionNum))); % fat is blue, fibro is red. Everything else is "grayscale" between these two values

    greenVal=0;
    blueVal=round(255*(1-thresholdPercentages(regionNum)));

    rho = (thresholdPercentages(regionNum)*rhoFibro+(1-thresholdPercentages(regionNum))*rhoAdipose);
    K = (thresholdPercentages(regionNum)*KFibro+(1-thresholdPercentages(regionNum))*KAdipose);
    heatCap = (thresholdPercentages(regionNum)*heatCapFibro+(1-thresholdPercentages(regionNum))*heatCapAdipose);
    bloodFlow = (thresholdPercentages(regionNum)*bloodFlowFibro+(1-thresholdPercentages(regionNum))*bloodFlowAdipose);
    metabolicRate = (thresholdPercentages(regionNum)*metabolicRateFibro+(1-thresholdPercentages(regionNum))*metabolicRateAdipose);


    prmInfRegion = thresholdPercentages(regionNum)*prmInfMax+(1-thresholdPercentages(regionNum))*prmInfFat;
    prmDelRegion = thresholdPercentages(regionNum)*prmDelMax+(1-thresholdPercentages(regionNum))*prmDelFat;
    cndStaRegion = thresholdPercentages(regionNum)*cndStaMax+(1-thresholdPercentages(regionNum))*cndStaFat;


    eps_c = prmInfRegion+prmDelRegion./(1+1j*w*tau)+cndStaRegion./(1j*w*eps0);



    % Generate VBA file to re-define material properties in CST.
    %%%%%%%%%%%%%% Begin VBA File
    materialNameString=sprintf('%1.2fFibro',thresholdPercentages(regionNum)); %e.g. '0.29Fibro'

    fprintf(vbaScriptFileID,'''@ define material: Voxel Data/%s\n\n', materialNameString);
    fprintf(vbaScriptFileID,strcat('With Material\n', ...
                                   '.Reset\n', ...
                                   '.Name "%s"\n' , ...
                                   '.Folder "Voxel Data"\n' , ...
                                   '.FrqType "all"\n' , ...
                                   '.Type "Normal"\n' , ...
                                   '.MaterialUnit "Frequency", "GHz"\n' , ...
                                   '.MaterialUnit "Geometry", "mm"\n' , ...
                                   '.MaterialUnit "Time", "s"\n' , ...
                                   '.MaterialUnit "Temperature", "Celsius"\n' , ...
                                   '.Epsilon "1"\n' , ...
                                   '.Mue "1"\n' , ...
                                   '.Sigma "0"\n' , ...
                                   '.TanD "0.0"\n' , ...
                                   '.TanDFreq "0.0"\n' , ...
                                   '.TanDGiven "False"\n' , ...
                                   '.TanDModel "ConstTanD"\n' , ...
                                   '.EnableUserConstTanDModelOrderEps "True"\n' , ...
                                   '.ConstTanDModelOrderEps "1"\n' , ...
                                   '.SetElParametricConductivity "False"\n' , ...
                                   '.ReferenceCoordSystem "Global"\n' , ...
                                   '.CoordSystemType "Cartesian"\n' , ...
                                   '.SigmaM "0"\n' , ...
                                   '.TanDM "0.0"\n' , ...
                                   '.TanDMFreq "0.0"\n' , ...
                                   '.TanDMGiven "False"\n' , ...
                                   '.TanDMModel "ConstTanD"\n' , ...
                                   '.EnableUserConstTanDModelOrderMue "True"\n' , ...
                                   '.ConstTanDModelOrderMue "1"\n' , ...
                                   '.SetMagParametricConductivity "False"\n' , ...
                                   '.DispModelEps "None"\n' , ...
                                   '.DispModelMue "None"\n' , ...
                                   '.DispersiveFittingSchemeEps "Nth Order"\n' , ...
                                   '.MaximalOrderNthModelFitEps "10"\n' , ...
                                   '.ErrorLimitNthModelFitEps "0.001"\n' , ...
                                   '.UseOnlyDataInSimFreqRangeNthModelEps "False"\n' , ...
                                   '.DispersiveFittingSchemeMue "1st Order"\n' , ...
                                   '.DispersiveFittingFormatEps "Real_Imag"\n' , ...
                                   '\n'),materialNameString);

    for i=1:length(freqs) %Add entries for CST's dispersion curve fitting.
        fprintf(vbaScriptFileID,'.AddDispersionFittingValueEps "%2.3f", "%3.11f", "%3.11f", "1.0"\n',freqs(i),real(eps_c(i)),-1*imag(eps_c(i)));
    end
    %     .AddDispersionFittingValueEps "5.0e+005", "18.87595921469", "11.5217097589375", "1.0"

    fprintf(vbaScriptFileID,strcat('.UseGeneralDispersionEps "True"\n',...
                                   '.UseGeneralDispersionMue "False"\n',...
                                   '.NLAnisotropy "False"\n',...
                                   '.NLAStackingFactor "1"\n',...
                                   '.NLADirectionX "1"\n',...
                                   '.NLADirectionY "0"\n',...
                                   '.NLADirectionZ "0"\n',...
                                   '.Rho "%4.2f"\n',...
                                   '.ThermalType "Normal"\n',...
                                   '.ThermalConductivity "%2.3f"\n',...
                                   '.HeatCapacity "%2.3f"\n',...
                                   '.MetabolicRate "%3.2f"\n',...
                                   '.BloodFlow "%4.2f"\n',...
                                   '.VoxelConvection "0"\n',...
                                   '.MechanicsType "Unused"\n',...
                                   '.Colour "%1.6f", "%1.6f", "%1.6f"\n',...
                                   '.Wireframe "False"\n',...
                                   '.Reflection "False"\n',...
                                   '.Allowoutline "True"\n',...
                                   '.Transparentoutline "False"\n',...
                                   '.Transparency "0"\n',...
                                   '.Create\n',...
                                   'End With\n'),rho,K,heatCap,metabolicRate,bloodFlow,redVal/255,greenVal/255,blueVal/255);

    %%%%%%%%%%%%% Done with VBA file


    fprintf(materialFileID,'%s\t%d\t%2.6f\t1.0\t0.0\t%4.2f\t%1.3f\t%1.3f\t%4.2f\t%3.1f\t%d\t%d\t%d\n',tissName,tissNum,eps,rho,K,heatCap,bloodFlow,metabolicRate,redVal,greenVal,blueVal);
end

% Write VBA file entry for dispersive skin:

prmDelSkin=24.8;
prmInfSkin=15.3;
sigmaSkin=0.741;

eps_c = prmInfSkin+prmDelSkin./(1+1j*w*tau)+sigmaSkin./(1j*w*eps0);

rho = 1085;
K = 0.397;
heatCap = 3.765;
bloodFlow = 5929;
metabolicRate = 1620;

redVal = 252;
greenVal = 219;
blueVal = 10;

%%%%%%%%%%%%%% Begin VBA File
materialNameString='Skin';

fprintf(vbaScriptFileID,'''@ define material: Voxel Data/%s\n\n', materialNameString);
fprintf(vbaScriptFileID,strcat('With Material\n', ...
                               '.Reset\n', ...
                               '.Name "%s"\n' , ...
                               '.Folder "Voxel Data"\n' , ...
                               '.FrqType "all"\n' , ...
                               '.Type "Normal"\n' , ...
                               '.MaterialUnit "Frequency", "GHz"\n' , ...
                               '.MaterialUnit "Geometry", "mm"\n' , ...
                               '.MaterialUnit "Time", "s"\n' , ...
                               '.MaterialUnit "Temperature", "Celsius"\n' , ...
                               '.Epsilon "1"\n' , ...
                               '.Mue "1"\n' , ...
                               '.Sigma "0"\n' , ...
                               '.TanD "0.0"\n' , ...
                               '.TanDFreq "0.0"\n' , ...
                               '.TanDGiven "False"\n' , ...
                               '.TanDModel "ConstTanD"\n' , ...
                               '.EnableUserConstTanDModelOrderEps "True"\n' , ...
                               '.ConstTanDModelOrderEps "1"\n' , ...
                               '.SetElParametricConductivity "False"\n' , ...
                               '.ReferenceCoordSystem "Global"\n' , ...
                               '.CoordSystemType "Cartesian"\n' , ...
                               '.SigmaM "0"\n' , ...
                               '.TanDM "0.0"\n' , ...
                               '.TanDMFreq "0.0"\n' , ...
                               '.TanDMGiven "False"\n' , ...
                               '.TanDMModel "ConstTanD"\n' , ...
                               '.EnableUserConstTanDModelOrderMue "True"\n' , ...
                               '.ConstTanDModelOrderMue "1"\n' , ...
                               '.SetMagParametricConductivity "False"\n' , ...
                               '.DispModelEps "None"\n' , ...
                               '.DispModelMue "None"\n' , ...
                               '.DispersiveFittingSchemeEps "Nth Order"\n' , ...
                               '.MaximalOrderNthModelFitEps "10"\n' , ...
                               '.ErrorLimitNthModelFitEps "0.001"\n' , ...
                               '.UseOnlyDataInSimFreqRangeNthModelEps "False"\n' , ...
                               '.DispersiveFittingSchemeMue "1st Order"\n' , ...
                               '.DispersiveFittingFormatEps "Real_Imag"\n' , ...
                               '\n'),materialNameString);

for i=1:length(freqs) % add all dispersion entries
    fprintf(vbaScriptFileID,'.AddDispersionFittingValueEps "%2.3f", "%3.11f", "%3.11f", "1.0"\n',freqs(i),real(eps_c(i)),-1*imag(eps_c(i)));
end

fprintf(vbaScriptFileID,strcat('.UseGeneralDispersionEps "True"\n',...
                               '.UseGeneralDispersionMue "False"\n',...
                               '.NLAnisotropy "False"\n',...
                               '.NLAStackingFactor "1"\n',...
                               '.NLADirectionX "1"\n',...
                               '.NLADirectionY "0"\n',...
                               '.NLADirectionZ "0"\n',...
                               '.Rho "%4.2f"\n',...
                               '.ThermalType "Normal"\n',...
                               '.ThermalConductivity "%2.3f"\n',...
                               '.HeatCapacity "%2.3f"\n',...
                               '.MetabolicRate "%3.2f"\n',...
                               '.BloodFlow "%4.2f"\n',...
                               '.VoxelConvection "0"\n',...
                               '.MechanicsType "Unused"\n',...
                               '.Colour "%1.6f", "%1.6f", "%1.6f"\n',...
                               '.Wireframe "False"\n',...
                               '.Reflection "False"\n',...
                               '.Allowoutline "True"\n',...
                               '.Transparentoutline "False"\n',...
                               '.Transparency "0"\n',...
                               '.Create\n',...
                               'End With\n'),rho,K,heatCap,metabolicRate,bloodFlow,redVal/255,greenVal/255,blueVal/255);

%%%%%%%%%%%%% Done with VBA file



fprintf(vbaScriptFileID,'End Sub');

fclose(vbaScriptFileID);




levelMask(regions==regionMaskBackground) = 1; %Assigning skin and background to their own "levels"
levelMask(regions==regionMaskSkin) = 2;

fwrite(phantomSpecFileID,levelMask);
fclose(phantomSpecFileID);

fclose(materialFileID);

%%% Build VOX File %%%
[nx,ny,nz] = size(levelMask);

fprintf(voxFileID,'[Version]\n1.0\n\n[Material]\n//f [MHz] filename\n1000\t%s\n[Background]\n0\n[Voxel]\n//type\tnx\tny\tnz\tdx[mm]\tdy[mm]\tdz[mm]\toffset\tfilename\n',MaterialFileString);
fprintf(voxFileID,'char\t%d\t%d\t%d\t%1.1f\t%1.1f\t%1.1f\t%d\t%s\n',nx,ny,nz,dx,dy,dz,0,PhantomSpecString);
fprintf(voxFileID,'[Bitmap]\n\n');

% Set up a working coordinate system for the breast model that agrees with the matlab 0,0,0 in the bottom left corner for the phantom
originX=(-nx/2*dx)*1e-3; % Based on experimentation, it looks like CST interprets these as being in meters.
originY=(-ny/2*dx)*1e-3;
originZ=(-nz/2*dx)*1e-3;

fprintf(voxFileID,'[WCS]\n//name\t originx\t originy\t originz\t ux\t uy\t uz\t nx\t ny\t nz\n');
fprintf(voxFileID,'breastWCS\t%3.6f\t%3.6f\t%3.6f\t%d\t%d\t%d\t%d\t%d\t%d\n',originX,originY,originZ,1,0,0,0,0,1);
fclose(voxFileID);
