- TImage::IndexType idx = it.GetIndex( );
- TScalar vidx = TScalar( it.Get( ) );
- unsigned long iidx = GetIndex( idx, region );
-
- // Neighbors
- TScalar s = TScalar( 0 );
- for( unsigned int d = 0; d < Dim; ++d )
- {
- TImage::IndexType jdx;
- for( int l = -1; l <= 1; l += 2 )
- {
- jdx = idx;
- jdx[ d ] += l;
- if( region.IsInside( jdx ) )
- {
- TScalar vjdx = TScalar( input->GetPixel( jdx ) );
- unsigned long ijdx = GetIndex( jdx, region );
- TScalar v = std::fabs( vidx - vjdx );
- Lt.push_back( TTriplet( iidx, ijdx, v ) );
- if( maxV < v )
- maxV = v;
-
- } // fi
-
- } // rof
-
- } // rof
-
- } // rof
-
- std::vector< TScalar > diag( N, TScalar( 0 ) );
- for( TTriplets::iterator tIt = Lt.begin( ); tIt != Lt.end( ); ++tIt )
- {
- TScalar v = std::exp( -beta * tIt->value( ) / maxV );
- if( v < eps )
- v = eps;
- *tIt = TTriplet( tIt->row( ), tIt->col( ), -v );
- diag[ tIt->col( ) ] += v;
-
- } // rof
- for( unsigned long i = 0; i < diag.size( ); ++i )
- Lt.push_back( TTriplet( i, i, diag[ i ] ) );
- TSparseMatrix L( N, N );
- L.setFromTriplets( Lt.begin( ), Lt.end( ) );
- L.makeCompressed( );
-
- // Boundary
- TTriplets boundaryt;
- // TODO: for( unsigned int i = 0; i < seeds.size( ); ++i )
- std::map< unsigned long, unsigned long >::const_iterator mIt = seeds.begin( );
- for( unsigned long i = 0; mIt != seeds.end( ); ++mIt, ++i )
- boundaryt.push_back( TTriplet( i, mIt->second - 1, TScalar( 1 ) ) );
- TSparseMatrix boundary( seeds.size( ), nLabels );
- boundary.setFromTriplets( boundaryt.begin( ), boundaryt.end( ) );
- boundary.makeCompressed( );
-
- // Compute RHS
- TTriplets RHSt;
- unsigned long x = 0;
- std::map< unsigned long, unsigned long >::const_iterator sIt = seeds.begin( );
- for( ; sIt != seeds.end( ); ++sIt )
- {
- unsigned long y = 0;
- for( unsigned long n = 0; n < N; ++n )
- {
- std::map< unsigned long, unsigned long >::const_iterator nIt =
- seeds.find( n );
- if( nIt == seeds.end( ) )
- RHSt.push_back( TTriplet( y++, x, L.coeff( sIt->first, n ) ) );
-
- } // rof
- x++;
-
- } // rof
- TSparseMatrix RHS( N - seeds.size( ), seeds.size( ) );
- RHS.setFromTriplets( RHSt.begin( ), RHSt.end( ) );
- RHS.makeCompressed( );
-
- // Compute A
- TTriplets At;
- x = 0;
- for( unsigned long m = 0; m < N; ++m )
- {
- if( seeds.find( m ) == seeds.end( ) )
- {
- unsigned long y = 0;
- for( unsigned long n = 0; n < N; ++n )
- if( seeds.find( n ) == seeds.end( ) )
- At.push_back( TTriplet( y++, x, L.coeff( m, n ) ) );
- x++;
-
- } // fi
-
- } // rof
- TSparseMatrix A( N - seeds.size( ), N - seeds.size( ) );
- A.setFromTriplets( At.begin( ), At.end( ) );
- A.makeCompressed( );
-
- // Solve dirichlet problem
- TSparseSolver solver;
- solver.compute( A );
- if( solver.info( ) != Eigen::Success )