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 
00321     double relative_position;
00322 
00323     if ( (dist_endo + dist_epi) != 0 )
00324     {
00325        relative_position = dist_endo / (dist_endo + dist_epi);
00326     }
00327     else
00328     {
00329         /*
00330          *  A node contained on both epicardium and lv (or rv) surfaces has wall thickness 0/0.
00331          *  By setting its value to 0 we consider it contained only on the lv (or rv) surface.
00332          */
00333         relative_position = 0;
00334     }
00335     return relative_position;
00336 }
00337 
00338 template<unsigned SPACE_DIM>
00339 void HeartGeometryInformation<SPACE_DIM>::DetermineLayerForEachNode(double epiFraction, double endoFraction)
00340 {
00341     if (epiFraction+endoFraction>1)
00342     {
00343         EXCEPTION("The sum of fractions of epicardial and endocardial layers must be lesser than 1");
00344     }
00345 
00346     if ((endoFraction<0) || (epiFraction<0))
00347     {
00348         EXCEPTION("A fraction of a layer must be positive");
00349     }
00350 
00351     mLayerForEachNode.resize(mpMesh->GetNumNodes());
00352     for(unsigned i=0; i<mpMesh->GetNumNodes(); i++)
00353     {
00354         double position = CalculateRelativeWallPosition(i);
00355         if (position<endoFraction)
00356         {
00357             mLayerForEachNode[i] = ENDO;
00358         }
00359         else if (position<(1-epiFraction))
00360         {
00361             mLayerForEachNode[i] = MID;
00362         }
00363         else
00364         {
00365             mLayerForEachNode[i] = EPI;
00366         }
00367     }
00368 }
00369 
00370 
00371 template<unsigned SPACE_DIM>
00372 void HeartGeometryInformation<SPACE_DIM>::WriteLayerForEachNode(std::string outputDir, std::string file)
00373 {
00374     OutputFileHandler handler(outputDir,false);
00375     if (PetscTools::AmMaster())
00376     {
00377         out_stream p_file = handler.OpenOutputFile(file);
00378 
00379         assert(mLayerForEachNode.size()>0);
00380         for(unsigned i=0; i<mpMesh->GetNumNodes(); i++)
00381         {
00382             if(mLayerForEachNode[i]==EPI)
00383             {
00384                 *p_file << "0\n";
00385             }
00386             else if(mLayerForEachNode[i]==MID)
00387             {
00388                 *p_file << "1\n";
00389             }
00390             else // endo
00391             {
00392                 *p_file << "2\n";
00393             }
00394         }
00395 
00396         p_file->close();
00397     }
00398     PetscTools::Barrier("HeartGeometryInformation::WriteLayerForEachNode"); // Make other processes wait until we're done
00399 }
00400 
00401 
00402 template<unsigned SPACE_DIM>
00403 ChasteCuboid<SPACE_DIM> HeartGeometryInformation<SPACE_DIM>::CalculateBoundingBoxOfSurface(const std::vector<unsigned>& rSurfaceNodes)
00404 {
00405 
00406     assert(rSurfaceNodes.size()>0);
00407     //Set min to DBL_MAX etc.
00408     c_vector<double, SPACE_DIM> my_minimum_point;
00409     for (unsigned i=0; i<SPACE_DIM; i++)
00410     {
00411         my_minimum_point[i]=DBL_MAX; //Start with max and work down to actual
00412     }
00413     c_vector<double, SPACE_DIM> my_maximum_point=-my_minimum_point;
00414 
00415     //Iterate through the set of points on the surface
00416     for (unsigned surface_index=0; surface_index<rSurfaceNodes.size(); surface_index++)
00417     {
00418         unsigned global_index=rSurfaceNodes[surface_index];
00419         if (mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(global_index) )
00420         {
00421             c_vector<double, SPACE_DIM> position = mpMesh->GetNode(global_index)->rGetLocation();
00422             //Update max/min
00423             for (unsigned i=0; i<SPACE_DIM; i++)
00424             {
00425                 if (position[i] < my_minimum_point[i])
00426                 {
00427                     my_minimum_point[i] = position[i];
00428                 }
00429                 if (position[i] > my_maximum_point[i])
00430                 {
00431                     my_maximum_point[i] = position[i];
00432                 }
00433             }
00434         }
00435     }
00436 
00437     //Share the local data and reduce over all processes
00438     c_vector<double, SPACE_DIM> global_minimum_point;
00439     c_vector<double, SPACE_DIM> global_maximum_point;
00440     MPI_Allreduce(&my_minimum_point[0], &global_minimum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MIN, PETSC_COMM_WORLD);
00441     MPI_Allreduce(&my_maximum_point[0], &global_maximum_point[0], SPACE_DIM, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);
00442 
00443 
00444     ChastePoint<SPACE_DIM> min(global_minimum_point);
00445     ChastePoint<SPACE_DIM> max(global_maximum_point);
00446 
00447     return ChasteCuboid<SPACE_DIM>(min, max);
00448 }
00449 
00450 
00452 // Explicit instantiation
00454 
00455 //template class HeartGeometryInformation<1>;
00456 template class HeartGeometryInformation<2>;
00457 template class HeartGeometryInformation<3>;

Generated by  doxygen 1.6.2