Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
HeartGeometryInformation.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "HeartGeometryInformation.hpp"
37
38#include <cmath>
39#include <fstream>
40#include <sstream>
41#include "OutputFileHandler.hpp"
42#include "Exception.hpp"
43#include "PetscTools.hpp"
44
45// Area of the septum considered to belong to the each ventricle (relative to 1)
46template<unsigned SPACE_DIM>
48
49template<unsigned SPACE_DIM>
51
52template<unsigned SPACE_DIM>
54 const std::string& rEpiFile,
55 const std::string& rEndoFile,
56 bool indexFromZero)
57 : mpMesh(&rMesh)
58{
60
61 // Get nodes defining each surface
62 GetNodesAtSurface(rEpiFile, mEpiSurface, indexFromZero);
63 GetNodesAtSurface(rEndoFile, mEndoSurface, indexFromZero);
64
65 // Compute the distance map of each surface
69}
70
71template<unsigned SPACE_DIM>
73 const std::string& rEpiFile,
74 const std::string& rLVFile,
75 const std::string& rRVFile,
76 bool indexFromZero)
77 : mpMesh(&rMesh)
78{
80
81 // Get nodes defining each surface and compute the distance map of each surface
82
83 GetNodesAtSurface(rEpiFile, mEpiSurface, indexFromZero);
84
85 if (rLVFile != "")
86 {
87 GetNodesAtSurface(rLVFile, mLVSurface, indexFromZero);
88 }
89 else
90 {
91 if (rRVFile == "")
92 {
93 EXCEPTION("At least one of left ventricle or right ventricle files is required");
94 }
95 }
96 if (rRVFile != "")
97 {
98 GetNodesAtSurface(rRVFile, mRVSurface, indexFromZero);
99 }
103
105}
106
107template<unsigned SPACE_DIM>
109{
110 mpMesh = NULL;
111 std::ifstream heterogeneity_file;
112 heterogeneity_file.open(nodeHeterogeneityFileName.c_str());
113
114 if (!(heterogeneity_file.is_open()))
115 {
116 heterogeneity_file.close();
117 EXCEPTION("Could not open heterogeneities file (" << nodeHeterogeneityFileName << ")");
118 }
119
120 while (!heterogeneity_file.eof())
121 {
122 int value;
123
124 heterogeneity_file >> value;
125
126 // Format error (for example read a double), or value not equal to 0, 1, or 2.
127 if ((heterogeneity_file.fail() && !heterogeneity_file.eof()) || value < 0 || value > 2)
128 {
129 heterogeneity_file.close();
130 EXCEPTION("A value in the heterogeneities file (" << nodeHeterogeneityFileName
131 << ") is out of range (or not an integer). It should be epi = 0, mid = 1, endo = 2");
132 }
133
134 if (!heterogeneity_file.eof())
135 {
136 if (value==0)
137 {
138 mLayerForEachNode.push_back(EPI);
139 }
140 else if (value==1)
141 {
142 mLayerForEachNode.push_back(MID);
143 }
144 else
145 {
146 assert(value==2);
147 mLayerForEachNode.push_back(ENDO);
148 }
149 }
150 }
151
152 heterogeneity_file.close();
153}
154
155template<unsigned SPACE_DIM>
157 const std::string& rLineFromFile,
158 std::set<unsigned>& rSurfaceNodeIndexSet,
159 unsigned offset) const
160{
161 std::stringstream line_stream(rLineFromFile);
162 while (!line_stream.eof())
163 {
164 unsigned item;
165 line_stream >> item;
166 // If offset==1 then shift the nodes, since we are assuming MEMFEM format (numbered from 1 on)
167 if (item == 0 && offset != 0)
168 {
169 EXCEPTION("Error when reading surface file. It was assumed not to be indexed from zero, but zero appeared in the list.");
170 }
171 rSurfaceNodeIndexSet.insert(item-offset);
172 }
173}
174
175template<unsigned SPACE_DIM>
177 const std::string& rSurfaceFileName,
178 std::vector<unsigned>& rSurfaceNodes,
179 bool indexFromZero) const
180{
181 // Open the file defining the surface
182 std::ifstream file_stream;
183 unsigned offset=0;
184 if (indexFromZero == false)
185 {
186 offset=1;
187 }
188
189 file_stream.open(rSurfaceFileName.c_str());
190 if (!file_stream.is_open())
191 {
192 EXCEPTION("Wrong surface definition file name " + rSurfaceFileName);
193 }
194
195 // Temporary storage for the nodes, helps discarding repeated values
196 std::set<unsigned> surface_original_node_index_set;
197
198 // Loop over all the triangles and add node indexes to the set
199 std::string line;
200 getline(file_stream, line);
201 do
202 {
203 ProcessLine(line, surface_original_node_index_set, offset);
204
205 getline(file_stream, line);
206 }
207 while (!file_stream.eof());
208 file_stream.close();
209
210 // Make vector big enough
211 rSurfaceNodes.reserve(surface_original_node_index_set.size());
212
213 if (mpMesh->rGetNodePermutation().empty())
214 {
215 // Copy the node indexes from the set to the vector as they are
216 for (std::set<unsigned>::iterator node_index_it=surface_original_node_index_set.begin();
217 node_index_it != surface_original_node_index_set.end();
218 node_index_it++)
219 {
220 rSurfaceNodes.push_back(*node_index_it);
221 }
222 }
223 else
224 {
225 // Copy the original node indices from the set to the vector applying the permutation
226 for (std::set<unsigned>::iterator node_index_it=surface_original_node_index_set.begin();
227 node_index_it != surface_original_node_index_set.end();
228 node_index_it++)
229 {
230 rSurfaceNodes.push_back(mpMesh->rGetNodePermutation()[*node_index_it]);
231 }
232 }
233}
234
235template<unsigned SPACE_DIM>
237{
238
239 if (mDistMapRightVentricle[nodeIndex] >= mDistMapEpicardium[nodeIndex] &&
240 mDistMapRightVentricle[nodeIndex] >= mDistMapLeftVentricle[nodeIndex])
241 {
242 return LEFT_VENTRICLE_WALL;
243 }
244
245 if (mDistMapLeftVentricle[nodeIndex] >= mDistMapEpicardium[nodeIndex] &&
246 mDistMapLeftVentricle[nodeIndex] >= mDistMapRightVentricle[nodeIndex])
247 {
248 return RIGHT_VENTRICLE_WALL;
249 }
250
251 if (mDistMapEpicardium[nodeIndex] >= mDistMapLeftVentricle[nodeIndex] &&
252 mDistMapEpicardium[nodeIndex] >= mDistMapRightVentricle[nodeIndex])
253 {
254 if (mDistMapLeftVentricle[nodeIndex]
255 < LEFT_SEPTUM_SIZE*(mDistMapLeftVentricle[nodeIndex] + mDistMapRightVentricle[nodeIndex]))
256 {
257 return LEFT_SEPTUM;
258 }
259 else
260 {
261 return RIGHT_SEPTUM;
262 }
263 }
264
266 return UNKNOWN; // LCOV_EXCL_LINE
267}
268
269template<unsigned SPACE_DIM>
271{
272 // General case where you provide 3 surfaces: LV, RV, epicardium
273 if (mNumberOfSurfacesProvided == 3)
274 {
275 HeartRegionType node_region = GetHeartRegion(nodeIndex);
276 switch(node_region)
277 {
278 case LEFT_VENTRICLE_WALL:
279 case LEFT_VENTRICLE_SURFACE:
280 return mDistMapLeftVentricle[nodeIndex];
281 break;
282
283 case RIGHT_VENTRICLE_WALL:
284 case RIGHT_VENTRICLE_SURFACE:
285 return mDistMapRightVentricle[nodeIndex];
286 break;
287
288 case LEFT_SEPTUM:
289 return mDistMapLeftVentricle[nodeIndex];
290 break;
291
292 case RIGHT_SEPTUM:
293 return mDistMapRightVentricle[nodeIndex] ;
294 break;
295
296 // LCOV_EXCL_START
297 case UNKNOWN:
298 std::cerr << "Wrong distances node: " << nodeIndex << "\t"
299 << "Epi " << mDistMapEpicardium[nodeIndex] << "\t"
300 << "RV " << mDistMapRightVentricle[nodeIndex] << "\t"
301 << "LV " << mDistMapLeftVentricle[nodeIndex]
302 << std::endl;
303
304 // Make wall_thickness=0 as in Martin's code
305 return 0.0;
306 break;
307 // LCOV_EXCL_STOP
308
309 default:
311 }
312 }
313 // Simplified case where you only provide epi and endo surface definitions
314 else
315 {
316 return mDistMapEndocardium[nodeIndex];
317 }
318
319 // gcc wants to see a return statement at the end of the method.
321 return 0.0; // LCOV_EXCL_LINE
322}
323
324template<unsigned SPACE_DIM>
326{
327 return mDistMapEpicardium[nodeIndex];
328}
329
330template<unsigned SPACE_DIM>
332{
333
334 double dist_endo = GetDistanceToEndo(nodeIndex);
335 double dist_epi = GetDistanceToEpi(nodeIndex);
336 assert( (dist_endo + dist_epi) != 0 );
337 /*
338 * A node contained on both epicardium and lv (or rv) surfaces has wall thickness 0/0.
339 * By setting its value to 0 we consider it contained only on the lv (or rv) surface.
340 */
341 double relative_position = dist_endo / (dist_endo + dist_epi);
342 return relative_position;
343}
344
345template<unsigned SPACE_DIM>
346void HeartGeometryInformation<SPACE_DIM>::DetermineLayerForEachNode(double epiFraction, double endoFraction)
347{
348 if (epiFraction+endoFraction>1)
349 {
350 EXCEPTION("The sum of fractions of epicardial and endocardial layers must be lesser than 1");
351 }
352
353 if ((endoFraction<0) || (epiFraction<0))
354 {
355 EXCEPTION("A fraction of a layer must be positive");
356 }
357
358 mLayerForEachNode.resize(mpMesh->GetNumNodes());
359 for (unsigned i=0; i<mpMesh->GetNumNodes(); i++)
360 {
361 double position = CalculateRelativeWallPosition(i);
362 if (position<endoFraction)
363 {
364 mLayerForEachNode[i] = ENDO;
365 }
366 else if (position<(1-epiFraction))
367 {
368 mLayerForEachNode[i] = MID;
369 }
370 else
371 {
372 mLayerForEachNode[i] = EPI;
373 }
374 }
375}
376
377template<unsigned SPACE_DIM>
378void HeartGeometryInformation<SPACE_DIM>::WriteLayerForEachNode(std::string outputDir, std::string file)
379{
380 OutputFileHandler handler(outputDir,false);
382 {
383 out_stream p_file = handler.OpenOutputFile(file);
384
385 assert(mLayerForEachNode.size()>0);
386 for (unsigned i=0; i<mpMesh->GetNumNodes(); i++)
387 {
388 if (mLayerForEachNode[i]==EPI)
389 {
390 *p_file << "0\n";
391 }
392 else if (mLayerForEachNode[i]==MID)
393 {
394 *p_file << "1\n";
395 }
396 else // endo
397 {
398 *p_file << "2\n";
399 }
400 }
401
402 p_file->close();
403 }
404 PetscTools::Barrier("HeartGeometryInformation::WriteLayerForEachNode"); // Make other processes wait until we're done
405}
406
407template<unsigned SPACE_DIM>
409 const std::vector<unsigned>& rSurfaceNodes)
410{
411
412 assert(rSurfaceNodes.size()>0);
413 //Set min to DBL_MAX etc.
414 c_vector<double, SPACE_DIM> my_minimum_point = scalar_vector<double>(SPACE_DIM, DBL_MAX);
415 //Set max to -DBL_MAX etc.
416 c_vector<double, SPACE_DIM> my_maximum_point=-my_minimum_point;
417
418 //Iterate through the set of points on the surface
419 for (unsigned surface_index=0; surface_index<rSurfaceNodes.size(); surface_index++)
420 {
421 unsigned global_index=rSurfaceNodes[surface_index];
422 if (mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(global_index) )
423 {
424 const c_vector<double, SPACE_DIM>& r_position = mpMesh->GetNode(global_index)->rGetLocation();
425 //Update max/min
426 for (unsigned i=0; i<SPACE_DIM; i++)
427 {
428 if (r_position[i] < my_minimum_point[i])
429 {
430 my_minimum_point[i] = r_position[i];
431 }
432 if (r_position[i] > my_maximum_point[i])
433 {
434 my_maximum_point[i] = r_position[i];
435 }
436 }
437 }
438 }
439
440 //Share the local data and reduce over all processes
441 c_vector<double, SPACE_DIM> global_minimum_point;
442 c_vector<double, SPACE_DIM> global_maximum_point;
443 MPI_Allreduce(&my_minimum_point[0], &global_minimum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MIN, PETSC_COMM_WORLD);
444 MPI_Allreduce(&my_maximum_point[0], &global_maximum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
445
446 ChastePoint<SPACE_DIM> min(global_minimum_point);
447 ChastePoint<SPACE_DIM> max(global_maximum_point);
448
449 return ChasteCuboid<SPACE_DIM>(min, max);
450}
451
452// Explicit instantiation
453template class HeartGeometryInformation<2>;
454template class HeartGeometryInformation<3>;
#define EXCEPTION(message)
#define NEVER_REACHED
void ComputeDistanceMap(const std::vector< unsigned > &rSourceNodeIndices, std::vector< double > &rNodeDistances)
HeartGeometryInformation(AbstractTetrahedralMesh< SPACE_DIM, SPACE_DIM > &rMesh, const std::string &rEpiFile, const std::string &rEndoFile, bool indexFromZero)
AbstractTetrahedralMesh< SPACE_DIM, SPACE_DIM > * mpMesh
double CalculateRelativeWallPosition(unsigned nodeIndex)
void ProcessLine(const std::string &rLineFromFile, std::set< unsigned > &rSurfaceNodeIndexSet, unsigned offset) const
std::vector< unsigned > mEpiSurface
std::vector< unsigned > mLVSurface
ChasteCuboid< SPACE_DIM > CalculateBoundingBoxOfSurface(const std::vector< unsigned > &rSurfaceNodes)
std::vector< double > mDistMapLeftVentricle
double GetDistanceToEndo(unsigned nodeIndex)
HeartRegionType GetHeartRegion(unsigned nodeIndex) const
std::vector< double > mDistMapRightVentricle
double GetDistanceToEpi(unsigned nodeIndex)
void GetNodesAtSurface(const std::string &rSurfaceFileName, std::vector< unsigned > &rSurfaceNodes, bool indexFromZero=true) const
void DetermineLayerForEachNode(double epiFraction, double endoFraction)
std::vector< unsigned > mEndoSurface
std::vector< double > mDistMapEpicardium
void WriteLayerForEachNode(std::string outputDir, std::string file)
std::vector< double > mDistMapEndocardium
std::vector< unsigned > mRVSurface
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static bool AmMaster()
static void Barrier(const std::string callerId="")