00001 ///////////////////////////////////////////////////////////////////////////// 00002 /// 00003 /// \file 00004 /// 00005 /// Homology gradient vector field computation: Algorithm 1, new version, 00006 /// but still without computing the transpose of the projection. 00007 /// 00008 ///////////////////////////////////////////////////////////////////////////// 00009 00010 // Copyright (C) 2009-2011 by Pawel Pilarczyk. 00011 // 00012 // This file is part of my research software package. This is free software: 00013 // you can redistribute it and/or modify it under the terms of the GNU 00014 // General Public License as published by the Free Software Foundation, 00015 // either version 3 of the License, or (at your option) any later version. 00016 // 00017 // This software is distributed in the hope that it will be useful, 00018 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00020 // GNU General Public License for more details. 00021 // 00022 // You should have received a copy of the GNU General Public License 00023 // along with this software; see the file "license.txt". If not, 00024 // please, see <http://www.gnu.org/licenses/>. 00025 00026 // Started on March 24, 2009. Last revision: March 24, 2011. 00027 00028 00029 #ifndef _CHAINCON_HOMGVF1_H_ 00030 #define _CHAINCON_HOMGVF1_H_ 00031 00032 00033 // include some standard C++ header files 00034 #include <istream> 00035 #include <ostream> 00036 00037 // include selected header files from the CHomP library 00038 #include "chomp/system/config.h" 00039 #include "chomp/system/timeused.h" 00040 #include "chomp/system/textfile.h" 00041 #include "chomp/struct/hashsets.h" 00042 00043 // include relevant local header files 00044 #include "chaincon/combchain.h" 00045 #include "chaincon/comblinmap.h" 00046 00047 00048 // -------------------------------------------------- 00049 // ------------------ Homology GVF ------------------ 00050 // -------------------------------------------------- 00051 00052 /// Computes the homology gradient vector field for the given 00053 /// filtered finite cell complex "K". 00054 /// Cells that represent homology generators are appended to the vector "H". 00055 /// The projection map "pi", the inclusion from "H" to the complex "K", 00056 /// and the homology gradient vector field "phi" 00057 /// are assumed to be initially zero and are constructed. 00058 template <class CellT, class CellArray1, class CellArray2> 00059 void homologyGVF1 (const CellArray1 &K, CellArray2 &H, 00060 tCombLinMap<CellT,CellT> &pi, 00061 tCombLinMap<CellT,CellT> &incl, 00062 tCombLinMap<CellT,CellT> &phi) 00063 { 00064 using chomp::homology::sout; 00065 using chomp::homology::scon; 00066 00067 // don't do anything for an empty complex 00068 if (K. empty ()) 00069 return; 00070 00071 // show a message indicating what is being computed 00072 sout << "Computing a homology gradient vector field...\n"; 00073 00074 // start measuring computation time 00075 chomp::homology::timeused compTime; 00076 compTime = 0; 00077 00078 // prepare a set of cells that represent homology generators 00079 chomp::homology::hashedset<CellT> homGener; 00080 00081 // process all the cells in the filter 00082 int_t KSize = K. size (); 00083 for (int_t i = 0; i < KSize; ++ i) 00084 { 00085 // show progress indicator 00086 if (!(i % 17)) 00087 scon << i << " out of " << KSize << ", " << 00088 (i * 100 / KSize) << "% done.\r"; 00089 00090 // retrieve the i-th cell in the filter 00091 const CellT &c = K [i]; 00092 00093 // compute the modified cell 00094 tCombChain<CellT> ch (c); 00095 ch. add (phi (boundary (c))); 00096 00097 // compute the boundary of the modified cell 00098 tCombChain<CellT> bd_ch (boundary (ch)); 00099 00100 // if the boundary is zero 00101 if (bd_ch. empty ()) 00102 { 00103 // add this cell to the set of homology representants 00104 homGener. add (c); 00105 00106 // define the inclusion for this cell 00107 incl. add (c, ch); 00108 00109 // define the projection on this cell 00110 pi. add (c, c); 00111 00112 // skip the rest 00113 continue; 00114 } 00115 00116 // project the boundary onto the homology generators 00117 tCombChain<CellT> d (pi (bd_ch)); 00118 00119 // make sure that this projection is nontrivial 00120 if (d. empty ()) 00121 { 00122 sout << "Debug: c = " << c << ".\n"; 00123 sout << "Debug: bd_c = " << boundary (c) << 00124 ".\n"; 00125 sout << "Debug: ch = " << ch << ".\n"; 00126 sout << "Debug: bd_ch = " << bd_ch << ".\n"; 00127 throw "Empty boundary in the homology algorithm."; 00128 } 00129 00130 // choose any element of this boundary 00131 const CellT &u = d [0]; 00132 00133 // go through all the previously processed cells 00134 for (int_t j = 0; j < i; ++ j) 00135 { 00136 // retrieve the j-th cell 00137 const CellT &cj = K [j]; 00138 00139 // compute the chain pi (cj) 00140 tCombChain<CellT> cjh (pi (cj)); 00141 00142 // if u does not appear in this chain 00143 // then skip the remaining computations for cj 00144 if (!cjh. contains (u)) 00145 continue; 00146 00147 // compute new phi (cj) by adding ch to it 00148 phi. add (cj, ch); 00149 00150 // compute the new image of cj by the projection 00151 pi. add (cj, d); 00152 } 00153 00154 // remove the cell from the set of homology representants 00155 homGener. remove (u); 00156 incl. remove (u); 00157 } 00158 00159 // copy the computed homology generators 00160 int_t genSize = homGener. size (); 00161 for (int_t i = 0; i < genSize; ++ i) 00162 H. push_back (homGener [i]); 00163 00164 // show information on the computation time 00165 sout << "Homology gvf computed in " << compTime << ".\n"; 00166 00167 return; 00168 } /* homologyGVF1 */ 00169 00170 00171 #endif // _CHAINCON_HOMGVF1_H_ 00172
1.7.2