examples/pl-functions/combustion-vineyard.cpp
author Dmitriy Morozov <dmitriy@mrzv.org>
Tue, 05 Jun 2012 17:53:08 -0700
branchdev
changeset 258 bb5bc5eff779
parent 199 2bde4c56101c
permissions -rw-r--r--
Added viewer.show_complex_2D + center points in the persistence diagram

/*
 * Author: Dmitriy Morozov
 * Department of Computer Science, Duke University, 2005
 */

#include <utilities/log.h>

#include <iostream>
#include <fstream>
#include <algorithm>

#include "grid2D.h"
#include <topology/lsvineyard.h>

const int xsize = 600;
const int ysize = 600;
const int var = 0;          // which variable to use out of nc of them in each record in the file

template<typename T>
void read(std::ifstream& ifs, T& var)
{
    ifs.read(reinterpret_cast<char*>(&var), sizeof(T));
}

int main(int argc, char** argv)
{
    if (argc < 3)
    {
        std::cout << "Usage: combustion-vineyard FRAME1 FRAME2" << std::endl;
        exit(0);
    }
    
    int size0, nc0;
    int size1, nc1;

    std::cout << "Reading: " << argv[1] << std::endl;
    std::ifstream ifs0(argv[1], std::ios::binary);
    std::cout << "Reading: " << argv[2] << std::endl;
    std::ifstream ifs1(argv[2], std::ios::binary);

    if (!ifs0 || !ifs1)
    {
        std::cout << "Could not open the frames" << std::endl;
        exit(0);
    }

    read(ifs0, size0); read(ifs0, nc0);
    read(ifs1, size1); read(ifs1, nc1);

    assert(size0 == size1); assert(nc0 == nc1);
    assert(size0 == xsize*ysize);
    
    Grid2D g0(xsize, ysize), g1(xsize, ysize);
    
    for (int y = 0; y < ysize; ++y)
        for (int x = 0; x < xsize; ++x)
            for (int d = 0; d < nc0; ++d)
            {
                float val0, val1;
                read(ifs0, val0);
                read(ifs1, val1);
                if (d == var)
                {
                    g0(x,y) = val0;
                    g1(x,y) = val1;
                }
            }
    std::cout << "Grids read" << std::endl;
    
    // Generate the complex, initialize the vineyard (which also computes the pairing)
    typedef                     LSVineyard<Grid2D::CoordinateIndex, Grid2D>             Grid2DVineyard;

    Grid2DVineyard::LSFiltration        simplices;    
    g0.complex_generator(make_push_back_functor(simplices));
    Grid2DVineyard::VertexComparison    vcmp(g0);
    Grid2DVineyard::SimplexComparison   scmp(vcmp);
    simplices.sort(scmp);

    Grid2DVineyard              v(g0.begin(), g0.end(), simplices, g0);
    std::cout << "Filtration generated, size: " << v.filtration().size() << std::endl;
    std::cout << "Pairing computed" << std::endl;
    
    // Compute vineyard
    v.compute_vineyard(g1);
    std::cout << "Vineyard computed" << std::endl;

    v.vineyard().save_edges("combustion");
}