1 /*
2 * Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18
19 #include <algorithm>
20 #include "sar_calculation.h"
21 #include "cfloat"
22 #include "array_ops.h"
23 #include "global.h"
24
25 using namespace std;
26
SAR_Calculation()27 SAR_Calculation::SAR_Calculation()
28 {
29 m_Vx_Used = NULL;
30 m_Vx_Valid = NULL;
31 m_DebugLevel = 0;
32 SetAveragingMethod(SIMPLE, true);
33 Reset();
34 }
35
Reset()36 void SAR_Calculation::Reset()
37 {
38 Delete3DArray(m_Vx_Used,m_numLines);
39 m_Vx_Used = NULL;
40 Delete3DArray(m_Vx_Valid,m_numLines);
41 m_Vx_Valid = NULL;
42
43 m_avg_mass = 0;
44 m_numLines[0]=m_numLines[1]=m_numLines[2]=0;
45 m_cellWidth[0]=m_cellWidth[1]=m_cellWidth[2]=NULL;
46
47 m_cell_volume = NULL;
48 m_cell_density = NULL;
49 m_cell_conductivity = NULL;
50 m_E_field = NULL;
51 m_J_field = NULL;
52
53 Delete3DArray(m_Vx_Used,m_numLines);
54 m_Vx_Used = NULL;
55 Delete3DArray(m_Vx_Valid,m_numLines);
56 m_Vx_Valid = NULL;
57 }
58
SetAveragingMethod(string method,bool silent)59 void SAR_Calculation::SetAveragingMethod(string method, bool silent)
60 {
61 if (method.compare("IEEE_C95_3")==0)
62 return SetAveragingMethod(IEEE_C95_3, silent);
63 if (method.compare("IEEE_62704")==0)
64 return SetAveragingMethod(IEEE_62704, silent);
65 if (method.compare("Simple")==0)
66 return SetAveragingMethod(SIMPLE, silent);
67
68 cerr << __func__ << ": Error, """ << method << """ is an unknown averaging method..." << endl;
69 // unknown method, fallback to simple
70 SetAveragingMethod(SIMPLE, false);
71 }
72
SetAveragingMethod(SARAveragingMethod method,bool silent)73 void SAR_Calculation::SetAveragingMethod(SARAveragingMethod method, bool silent)
74 {
75 if (method==IEEE_62704)
76 {
77 m_massTolerance = 0.0001;
78 m_maxMassIterations = 100;
79 m_maxBGRatio = 0.1;
80 m_markPartialAsUsed = false;
81 m_UnusedRelativeVolLimit = 1.05;
82 m_IgnoreFaceValid = false;
83 if (!silent)
84 cerr << __func__ << ": Setting averaging method to IEEE_62704" << endl;
85 return;
86 }
87 else if (method==IEEE_C95_3)
88 {
89 m_massTolerance = 0.05;
90 m_maxMassIterations = 100;
91 m_maxBGRatio = 1;
92 m_markPartialAsUsed = true;
93 m_UnusedRelativeVolLimit = 1;
94 m_IgnoreFaceValid = false;
95 if (!silent)
96 cerr << __func__ << ": Setting averaging method to IEEE_C95_3" << endl;
97 return;
98 }
99 else if (method==SIMPLE)
100 {
101 m_massTolerance = 0.05;
102 m_maxMassIterations = 100;
103 m_maxBGRatio = 1;
104 m_markPartialAsUsed = true;
105 m_UnusedRelativeVolLimit = 1;
106 m_IgnoreFaceValid = true;
107 if (!silent)
108 cerr << __func__ << ": Setting averaging method to SIMPLE" << endl;
109 return;
110 }
111
112 cerr << __func__ << ": Error, unknown averaging method..." << endl;
113 // unknown method, fallback to simple
114 SetAveragingMethod(SIMPLE, false);
115 }
116
SetNumLines(unsigned int numLines[3])117 void SAR_Calculation::SetNumLines(unsigned int numLines[3])
118 {
119 Delete3DArray(m_Vx_Used,m_numLines);
120 m_Vx_Used = NULL;
121 Delete3DArray(m_Vx_Valid,m_numLines);
122 m_Vx_Valid = NULL;
123
124 for (int n=0;n<3;++n)
125 m_numLines[n]=numLines[n];
126 }
127
SetCellWidth(float * cellWidth[3])128 void SAR_Calculation::SetCellWidth(float* cellWidth[3])
129 {
130 for (int n=0;n<3;++n)
131 m_cellWidth[n]=cellWidth[n];
132 }
133
CalcSAR(float *** SAR)134 float*** SAR_Calculation::CalcSAR(float*** SAR)
135 {
136 if (CheckValid()==false)
137 {
138 cerr << "SAR_Calculation::CalcSAR: SAR calculation is invalid due to missing values... Abort..." << endl;
139 return NULL;
140 }
141 if (m_avg_mass<=0)
142 return CalcLocalSAR(SAR);
143 return CalcAveragedSAR(SAR);
144 }
145
CheckValid()146 bool SAR_Calculation::CheckValid()
147 {
148 for (int n=0;n<3;++n)
149 if (m_cellWidth[n]==NULL)
150 return false;
151 if (m_E_field==NULL)
152 return false;
153 if ((m_J_field==NULL) && (m_cell_conductivity==NULL))
154 return false;
155 if (m_cell_density==NULL)
156 return false;
157 if (m_avg_mass<0)
158 return false;
159 return true;
160 }
161
CalcSARPower()162 double SAR_Calculation::CalcSARPower()
163 {
164 if (CheckValid()==false)
165 {
166 cerr << "SAR_Calculation::CalcSARPower: SAR calculation is invalid due to missing values... Abort..." << endl;
167 return 0;
168 }
169 double power=0;
170 unsigned int pos[3];
171 for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
172 {
173 for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
174 {
175 for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
176 {
177 power += CalcLocalPowerDensity(pos)*CellVolume(pos);
178 }
179 }
180 }
181 return power;
182 }
183
CalcLocalPowerDensity(unsigned int pos[3])184 double SAR_Calculation::CalcLocalPowerDensity(unsigned int pos[3])
185 {
186 double l_pow=0;
187 if (m_cell_conductivity==NULL)
188 {
189 l_pow = abs(m_E_field[0][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[0][pos[0]][pos[1]][pos[2]]);
190 l_pow += abs(m_E_field[1][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[1][pos[0]][pos[1]][pos[2]]);
191 l_pow += abs(m_E_field[2][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[2][pos[0]][pos[1]][pos[2]]);
192 }
193 else
194 {
195 l_pow = m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[0][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[0][pos[0]][pos[1]][pos[2]]);
196 l_pow += m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[1][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[1][pos[0]][pos[1]][pos[2]]);
197 l_pow += m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[2][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[2][pos[0]][pos[1]][pos[2]]);
198 }
199 l_pow*=0.5;
200 return l_pow;
201 }
202
CalcLocalSAR(float *** SAR)203 float*** SAR_Calculation::CalcLocalSAR(float*** SAR)
204 {
205 unsigned int pos[3];
206 m_Valid=0;
207 m_Used=0;
208 m_Unused=0;
209 m_AirVoxel=0;
210 for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
211 {
212 for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
213 {
214 for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
215 {
216 if (m_cell_density[pos[0]][pos[1]][pos[2]]>0)
217 {
218 ++m_Valid;
219 SAR[pos[0]][pos[1]][pos[2]] = CalcLocalPowerDensity(pos)/m_cell_density[pos[0]][pos[1]][pos[2]];
220 }
221 else
222 {
223 ++m_AirVoxel;
224 SAR[pos[0]][pos[1]][pos[2]] = 0;
225 }
226 }
227 }
228 }
229 return SAR;
230 }
231
FindFittingCubicalMass(unsigned int pos[3],float box_size,unsigned int start[3],unsigned int stop[3],float partial_start[3],float partial_stop[3],double & mass,double & volume,double & bg_ratio,int disabledFace,bool ignoreFaceValid)232 int SAR_Calculation::FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
233 float partial_start[3], float partial_stop[3], double &mass, double &volume, double &bg_ratio, int disabledFace, bool ignoreFaceValid)
234 {
235 unsigned int mass_iterations = 0;
236 double old_mass=0;
237 double old_box_size=0;
238 bool face_valid;
239 bool mass_valid;
240 bool voxel_valid;
241
242 //iterate over cubical sizes to find fitting volume to mass
243 while (mass_iterations<m_maxMassIterations)
244 {
245 // calculate cubical mass
246 face_valid = GetCubicalMass(pos, box_size/2, start, stop, partial_start, partial_stop, mass, volume, bg_ratio, disabledFace);
247
248 // check if found mass is valid
249 mass_valid = abs(mass-m_avg_mass)<=m_massTolerance*m_avg_mass;
250 voxel_valid = mass_valid && (face_valid==true) && (bg_ratio<m_maxBGRatio);
251
252 if ((face_valid==false) && (mass<m_avg_mass*(1.0-m_massTolerance)) && (ignoreFaceValid==false))
253 {
254 // this is an invalid cube with a too small total mass --> increasing the box will not yield a valid cube
255 return 1;
256 }
257 else if ((face_valid==false || bg_ratio>=m_maxBGRatio) && (mass_valid))
258 {
259 // this is an invalid cube with a valid total mass
260 if (ignoreFaceValid)
261 return 0;
262 return 2;
263 }
264 if (voxel_valid)
265 {
266 // valid cube found
267 return 0;
268 }
269
270 // if no valid or finally invalid cube is found, calculate an alternaive cube size
271 if (mass_iterations==0)
272 {
273 // on first interation, try a relative resize
274 old_box_size=box_size;
275 box_size*=pow(m_avg_mass/mass,1.0/3.0);
276 }
277 else
278 {
279 // on later interations, try a newton approach
280 float new_box_size = box_size - (mass-m_avg_mass)/(mass-old_mass)*(box_size-old_box_size);
281 old_box_size = box_size;
282 box_size = new_box_size;
283 }
284 old_mass=mass;
285
286 ++mass_iterations;
287 }
288
289 // m_maxMassIterations iterations are exhausted...
290 return -1;
291 }
292
GetCubicalMass(unsigned int pos[3],double box_size,unsigned int start[3],unsigned int stop[3],float partial_start[3],float partial_stop[3],double & mass,double & volume,double & bg_ratio,int disabledFace)293 bool SAR_Calculation::GetCubicalMass(unsigned int pos[3], double box_size, unsigned int start[3], unsigned int stop[3],
294 float partial_start[3], float partial_stop[3], double &mass, double &volume, double &bg_ratio, int disabledFace)
295 {
296 if ((box_size<=0) || std::isnan(box_size) || std::isinf(box_size))
297 {
298 cerr << "SAR_Calculation::GetCubicalMass: critical error: invalid averaging box size!! EXIT" << endl;
299 exit(-1);
300 }
301 bool face_valid=true;
302 for (int n=0;n<3;++n)
303 {
304 // check start position
305 start[n]=pos[n];
306 partial_start[n]=1;
307 float dist=m_cellWidth[n][pos[n]]/2;
308 if (disabledFace==2*n)
309 dist=box_size;
310 else
311 {
312 while (dist<box_size)
313 {
314 // box is outside of domain
315 if (start[n]==0)
316 {
317 partial_start[n]=-1;
318 break;
319 }
320 --start[n];
321 dist+=m_cellWidth[n][start[n]];
322 }
323
324 if (partial_start[n]!=-1)
325 { // box ends inside stop[n] voxel
326 partial_start[n]=1-(dist-box_size)/m_cellWidth[n][start[n]];
327 }
328 if ((partial_start[n]!=-1) && (pos[n]==start[n]))
329 partial_start[n]=2*box_size/m_cellWidth[n][start[n]];
330 }
331
332 // check stop position
333 stop[n]=pos[n];
334 partial_stop[n]=1;
335 dist=m_cellWidth[n][pos[n]]/2;
336 if (disabledFace==2*n+1)
337 dist=box_size;
338 else
339 {
340 while (dist<box_size)
341 {
342 // box is outside of domain
343 if (stop[n]==m_numLines[n]-1)
344 {
345 partial_stop[n]=-1;
346 break;
347 }
348 ++stop[n];
349 dist+=m_cellWidth[n][stop[n]];
350 }
351
352 if (partial_stop[n]!=-1)
353 { // box ends inside stop[n] voxel
354 partial_stop[n]=1-(dist-box_size)/m_cellWidth[n][stop[n]];
355 }
356 if ((partial_stop[n]!=-1) && (pos[n]==stop[n]))
357 partial_stop[n]=2*box_size/m_cellWidth[n][stop[n]];
358 }
359 }
360
361 for (int n=0;n<3;++n)
362 {
363 if (partial_start[n]==-1)
364 face_valid=false;
365 if (partial_stop[n]==-1)
366 face_valid=false;
367 }
368
369 mass = 0;
370 volume = 0;
371 double bg_volume=0;
372 double weight[3];
373 unsigned int f_pos[3];
374 bool face_intersect[6] = {false,false,false,false,false,false};
375 for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
376 {
377 weight[0]=1;
378 if (f_pos[0]==start[0])
379 weight[0]*=abs(partial_start[0]);
380 if (f_pos[0]==stop[0])
381 weight[0]*=abs(partial_stop[0]);
382
383 for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
384 {
385 weight[1]=1;
386 if (f_pos[1]==start[1])
387 weight[1]*=abs(partial_start[1]);
388 if (f_pos[1]==stop[1])
389 weight[1]*=abs(partial_stop[1]);
390
391 for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
392 {
393 weight[2]=1;
394 if (f_pos[2]==start[2])
395 weight[2]*=abs(partial_start[2]);
396 if (f_pos[2]==stop[2])
397 weight[2]*=abs(partial_stop[2]);
398
399 mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2];
400 volume += CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
401
402 if (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]==0)
403 bg_volume += CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
404 else
405 {
406 for (int n=0;n<3;++n)
407 {
408 if (start[n]==f_pos[n])
409 face_intersect[2*n]=true;
410 if (stop[n]==f_pos[n])
411 face_intersect[2*n+1]=true;
412 }
413 }
414 }
415 }
416 }
417 //check if all bounds have intersected a material boundary
418 for (int n=0;n<6;++n)
419 face_valid *= face_intersect[n];
420
421 bg_ratio = bg_volume/volume;
422
423 return face_valid;
424 }
425
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)426 float SAR_Calculation::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)
427 {
428 double power_mass=0;
429 double mass=0;
430 double weight[3];
431 unsigned int f_pos[3];
432 for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
433 {
434 weight[0]=1;
435 if (f_pos[0]==start[0])
436 weight[0]*=abs(partial_start[0]);
437 if (f_pos[0]==stop[0])
438 weight[0]*=abs(partial_stop[0]);
439
440 for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
441 {
442 weight[1]=1;
443 if (f_pos[1]==start[1])
444 weight[1]*=abs(partial_start[1]);
445 if (f_pos[1]==stop[1])
446 weight[1]*=abs(partial_stop[1]);
447
448 for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
449 {
450 weight[2]=1;
451 if (f_pos[2]==start[2])
452 weight[2]*=abs(partial_start[2]);
453 if (f_pos[2]==stop[2])
454 weight[2]*=abs(partial_stop[2]);
455
456 if (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]>=0)
457 {
458 mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2];
459 power_mass += CalcLocalPowerDensity(f_pos) * CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
460 }
461 }
462 }
463 }
464 float vx_SAR = power_mass/mass;
465 if (SAR==NULL)
466 return vx_SAR;
467
468 SAR[pos[0]][pos[1]][pos[2]]=vx_SAR;
469
470 if (assignUsed==false)
471 return vx_SAR;
472
473 // assign SAR to all used voxel
474 bool is_partial[3];
475 for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
476 {
477 if ( ((f_pos[0]==start[0]) && (partial_start[0]!=1)) || ((f_pos[0]==stop[0]) && (partial_stop[0]!=1)) )
478 is_partial[0]=true;
479 else
480 is_partial[0]=false;
481
482 for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
483 {
484 if ( ((f_pos[1]==start[1]) && (partial_start[1]!=1)) || ((f_pos[1]==stop[1]) && (partial_stop[1]!=1)) )
485 is_partial[1]=true;
486 else
487 is_partial[1]=false;
488
489 for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
490 {
491 if ( ((f_pos[2]==start[2]) && (partial_start[2]!=1)) || ((f_pos[2]==stop[2]) && (partial_stop[2]!=1)) )
492 is_partial[2]=true;
493 else
494 is_partial[2]=false;
495
496 if ( (!is_partial[0] && !is_partial[1] && !is_partial[2]) || m_markPartialAsUsed)
497 {
498 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))
499 {
500 m_Vx_Used[f_pos[0]][f_pos[1]][f_pos[2]]=true;
501 SAR[f_pos[0]][f_pos[1]][f_pos[2]]=max(SAR[f_pos[0]][f_pos[1]][f_pos[2]], vx_SAR);
502 }
503 }
504 }
505 }
506 }
507 return vx_SAR;
508 }
509
510
CalcAveragedSAR(float *** SAR)511 float*** SAR_Calculation::CalcAveragedSAR(float*** SAR)
512 {
513 unsigned int pos[3];
514 m_Vx_Used = Create3DArray<bool>(m_numLines);
515 m_Vx_Valid = Create3DArray<bool>(m_numLines);
516
517 double voxel_volume;
518 double total_mass;
519 unsigned int start[3];
520 unsigned int stop[3];
521 float partial_start[3];
522 float partial_stop[3];
523 double bg_ratio;
524 int EC=0;
525
526 // debug counter
527 unsigned int cnt_case1=0;
528 unsigned int cnt_case2=0;
529 unsigned int cnt_NoConvergence=0;
530
531 m_Valid=0;
532 m_Used=0;
533 m_Unused=0;
534 m_AirVoxel=0;
535
536 for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
537 {
538 for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
539 {
540 for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
541 {
542 if (m_cell_density[pos[0]][pos[1]][pos[2]]==0)
543 {
544 SAR[pos[0]][pos[1]][pos[2]] = 0;
545 ++m_AirVoxel;
546 continue;
547 }
548
549 // guess an initial box size and find a fitting cube
550 EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop,
551 partial_start, partial_stop, total_mass, voxel_volume, bg_ratio, -1, m_IgnoreFaceValid);
552
553 if (EC==0)
554 {
555 m_Vx_Valid[pos[0]][pos[1]][pos[2]] = true;
556 m_Vx_Used[pos[0]][pos[1]][pos[2]] = true;
557 ++m_Valid;
558 CalcCubicalSAR(SAR, pos, start, stop, partial_start, partial_stop, true);
559 }
560 else if (EC==1)
561 ++cnt_case1;
562 else if (EC==2)
563 ++cnt_case2;
564 else if (EC==-1)
565 ++cnt_NoConvergence;
566 else
567 cerr << "other EC" << EC << endl;
568 }
569 }
570 }
571 if (cnt_NoConvergence>0)
572 {
573 cerr << "SAR_Calculation::CalcAveragedSAR: Warning, for some voxel a valid averaging cube could not be found (no convergence)... " << endl;
574 }
575 if (m_DebugLevel>0)
576 {
577 cerr << "Number of invalid cubes (case 1): " << cnt_case1 << endl;
578 cerr << "Number of invalid cubes (case 2): " << cnt_case2 << endl;
579 cerr << "Number of invalid cubes (failed to converge): " << cnt_NoConvergence << endl;
580 }
581
582 // count all used and unused etc. + special handling of unused voxels!!
583 m_Used=0;
584 m_Unused=0;
585 for (pos[0]=0;pos[0]<m_numLines[0];++pos[0])
586 {
587 for (pos[1]=0;pos[1]<m_numLines[1];++pos[1])
588 {
589 for (pos[2]=0;pos[2]<m_numLines[2];++pos[2])
590 {
591 if (!m_Vx_Valid[pos[0]][pos[1]][pos[2]] && m_Vx_Used[pos[0]][pos[1]][pos[2]])
592 ++m_Used;
593 if ((m_cell_density[pos[0]][pos[1]][pos[2]]>0) && !m_Vx_Valid[pos[0]][pos[1]][pos[2]] && !m_Vx_Used[pos[0]][pos[1]][pos[2]])
594 {
595 ++m_Unused;
596
597 SAR[pos[0]][pos[1]][pos[2]] = 0;
598 double unused_volumes[6];
599 float unused_SAR[6];
600
601 double min_Vol=FLT_MAX;
602
603 // special handling of unused voxels:
604 for (int n=0;n<6;++n)
605 {
606 EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop,
607 partial_start, partial_stop, total_mass, unused_volumes[n], bg_ratio, n, true);
608 if ((EC!=0) && (EC!=2))
609 {
610 // this should not happen
611 cerr << "SAR_Calculation::CalcAveragedSAR: Error handling unused voxels, can't find fitting cubical averaging volume' " << endl;
612 unused_SAR[n]=0;
613 }
614 else
615 {
616 unused_SAR[n]=CalcCubicalSAR(NULL, pos, start, stop, partial_start, partial_stop, false);
617 min_Vol = min(min_Vol,unused_volumes[n]);
618 }
619 }
620 for (int n=0;n<6;++n)
621 {
622 if (unused_volumes[n]<=m_UnusedRelativeVolLimit*min_Vol)
623 SAR[pos[0]][pos[1]][pos[2]] = max(SAR[pos[0]][pos[1]][pos[2]],unused_SAR[n]);
624 }
625 }
626 }
627 }
628 }
629
630 if (m_Valid+m_Used+m_Unused+m_AirVoxel!=m_numLines[0]*m_numLines[1]*m_numLines[2])
631 {
632 cerr << "SAR_Calculation::CalcAveragedSAR: critical error, mismatch in voxel status count... EXIT" << endl;
633 exit(1);
634 }
635
636 if (m_DebugLevel>0)
637 cerr << "SAR_Calculation::CalcAveragedSAR: Stats: Valid=" << m_Valid << " Used=" << m_Used << " Unused=" << m_Unused << " Air-Voxel=" << m_AirVoxel << endl;
638
639 return SAR;
640 }
641
CellVolume(unsigned int pos[3])642 double SAR_Calculation::CellVolume(unsigned int pos[3])
643 {
644 if (m_cell_volume)
645 return m_cell_volume[pos[0]][pos[1]][pos[2]];
646
647 double volume=1;
648 for (int n=0;n<3;++n)
649 volume*=m_cellWidth[n][pos[n]];
650 return volume;
651 }
652
CellMass(unsigned int pos[3])653 double SAR_Calculation::CellMass(unsigned int pos[3])
654 {
655 return m_cell_density[pos[0]][pos[1]][pos[2]]*CellVolume(pos);
656 }
657
658