From c536e1f34440622d919decf0cf27fb65292cea66 Mon Sep 17 00:00:00 2001 From: Thorsten Liebig Date: Thu, 29 Nov 2012 16:45:48 +0100 Subject: [PATCH] process fields: new SAR calculation todo: needs much testing and evaluation Signed-off-by: Thorsten Liebig --- Common/processfields.cpp | 18 + Common/processfields.h | 6 +- Common/processfields_sar.cpp | 263 ++++++++++----- Common/processfields_sar.h | 9 + openEMS.pro | 2 + openems.cpp | 12 +- tools/sar_calculation.cpp | 632 +++++++++++++++++++++++++++++++++++ tools/sar_calculation.h | 119 +++++++ 8 files changed, 973 insertions(+), 88 deletions(-) create mode 100644 tools/sar_calculation.cpp create mode 100644 tools/sar_calculation.h diff --git a/Common/processfields.cpp b/Common/processfields.cpp index cff34d3..f52fc14 100644 --- a/Common/processfields.cpp +++ b/Common/processfields.cpp @@ -73,10 +73,28 @@ string ProcessFields::GetFieldNameByType(DumpType type) return "RotH-Field"; case SAR_LOCAL_DUMP: return "SAR-local"; + case SAR_1G_DUMP: + return "SAR_1g"; + case SAR_10G_DUMP: + return "SAR_10g"; + case SAR_RAW_DATA: + return "SAR_raw_data"; } return "unknown field"; } +bool ProcessFields::NeedConductivity() const +{ + switch (m_DumpType) + { + case J_FIELD_DUMP: + return true; + default: + return false; + } + return false; +} + void ProcessFields::InitProcess() { if (Enabled==false) return; diff --git a/Common/processfields.h b/Common/processfields.h index 82eaa21..b628d4f 100644 --- a/Common/processfields.h +++ b/Common/processfields.h @@ -40,7 +40,7 @@ public: Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP), (conduction) electric current density (kappa*E) (J_FIELD_DUMP) and total current density (rotH) */ - enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP, J_FIELD_DUMP, ROTH_FIELD_DUMP, SAR_LOCAL_DUMP}; + enum DumpType { E_FIELD_DUMP=0, H_FIELD_DUMP=1, J_FIELD_DUMP=2, ROTH_FIELD_DUMP=3, SAR_LOCAL_DUMP=20, SAR_1G_DUMP=21, SAR_10G_DUMP=22, SAR_RAW_DATA=29}; virtual string GetProcessingName() const {return "common field processing";} @@ -66,7 +66,7 @@ public: void SetDumpMode2Cell() {SetDumpMode(Engine_Interface_Base::CELL_INTERPOLATE);} //! Set dump type: 0 for E-fields, 1 for H-fields, 2 for D-fields, 3 for B-fields, 4 for J-fields, etc... - void SetDumpType(DumpType type) {m_DumpType=type;} + virtual void SetDumpType(DumpType type) {m_DumpType=type;} double CalcTotalEnergyEstimate() const; @@ -74,6 +74,8 @@ public: static string GetFieldNameByType(DumpType type); + virtual bool NeedConductivity() const; + protected: DumpType m_DumpType; FileType m_fileType; diff --git a/Common/processfields_sar.cpp b/Common/processfields_sar.cpp index 076390c..66428aa 100644 --- a/Common/processfields_sar.cpp +++ b/Common/processfields_sar.cpp @@ -19,37 +19,53 @@ #include "operator_base.h" #include "tools/vtk_file_writer.h" #include "tools/hdf5_file_writer.h" +#include "tools/sar_calculation.h" ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessFieldsFD(eng_if) { + m_UseCellKappa = false; } ProcessFieldsSAR::~ProcessFieldsSAR() { - for (size_t n = 0; n >(numLines)); - m_J_FD_Fields.push_back(Create_N_3DArray >(numLines)); + if (!m_UseCellKappa) + m_J_FD_Fields.push_back(Create_N_3DArray >(numLines)); } } @@ -85,6 +111,9 @@ int ProcessFieldsSAR::Process() double T; FDTD_FLOAT**** field_td=NULL; + //save dump type + DumpType save_dump_type = m_DumpType; + // calc E-field m_DumpType = E_FIELD_DUMP; field_td = CalcField(); @@ -111,32 +140,35 @@ int ProcessFieldsSAR::Process() Delete_N_3DArray(field_td,numLines); // calc J-field - m_DumpType = J_FIELD_DUMP; - field_td = CalcField(); - T = m_Eng_Interface->GetTime(m_dualTime); - for (size_t n = 0; n exp_jwt_2_dt = std::exp( (std::complex)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) ); - exp_jwt_2_dt *= 2; // *2 for single-sided spectrum - exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval - field_fd = m_J_FD_Fields.at(n); - for (pos[0]=0; pos[0]GetTime(m_dualTime); + for (size_t n = 0; n exp_jwt_2_dt = std::exp( (std::complex)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) ); + exp_jwt_2_dt *= 2; // *2 for single-sided spectrum + exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval + field_fd = m_J_FD_Fields.at(n); + for (pos[0]=0; pos[0](field_td,numLines); } - Delete_N_3DArray(field_td,numLines); //reset dump type - m_DumpType = SAR_LOCAL_DUMP; + m_DumpType = save_dump_type; ++m_FD_SampleCount; return GetNextInterval(); @@ -147,85 +179,152 @@ void ProcessFieldsSAR::DumpFDData() if (Enabled==false) return; unsigned int pos[3]; unsigned int orig_pos[3]; - FDTD_FLOAT*** SAR = Create3DArray(numLines); - std::complex**** E_field_fd = NULL; - std::complex**** J_field_fd = NULL; + float*** SAR = Create3DArray(numLines); double coord[3]; - double density; ContinuousStructure* CSX = Op->GetGeometryCSX(); CSProperties* prop = NULL; CSPropMaterial* matProp = NULL; double power; - double l_pow; - for (size_t n = 0; n(numLines); + float*** cell_density = Create3DArray(numLines); + float*** cell_kappa = NULL; + if (m_UseCellKappa) + cell_kappa = Create3DArray(numLines); + + bool found_UnIsotropic=false; + + // calculate volumes and masses for all cells + for (pos[0]=0; pos[0]GetDiscLine(0,orig_pos[0],true); + for (pos[1]=0; pos[1]GetDiscLine(0,orig_pos[0],true); - for (pos[1]=0; pos[1]GetDiscLine(1,orig_pos[1],true); + for (pos[2]=0; pos[2]GetDiscLine(1,orig_pos[1],true); - for (pos[2]=0; pos[2]GetDiscLine(2,orig_pos[2],true); + + cell_volume[pos[0]][pos[1]][pos[2]] = Op->GetCellVolume(orig_pos); + cell_density[pos[0]][pos[1]][pos[2]] = 0.0; + + prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); + if (prop!=0) { - orig_pos[2] = posLines[2][pos[2]]; - coord[2] = Op->GetDiscLine(2,orig_pos[2],true); - l_pow = abs(E_field_fd[0][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[0][pos[0]][pos[1]][pos[2]]); - l_pow += abs(E_field_fd[1][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[1][pos[0]][pos[1]][pos[2]]); - l_pow += abs(E_field_fd[2][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[2][pos[0]][pos[1]][pos[2]]); - - power += 0.5*l_pow*Op->GetCellVolume(orig_pos); - - prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL); - SAR[pos[0]][pos[1]][pos[2]] = 0.0; - density=0.0; - if (prop!=0) + matProp = dynamic_cast(prop); + if (matProp) { - matProp = dynamic_cast(prop); - density = matProp->GetDensityWeighted(coord); - if (density>0) - { - SAR[pos[0]][pos[1]][pos[2]] = l_pow*0.5/density; - } + found_UnIsotropic |= !matProp->GetIsotropy(); + cell_density[pos[0]][pos[1]][pos[2]] = matProp->GetDensityWeighted(coord); + if (m_UseCellKappa) + cell_kappa[pos[0]][pos[1]][pos[2]] = matProp->GetKappaWeighted(0,coord); } } } } + } + if (found_UnIsotropic) + cerr << "ProcessFieldsSAR::DumpFDData(): Warning, found unisotropic material in SAR calculation... this is unsupported!" << endl; - if (m_fileType==VTK_FILETYPE) + float* cellWidth[3]; + for (int n=0;n<3;++n) + { + cellWidth[n]=new float[numLines[n]]; + for (unsigned int i=0;iGetDiscDelta(n,posLines[n][i])*Op->GetGridDelta(); + } + + if (m_DumpType == SAR_RAW_DATA) + { + if (m_fileType!=HDF5_FILETYPE) { - stringstream ss; - ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n); - - m_Vtk_Dump_File->SetFilename(ss.str()); - m_Vtk_Dump_File->ClearAllFields(); - m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR); - if (m_Vtk_Dump_File->Write()==false) - cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl; + cerr << "ProcessFieldsSAR::DumpFDData(): Error, wrong file type, this should not happen!!! skipped" << endl; + return; } - else if (m_fileType==HDF5_FILETYPE) + + size_t datasize[]={numLines[0],numLines[1],numLines[2]}; + for (size_t n = 0; nWriteScalarField(ss.str(), SAR, datasize)==false) - cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl; - float freq[1] = {(float)m_FD_Samples.at(n)}; - if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false) - cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl; - float pow[1] = {(float)power}; - if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false) - cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl; + if (m_HDF5_Dump_File->WriteVectorField(ss.str(), m_E_FD_Fields.at(n), datasize)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; } - else - cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl; - } + m_HDF5_Dump_File->SetCurrentGroup("/CellData"); + if (m_UseCellKappa==false) + cerr << "ProcessFieldsSAR::DumpFDData: Error, cell conductivity data not available, this should not happen... skipping! " << endl; + else if (m_HDF5_Dump_File->WriteScalarField("Conductivity", cell_kappa, datasize)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; + if (m_HDF5_Dump_File->WriteScalarField("Density", cell_density, datasize)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; + if (m_HDF5_Dump_File->WriteScalarField("Volume", cell_volume, datasize)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; + } + else + { + SAR_Calculation SAR_Calc; + SAR_Calc.SetNumLines(numLines); + if (m_DumpType == SAR_LOCAL_DUMP) + SAR_Calc.SetAveragingMass(0); + else if (m_DumpType == SAR_1G_DUMP) + SAR_Calc.SetAveragingMass(1e-3); + else if (m_DumpType == SAR_10G_DUMP) + SAR_Calc.SetAveragingMass(10e-3); + else + { + cerr << "ProcessFieldsSAR::DumpFDData: unknown SAR dump type...!" << endl; + } + SAR_Calc.SetCellDensities(cell_density); + SAR_Calc.SetCellWidth(cellWidth); + SAR_Calc.SetCellVolumes(cell_volume); + SAR_Calc.SetCellCondictivity(cell_kappa); + + for (size_t n = 0; nWriteScalarField(ss.str(), SAR, datasize)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; + float freq[1] = {(float)m_FD_Samples.at(n)}; + if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; + float pow[1] = {(float)power}; + if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false) + cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl; + } + else + cerr << "ProcessFieldsSAR::DumpFDData: unknown File-Type" << endl; + } + } + for (int n=0;n<3;++n) + delete[] cellWidth[n]; + Delete3DArray(cell_volume,numLines); + Delete3DArray(cell_density,numLines); + Delete3DArray(cell_kappa,numLines); Delete3DArray(SAR,numLines); } diff --git a/Common/processfields_sar.h b/Common/processfields_sar.h index cfa69ef..22fd0f6 100644 --- a/Common/processfields_sar.h +++ b/Common/processfields_sar.h @@ -26,6 +26,10 @@ public: ProcessFieldsSAR(Engine_Interface_Base* eng_if); virtual ~ProcessFieldsSAR(); + virtual void SetDumpType(DumpType type); + + virtual bool NeedConductivity() const; + virtual string GetProcessingName() const {return "SAR dump";} virtual void InitProcess(); @@ -36,9 +40,14 @@ public: virtual void SetOptResolution(double optRes, int dir=-1); + //! Set to true for using the conductivity found at the center of a cell, instead of default E*J + virtual void SetUseCellConductivity(bool val) {m_UseCellKappa=val;} + protected: virtual void DumpFDData(); + bool m_UseCellKappa; + //! frequency domain electric field storage vector****> m_E_FD_Fields; //! frequency domain current density storage diff --git a/openEMS.pro b/openEMS.pro index e746718..75cdefd 100644 --- a/openEMS.pro +++ b/openEMS.pro @@ -155,6 +155,7 @@ SOURCES += tools/global.cpp \ tools/array_ops.cpp \ tools/ErrorMsg.cpp \ tools/AdrOp.cpp \ + tools/sar_calculation.cpp \ tools/vtk_file_writer.cpp \ tools/hdf5_file_writer.cpp @@ -224,6 +225,7 @@ HEADERS += tools/ErrorMsg.h \ tools/global.h \ tools/useful.h \ tools/aligned_allocator.h \ + tools/sar_calculation.h \ tools/vtk_file_writer.h \ tools/hdf5_file_writer.h diff --git a/openems.cpp b/openems.cpp index d7e7ab8..b0c2a60 100644 --- a/openems.cpp +++ b/openems.cpp @@ -432,7 +432,7 @@ bool openEMS::SetupProcessing() ProcField = new ProcessFieldsTD(NewEngineInterface()); else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=13)) ProcField = new ProcessFieldsFD(NewEngineInterface()); - else if (db->GetDumpType()==20) + else if ( ((db->GetDumpType()>=20) && (db->GetDumpType()<=22)) || (db->GetDumpType()==29) ) ProcField = new ProcessFieldsSAR(NewEngineInterface()); else cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl; @@ -455,12 +455,16 @@ bool openEMS::SetupProcessing() ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType()); if (db->GetDumpType()==20) - { ProcField->SetDumpType(ProcessFields::SAR_LOCAL_DUMP); - } + if (db->GetDumpType()==21) + ProcField->SetDumpType(ProcessFields::SAR_1G_DUMP); + if (db->GetDumpType()==22) + ProcField->SetDumpType(ProcessFields::SAR_10G_DUMP); + if (db->GetDumpType()==29) + ProcField->SetDumpType(ProcessFields::SAR_RAW_DATA); //SetupMaterialStorages() has previewed storage needs... refresh here to prevent cleanup!!! - if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12) || (db->GetDumpType()==20)) && Enable_Dumps ) + if ( ProcField->NeedConductivity() && Enable_Dumps ) FDTD_Op->SetMaterialStoreFlags(1,true); ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode()); diff --git a/tools/sar_calculation.cpp b/tools/sar_calculation.cpp new file mode 100644 index 0000000..077ca2f --- /dev/null +++ b/tools/sar_calculation.cpp @@ -0,0 +1,632 @@ +/* +* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + + +#include "sar_calculation.h" +#include "cfloat" +#include "array_ops.h" + +SAR_Calculation::SAR_Calculation() +{ + m_Vx_Used = NULL; + m_Vx_Valid = NULL; + m_DebugLevel = 1; + SetAveragingMethod(SIMPLE); + Reset(); +} + +void SAR_Calculation::Reset() +{ + Delete3DArray(m_Vx_Used,m_numLines); + m_Vx_Used = NULL; + Delete3DArray(m_Vx_Valid,m_numLines); + m_Vx_Valid = NULL; + + m_avg_mass = 0; + m_numLines[0]=m_numLines[1]=m_numLines[2]=0; + m_cellWidth[0]=m_cellWidth[1]=m_cellWidth[2]=NULL; + + m_cell_volume = NULL; + m_cell_density = NULL; + m_cell_conductivity = NULL; + m_E_field = NULL; + m_J_field = NULL; + + Delete3DArray(m_Vx_Used,m_numLines); + m_Vx_Used = NULL; + Delete3DArray(m_Vx_Valid,m_numLines); + m_Vx_Valid = NULL; +} + +void SAR_Calculation::SetAveragingMethod(SARAveragingMethod method) +{ + if (method==IEEE_C95_3) + { + m_massTolerance = 0.0001; + m_maxMassIterations = 100; + m_maxBGRatio = 0.1; + m_markPartialAsUsed = false; + m_UnusedRelativeVolLimit = 1.05; + m_IgnoreFaceValid = false; + return; + } + else if (method==IEEE_62704) + { + m_massTolerance = 0.05; + m_maxMassIterations = 100; + m_maxBGRatio = 1; + m_markPartialAsUsed = true; + m_UnusedRelativeVolLimit = 1; + m_IgnoreFaceValid = false; + return; + } + else if (method==SIMPLE) + { + m_massTolerance = 0.05; + m_maxMassIterations = 100; + m_maxBGRatio = 1; + m_markPartialAsUsed = true; + m_UnusedRelativeVolLimit = 1; + m_IgnoreFaceValid = true; + return; + } + + // unknown method, fallback to simple + SetAveragingMethod(SIMPLE); +} + +void SAR_Calculation::SetNumLines(unsigned int numLines[3]) +{ + Delete3DArray(m_Vx_Used,m_numLines); + m_Vx_Used = NULL; + Delete3DArray(m_Vx_Valid,m_numLines); + m_Vx_Valid = NULL; + + for (int n=0;n<3;++n) + m_numLines[n]=numLines[n]; +} + +void SAR_Calculation::SetCellWidth(float* cellWidth[3]) +{ + for (int n=0;n<3;++n) + m_cellWidth[n]=cellWidth[n]; +} + +float*** SAR_Calculation::CalcSAR(float*** SAR) +{ + if (CheckValid()==false) + { + cerr << "SAR_Calculation::CalcSAR: SAR calculation is invalid due to missing values... Abort..." << endl; + return NULL; + } + if (m_avg_mass<=0) + return CalcLocalSAR(SAR); + return CalcAveragedSAR(SAR); +} + +bool SAR_Calculation::CheckValid() +{ + for (int n=0;n<3;++n) + if (m_cellWidth[n]==NULL) + return false; + if (m_E_field==NULL) + return false; + if ((m_J_field==NULL) && (m_cell_conductivity==NULL)) + return false; + if (m_cell_density==NULL) + return false; + if (m_avg_mass<0) + return false; + return true; +} + +double SAR_Calculation::CalcSARPower() +{ + if (CheckValid()==false) + { + cerr << "SAR_Calculation::CalcSARPower: SAR calculation is invalid due to missing values... Abort..." << endl; + return 0; + } + double power=0; + unsigned int pos[3]; + for (pos[0]=0; pos[0]0) + { + ++m_Valid; + SAR[pos[0]][pos[1]][pos[2]] = CalcLocalPowerDensity(pos)/m_cell_density[pos[0]][pos[1]][pos[2]]; + } + else + { + ++m_AirVoxel; + SAR[pos[0]][pos[1]][pos[2]] = 0; + } + } + } + } + return SAR; +} + +int SAR_Calculation::FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3], + float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace, bool ignoreFaceValid) +{ + unsigned int mass_iterations = 0; + float old_mass=0; + float old_box_size=0; + bool face_valid; + bool mass_valid; + bool voxel_valid; + + //iterate over cubical sizes to find fitting volume to mass + while (mass_iterations increasing the box will not yield a valid cube + return 1; + } + else if ((face_valid==false || bg_ratio>=m_maxBGRatio) && (mass_valid)) + { + // this is an invalid cube with a valid total mass + if (ignoreFaceValid) + return 0; + return 2; + } + if (voxel_valid) + { + // valid cube found + return 0; + } + + // if no valid or finally invalid cube is found, calculate an alternaive cube size + if (mass_iterations==0) + { + // on first interation, try a relative resize + old_box_size=box_size; + box_size*=pow(m_avg_mass/mass,1.0/3.0); + } + else + { + // on later interations, try a newton approach + float new_box_size = box_size - (mass-m_avg_mass)/(mass-old_mass)*(box_size-old_box_size); + old_box_size = box_size; + box_size = new_box_size; + } + old_mass=mass; + + ++mass_iterations; + } + + // m_maxMassIterations iterations are exhausted... + return -1; +} + +bool SAR_Calculation::GetCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3], + float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace) +{ + if ((box_size<=0) || isnan(box_size) || isinf(box_size)) + { + cerr << "SAR_Calculation::GetCubicalMass: critical error: invalid averaging box size!! EXIT" << endl; + exit(-1); + } + bool face_valid=true; + for (int n=0;n<3;++n) + { + // check start position + start[n]=pos[n]; + partial_start[n]=1; + float dist=m_cellWidth[n][pos[n]]/2; + if (disabledFace==2*n) + dist=box_size; + else + { + while (dist=0) + { + mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2]; + power_mass += CalcLocalPowerDensity(f_pos) * CellVolume(f_pos)*weight[0]*weight[1]*weight[2]; + } + } + } + } + float vx_SAR = power_mass/mass; + if (SAR==NULL) + return vx_SAR; + + SAR[pos[0]][pos[1]][pos[2]]=vx_SAR; + + if (assignUsed==false) + return vx_SAR; + + // assign SAR to all used voxel + bool is_partial[3]; + for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0]) + { + if ( ((f_pos[0]==start[0]) && (partial_start[0]!=1)) || ((f_pos[0]==stop[0]) && (partial_stop[0]!=1)) ) + is_partial[0]=true; + else + is_partial[0]=false; + + for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1]) + { + if ( ((f_pos[1]==start[1]) && (partial_start[1]!=1)) || ((f_pos[1]==stop[1]) && (partial_stop[1]!=1)) ) + is_partial[1]=true; + else + is_partial[1]=false; + + for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2]) + { + if ( ((f_pos[2]==start[2]) && (partial_start[2]!=1)) || ((f_pos[2]==stop[2]) && (partial_stop[2]!=1)) ) + is_partial[2]=true; + else + is_partial[2]=false; + + if ( (!is_partial[0] && !is_partial[1] && !is_partial[2]) || m_markPartialAsUsed) + { + if (!m_Vx_Valid[f_pos[0]][f_pos[1]][f_pos[2]] && (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]>0)) + { + m_Vx_Used[f_pos[0]][f_pos[1]][f_pos[2]]=true; + SAR[f_pos[0]][f_pos[1]][f_pos[2]]=max(SAR[f_pos[0]][f_pos[1]][f_pos[2]], vx_SAR); + } + } + } + } + } + return vx_SAR; +} + + +float*** SAR_Calculation::CalcAveragedSAR(float*** SAR) +{ + unsigned int pos[3]; + m_Vx_Used = Create3DArray(m_numLines); + m_Vx_Valid = Create3DArray(m_numLines); + + float voxel_volume; + float total_mass; + unsigned int start[3]; + unsigned int stop[3]; + float partial_start[3]; + float partial_stop[3]; + float bg_ratio; + int EC=0; + + // debug counter + unsigned int cnt_case1=0; + unsigned int cnt_case2=0; + unsigned int cnt_NoConvergence=0; + + m_Valid=0; + m_Used=0; + m_Unused=0; + m_AirVoxel=0; + + for (pos[0]=0; pos[0]0) + { + cerr << "SAR_Calculation::CalcAveragedSAR: Warning, for some voxel a valid averaging cube could not be found (no convergence)... " << endl; + } + if (m_DebugLevel>0) + { + cerr << "Number of invalid cubes (case 1): " << cnt_case1 << endl; + cerr << "Number of invalid cubes (case 2): " << cnt_case2 << endl; + cerr << "Number of invalid cubes (falied to converge): " << cnt_NoConvergence << endl; + } + + // count all used and unused etc. + special handling of unused voxels!! + m_Used=0; + m_Unused=0; + for (pos[0]=0;pos[0]0) && !m_Vx_Valid[pos[0]][pos[1]][pos[2]] && !m_Vx_Used[pos[0]][pos[1]][pos[2]]) + { + ++m_Unused; + + SAR[pos[0]][pos[1]][pos[2]] = 0; + float unused_volumes[6]; + float unused_SAR[6]; + + float min_Vol=FLT_MAX; + + // special handling of unused voxels: + for (int n=0;n<6;++n) + { + EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop, + partial_start, partial_stop, total_mass, unused_volumes[n], bg_ratio, n, true); + if ((EC!=0) && (EC!=2)) + { + // this should not happen + cerr << "SAR_Calculation::CalcAveragedSAR: Error handling unused voxels, can't find fitting cubical averaging volume' " << endl; + unused_SAR[n]=0; + } + else + { + unused_SAR[n]=CalcCubicalSAR(NULL, pos, start, stop, partial_start, partial_stop, false); + min_Vol = min(min_Vol,unused_volumes[n]); + } + } + for (int n=0;n<6;++n) + { + if (unused_volumes[n]<=m_UnusedRelativeVolLimit*min_Vol) + SAR[pos[0]][pos[1]][pos[2]] = max(SAR[pos[0]][pos[1]][pos[2]],unused_SAR[n]); + } + } + } + } + } + + if (m_Valid+m_Used+m_Unused+m_AirVoxel!=m_numLines[0]*m_numLines[1]*m_numLines[2]) + { + cerr << "SAR_Calculation::CalcAveragedSAR: critical error, mismatch in voxel status count... EXIT" << endl; + exit(1); + } + + if (m_DebugLevel>0) + cerr << "SAR_Calculation::CalcAveragedSAR: Stats: Valid=" << m_Valid << " Used=" << m_Used << " Unused=" << m_Unused << " Air-Voxel=" << m_AirVoxel << endl; + + return SAR; +} + +float SAR_Calculation::CellVolume(unsigned int pos[3]) +{ + if (m_cell_volume) + return m_cell_volume[pos[0]][pos[1]][pos[2]]; + + float volume=1; + for (int n=0;n<3;++n) + volume*=m_cellWidth[n][pos[n]]; + return volume; +} + +float SAR_Calculation::CellMass(unsigned int pos[3]) +{ + return m_cell_density[pos[0]][pos[1]][pos[2]]*CellVolume(pos); +} + diff --git a/tools/sar_calculation.h b/tools/sar_calculation.h new file mode 100644 index 0000000..7434e21 --- /dev/null +++ b/tools/sar_calculation.h @@ -0,0 +1,119 @@ +/* +* Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de) +* +* This program is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* This program is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with this program. If not, see . +*/ + +#ifndef SAR_CALCULATION_H +#define SAR_CALCULATION_H + +#include + +class SAR_Calculation +{ +public: + SAR_Calculation(); + + enum SARAveragingMethod { IEEE_C95_3, IEEE_62704, SIMPLE}; + + //! Reset and initialize all values (will keep all SAR settings) + void Reset(); + + //! Set the debug level + void SetDebugLevel(int level) {m_DebugLevel=level;} + + //! Set the used averaging method + void SetAveragingMethod(SARAveragingMethod method); + + //! Set number of lines in all direcitions. (mandatory information) + void SetNumLines(unsigned int numLines[3]); + //! Set cell width in all direcitions. (mandatory information for averaging) + void SetCellWidth(float* cellWidth[3]); + + //! Set the averaging mash. (mandatory information for averaging) + void SetAveragingMass(float mass) {m_avg_mass=mass;} + + //! Set the cell volumes (optional for speedup) + void SetCellVolumes(float*** cell_volume) {m_cell_volume=cell_volume;} + + //! Set the cell densities (mandatory information) + void SetCellDensities(float*** cell_density) {m_cell_density=cell_density;} + + //! Set the cell conductivities (mandatory if no current density field is given) + void SetCellCondictivity(float*** cell_conductivity) {m_cell_conductivity=cell_conductivity;} + + //! Set the electric field (mandatory information) + void SetEField(std::complex**** field) {m_E_field=field;} + //! Set the current density field (mandatory if no conductivity distribution is given) + void SetJField(std::complex**** field) {m_J_field=field;} + + //! Calculate the SAR, requires a preallocated 3D array + float*** CalcSAR(float*** SAR); + + //! Calculate the total power dumped + double CalcSARPower(); + +protected: + unsigned int m_numLines[3]; + float* m_cellWidth[3]; + + float m_avg_mass; + float*** m_cell_volume; + float*** m_cell_density; + float*** m_cell_conductivity; + std::complex**** m_E_field; + std::complex**** m_J_field; + + bool*** m_Vx_Used; + bool*** m_Vx_Valid; + + unsigned int m_Valid; + unsigned int m_Used; + unsigned int m_Unused; + unsigned int m_AirVoxel; + + int m_DebugLevel; + + /*********** SAR calculation parameter and settings ***********/ + float m_massTolerance; + unsigned int m_maxMassIterations; + float m_maxBGRatio; + bool m_markPartialAsUsed; + float m_UnusedRelativeVolLimit; + bool m_IgnoreFaceValid; + + /*********** SAR calculations methods ********/ + double CalcLocalPowerDensity(unsigned int pos[3]); + + //! Calculate the local SAR + float*** CalcLocalSAR(float*** SAR); + + /****** start SAR averaging and all necessary methods ********/ + //! Calculate the averaged SAR + float*** CalcAveragedSAR(float*** SAR); + + int FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3], + float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace=-1, bool ignoreFaceValid=false); + bool GetCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3], + float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace=-1); + + float CalcCubicalSAR(float*** SAR, unsigned int pos[3], unsigned int start[3], unsigned int stop[3], float partial_start[3], float partial_stop[3], bool assignUsed=false); + /****** end SAR averaging and all necessary methods ********/ + + bool CheckValid(); + float CellVolume(unsigned int pos[3]); + float CellMass(unsigned int pos[3]); +}; + +#endif // SAR_CALCULATION_H