HeartGeometryInformation.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "HeartGeometryInformation.hpp"
00030 
00031 #include <cmath>
00032 #include <fstream>
00033 #include <sstream>
00034 #include "OutputFileHandler.hpp"
00035 #include "Exception.hpp"
00036 #include "PetscTools.hpp"
00037 
00038 // Area of the septum considered to belong to the each ventricle (relative to 1)
00039 template<unsigned SPACE_DIM>
00040 const double HeartGeometryInformation<SPACE_DIM>::LEFT_SEPTUM_SIZE = 2.0/3.0;
00041 
00042 template<unsigned SPACE_DIM>
00043 const double HeartGeometryInformation<SPACE_DIM>::RIGHT_SEPTUM_SIZE = 1.0/3.0;
00044 
00045 template<unsigned SPACE_DIM>
00046 HeartGeometryInformation<SPACE_DIM>::HeartGeometryInformation(AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>& rMesh,
00047                                                               const std::string& rEpiFile,
00048                                                               const std::string& rEndoFile,
00049                                                               bool indexFromZero)
00050    : mpMesh(&rMesh)
00051 {
00052     DistanceMapCalculator<SPACE_DIM, SPACE_DIM> distance_calculator(*mpMesh);
00053 
00054     // Get nodes defining each surface
00055     GetNodesAtSurface(rEpiFile, mEpiSurface, indexFromZero);
00056     GetNodesAtSurface(rEndoFile, mEndoSurface, indexFromZero);
00057 
00058     // Compute the distance map of each surface
00059     distance_calculator.ComputeDistanceMap(mEpiSurface, mDistMapEpicardium);
00060     distance_calculator.ComputeDistanceMap(mEndoSurface, mDistMapEndocardium);
00061     mNumberOfSurfacesProvided = 2;
00062 }
00063 
00064 template<unsigned SPACE_DIM>
00065 HeartGeometryInformation<SPACE_DIM>::HeartGeometryInformation (AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>& rMesh,
00066                                                                const std::string& rEpiFile,
00067                                                                const std::string& rLVFile,
00068                                                                const std::string& rRVFile,
00069                                                                bool indexFromZero)
00070     : mpMesh(&rMesh)
00071 {
00072     DistanceMapCalculator<SPACE_DIM, SPACE_DIM> distance_calculator(*mpMesh);
00073 
00074     // Get nodes defining each surface
00075 
00076     GetNodesAtSurface(rEpiFile, mEpiSurface, indexFromZero);
00077     GetNodesAtSurface(rLVFile, mLVSurface, indexFromZero);
00078     GetNodesAtSurface(rRVFile, mRVSurface, indexFromZero);
00079 
00080     // Compute the distance map of each surface
00081     distance_calculator.ComputeDistanceMap(mEpiSurface, mDistMapEpicardium);
00082     distance_calculator.ComputeDistanceMap(mLVSurface, mDistMapLeftVentricle);
00083     distance_calculator.ComputeDistanceMap(mRVSurface, mDistMapRightVentricle);
00084 
00085     mNumberOfSurfacesProvided = 3;
00086 }
00087 
00088 
00089 template<unsigned SPACE_DIM>
00090 HeartGeometryInformation<SPACE_DIM>::HeartGeometryInformation (std::string nodeHeterogeneityFileName)
00091 {
00092     mpMesh = NULL;
00093     std::ifstream heterogeneity_file;
00094     heterogeneity_file.open(nodeHeterogeneityFileName.c_str());
00095 
00096     if (!(heterogeneity_file.is_open()))
00097     {
00098         heterogeneity_file.close();
00099         std::stringstream ss;
00100         ss << "Could not open heterogeneities file (" << nodeHeterogeneityFileName << ")";
00101         EXCEPTION(ss.str());
00102     }
00103 
00104     while(!heterogeneity_file.eof())
00105     {
00106         int value;
00107 
00108         heterogeneity_file >> value;
00109 
00110         // format error (for example read a double), or value not equal to 0, 1, or 2.
00111         if( (heterogeneity_file.fail() && !heterogeneity_file.eof()) || value < 0 || value > 2)
00112         {
00113             heterogeneity_file.close();
00114             std::stringstream ss;
00115             ss << "A value in the heterogeneities file (" << nodeHeterogeneityFileName
00116                << ") is out of range (or not an integer). It should be epi = 0, mid = 1, endo = 2";
00117             EXCEPTION(ss.str());
00118         }
00119 
00120         if(!heterogeneity_file.eof())
00121         {
00122             if(value==0)
00123             {
00124                 mLayerForEachNode.push_back(EPI);
00125             }
00126             else if(value==1)
00127             {
00128                 mLayerForEachNode.push_back(MID);
00129             }
00130             else
00131             {
00132                 assert(value==2);
00133                 mLayerForEachNode.push_back(ENDO);
00134             }
00135         }
00136     }
00137 
00138     heterogeneity_file.close();
00139 }
00140 
00141 template<unsigned SPACE_DIM>
00142 void HeartGeometryInformation<SPACE_DIM>::ProcessLine(
00143         const std::string& line, std::set<unsigned>& rSurfaceNodeIndexSet, unsigned offset) const
00144 {
00145     std::stringstream line_stream(line);
00146     while (!line_stream.eof())
00147     {
00148         unsigned item;
00149         line_stream >> item;
00150         // If offset==1 then shift the nodes, since we are assuming MEMFEM format (numbered from 1 on)
00151         if (item == 0 && offset != 0) 
00152         {
00153             EXCEPTION("Error when reading surface file.  It was assumed not to be indexed from zero, but zero appeared in the list.");
00154         }
00155         rSurfaceNodeIndexSet.insert(item-offset);
00156     }
00157 }
00158 
00159 
00160 template<unsigned SPACE_DIM>
00161 void HeartGeometryInformation<SPACE_DIM>::GetNodesAtSurface(
00162         const std::string& surfaceFile, std::vector<unsigned>& rSurfaceNodes, bool indexFromZero) const
00163 {
00164     // Open the file defining the surface
00165     std::ifstream file_stream;
00166     unsigned offset=0;
00167     if (indexFromZero == false)
00168     {
00169         offset=1;
00170     }
00171 
00172     file_stream.open(surfaceFile.c_str());
00173     if (!file_stream.is_open())
00174     {
00175         EXCEPTION("Wrong surface definition file name " + surfaceFile);
00176     }
00177 
00178     // Temporary storage for the nodes, helps discarding repeated values
00179     std::set<unsigned> surface_original_node_index_set;
00180 
00181     // Loop over all the triangles and add node indexes to the set
00182     std::string line;
00183     getline(file_stream, line);
00184     do
00185     {
00186         ProcessLine(line, surface_original_node_index_set, offset);
00187 
00188         getline(file_stream, line);
00189     }
00190     while(!file_stream.eof());
00191     file_stream.close();
00192 
00193     // Make vector big enough
00194     rSurfaceNodes.reserve(surface_original_node_index_set.size());
00195 
00196     if (mpMesh->rGetNodePermutation().empty())
00197     {
00198         // Copy the node indexes from the set to the vector as they are
00199         for(std::set<unsigned>::iterator node_index_it=surface_original_node_index_set.begin();
00200             node_index_it != surface_original_node_index_set.end();
00201             node_index_it++)
00202         {
00203             rSurfaceNodes.push_back(*node_index_it);
00204         }
00205     }
00206     else
00207     {
00208         // Copy the original node indices from the set to the vector applying the permutation
00209         for(std::set<unsigned>::iterator node_index_it=surface_original_node_index_set.begin();
00210             node_index_it != surface_original_node_index_set.end();
00211             node_index_it++)
00212         {
00213             rSurfaceNodes.push_back(mpMesh->rGetNodePermutation()[*node_index_it]);
00214         }
00215     }
00216 }
00217 
00218 
00219 
00220 template<unsigned SPACE_DIM>
00221 HeartRegionType HeartGeometryInformation<SPACE_DIM>::GetHeartRegion(unsigned nodeIndex) const
00222 {
00223 
00224     if (mDistMapRightVentricle[nodeIndex] >= mDistMapEpicardium[nodeIndex] &&
00225         mDistMapRightVentricle[nodeIndex] >= mDistMapLeftVentricle[nodeIndex])
00226     {
00227         return HeartRegionCode::LEFT_VENTRICLE_WALL;
00228     }
00229 
00230     if (mDistMapLeftVentricle[nodeIndex] >= mDistMapEpicardium[nodeIndex] &&
00231         mDistMapLeftVentricle[nodeIndex] >= mDistMapRightVentricle[nodeIndex])
00232     {
00233         return HeartRegionCode::RIGHT_VENTRICLE_WALL;
00234     }
00235 
00236     if (mDistMapEpicardium[nodeIndex] >= mDistMapLeftVentricle[nodeIndex] &&
00237         mDistMapEpicardium[nodeIndex] >= mDistMapRightVentricle[nodeIndex])
00238     {
00239         if (mDistMapLeftVentricle[nodeIndex]
00240             < LEFT_SEPTUM_SIZE*(mDistMapLeftVentricle[nodeIndex] + mDistMapRightVentricle[nodeIndex]))
00241         {
00242             return HeartRegionCode::LEFT_SEPTUM;
00243         }
00244         else
00245         {
00246             return HeartRegionCode::RIGHT_SEPTUM;
00247         }
00248     }
00249 
00250     return HeartRegionCode::UNKNOWN;
00251 }
00252 
00253 template<unsigned SPACE_DIM>
00254 double HeartGeometryInformation<SPACE_DIM>::GetDistanceToEndo(unsigned nodeIndex)
00255 {
00256     // General case where you provide 3 surfaces: LV, RV, epicardium
00257     if ( mNumberOfSurfacesProvided == 3)
00258     {
00259         HeartRegionType node_region = GetHeartRegion(nodeIndex);
00260         switch(node_region)
00261         {
00262             case HeartRegionCode::LEFT_VENTRICLE_WALL:
00263             case HeartRegionCode::LEFT_VENTRICLE_SURFACE:
00264                 return mDistMapLeftVentricle[nodeIndex];
00265                 break;
00266 
00267             case HeartRegionCode::RIGHT_VENTRICLE_WALL:
00268             case HeartRegionCode::RIGHT_VENTRICLE_SURFACE:
00269                 return mDistMapRightVentricle[nodeIndex];
00270                 break;
00271 
00272             case HeartRegionCode::LEFT_SEPTUM:
00273                 return mDistMapLeftVentricle[nodeIndex];
00274                 break;
00275 
00276             case HeartRegionCode::RIGHT_SEPTUM:
00277                 return mDistMapRightVentricle[nodeIndex] ;
00278                 break;
00279 
00280             case HeartRegionCode::UNKNOWN:
00281                 #define COVERAGE_IGNORE
00282                 std::cerr << "Wrong distances node: " << nodeIndex << "\t"
00283                           << "Epi " << mDistMapEpicardium[nodeIndex] << "\t"
00284                           << "RV " << mDistMapRightVentricle[nodeIndex] << "\t"
00285                           << "LV " << mDistMapLeftVentricle[nodeIndex]
00286                           << std::endl;
00287 
00288                 // Make wall_thickness=0 as in Martin's code
00289                 return 0.0;
00290                 break;
00291                 #undef COVERAGE_IGNORE
00292 
00293             default:
00294                 NEVER_REACHED;
00295         }
00296     }
00297     // Simplified case where you only provide epi and endo surface definitions
00298     else
00299     {
00300         return mDistMapEndocardium[nodeIndex];
00301     }
00302 
00303     // gcc wants to see a return statement at the end of the method.
00304     NEVER_REACHED;
00305     return 0.0;
00306 }
00307 
00308 template<unsigned SPACE_DIM>
00309 double HeartGeometryInformation<SPACE_DIM>::GetDistanceToEpi(unsigned nodeIndex)
00310 {
00311     return mDistMapEpicardium[nodeIndex];
00312 }
00313 
00314 template<unsigned SPACE_DIM>
00315 double HeartGeometryInformation<SPACE_DIM>::CalculateRelativeWallPosition(unsigned nodeIndex)
00316 {
00317 
00318     double dist_endo = GetDistanceToEndo(nodeIndex);
00319     double dist_epi = GetDistanceToEpi(nodeIndex);
00320     assert( (dist_endo + dist_epi) != 0 );
00321     /*
00322      *  A node contained on both epicardium and lv (or rv) surfaces has wall thickness 0/0.
00323      *  By setting its value to 0 we consider it contained only on the lv (or rv) surface.
00324      */
00325     double relative_position = dist_endo / (dist_endo + dist_epi);
00326     return relative_position;
00327 }
00328 
00329 template<unsigned SPACE_DIM>
00330 void HeartGeometryInformation<SPACE_DIM>::DetermineLayerForEachNode(double epiFraction, double endoFraction)
00331 {
00332     if (epiFraction+endoFraction>1)
00333     {
00334         EXCEPTION("The sum of fractions of epicardial and endocardial layers must be lesser than 1");
00335     }
00336 
00337     if ((endoFraction<0) || (epiFraction<0))
00338     {
00339         EXCEPTION("A fraction of a layer must be positive");
00340     }
00341 
00342     mLayerForEachNode.resize(mpMesh->GetNumNodes());
00343     for(unsigned i=0; i<mpMesh->GetNumNodes(); i++)
00344     {
00345         double position = CalculateRelativeWallPosition(i);
00346         if (position<endoFraction)
00347         {
00348             mLayerForEachNode[i] = ENDO;
00349         }
00350         else if (position<(1-epiFraction))
00351         {
00352             mLayerForEachNode[i] = MID;
00353         }
00354         else
00355         {
00356             mLayerForEachNode[i] = EPI;
00357         }
00358     }
00359 }
00360 
00361 
00362 template<unsigned SPACE_DIM>
00363 void HeartGeometryInformation<SPACE_DIM>::WriteLayerForEachNode(std::string outputDir, std::string file)
00364 {
00365     OutputFileHandler handler(outputDir,false);
00366     if (PetscTools::AmMaster())
00367     {
00368         out_stream p_file = handler.OpenOutputFile(file);
00369 
00370         assert(mLayerForEachNode.size()>0);
00371         for(unsigned i=0; i<mpMesh->GetNumNodes(); i++)
00372         {
00373             if(mLayerForEachNode[i]==EPI)
00374             {
00375                 *p_file << "0\n";
00376             }
00377             else if(mLayerForEachNode[i]==MID)
00378             {
00379                 *p_file << "1\n";
00380             }
00381             else // endo
00382             {
00383                 *p_file << "2\n";
00384             }
00385         }
00386 
00387         p_file->close();
00388     }
00389     PetscTools::Barrier("HeartGeometryInformation::WriteLayerForEachNode"); // Make other processes wait until we're done
00390 }
00391 
00392 
00393 template<unsigned SPACE_DIM>
00394 ChasteCuboid<SPACE_DIM> HeartGeometryInformation<SPACE_DIM>::CalculateBoundingBoxOfSurface(const std::vector<unsigned>& rSurfaceNodes)
00395 {
00396 
00397     assert(rSurfaceNodes.size()>0);
00398     //Set min to DBL_MAX etc.
00399     c_vector<double, SPACE_DIM> my_minimum_point = scalar_vector<double>(SPACE_DIM, DBL_MAX);
00400     //Set max to -DBL_MAX etc.
00401     c_vector<double, SPACE_DIM> my_maximum_point=-my_minimum_point;
00402 
00403     //Iterate through the set of points on the surface
00404     for (unsigned surface_index=0; surface_index<rSurfaceNodes.size(); surface_index++)
00405     {
00406         unsigned global_index=rSurfaceNodes[surface_index];
00407         if (mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(global_index) )
00408         {
00409             c_vector<double, SPACE_DIM> position = mpMesh->GetNode(global_index)->rGetLocation();
00410             //Update max/min
00411             for (unsigned i=0; i<SPACE_DIM; i++)
00412             {
00413                 if (position[i] < my_minimum_point[i])
00414                 {
00415                     my_minimum_point[i] = position[i];
00416                 }
00417                 if (position[i] > my_maximum_point[i])
00418                 {
00419                     my_maximum_point[i] = position[i];
00420                 }
00421             }
00422         }
00423     }
00424 
00425     //Share the local data and reduce over all processes
00426     c_vector<double, SPACE_DIM> global_minimum_point;
00427     c_vector<double, SPACE_DIM> global_maximum_point;
00428     MPI_Allreduce(&my_minimum_point[0], &global_minimum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MIN, PETSC_COMM_WORLD);
00429     MPI_Allreduce(&my_maximum_point[0], &global_maximum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
00430 
00431 
00432     ChastePoint<SPACE_DIM> min(global_minimum_point);
00433     ChastePoint<SPACE_DIM> max(global_maximum_point);
00434 
00435     return ChasteCuboid<SPACE_DIM>(min, max);
00436 }
00437 
00438 
00440 // Explicit instantiation
00442 
00443 //template class HeartGeometryInformation<1>;
00444 template class HeartGeometryInformation<2>;
00445 template class HeartGeometryInformation<3>;

Generated on Mon Nov 1 12:35:17 2010 for Chaste by  doxygen 1.5.5