isosurface_visualiser_proj.cpp

#include <iostream>
#include <cmath>
#include <string>
#include <fstream>
#include <vector>
#include <algorithm>

const int grid = 50;
const double dx = 1.0 / grid;
const double dy = dx;

const double iso_val = 0.5;

const double tol = 0.001;

double x[grid][grid];
double y[grid][grid];
double xv[grid+1][grid+1];
double yv[grid+1][grid+1];
double r[grid][grid];

double f(double x, double y, double r, double xc, double yc){
  return (r - (std::sqrt(std::pow((x-xc),2) + std::pow((y-yc),2))));
}


int main(){

  // 4 (x,y) pairs for each cell
  double corners[4][2];
  double corner_vals[4];
  double cell_vals[4][4];
  double total = 0.0;

  // store location of cell centres
  for (int i=0; i<grid; ++i){
    for (int j=0; j<grid; ++j){
      x[i][j] = i*dx+0.5*dx;
      y[i][j] = j*dy+0.5*dy;
    }
  }

  // store location of cell corners
  for (int i=0; i<grid+1; ++i){
    for (int j=0; j<grid+1; ++j){
      xv[i][j] = i*dx;
      yv[i][j] = j*dy;
    }
  }


  // create empty .vtk file in text mode
  std::ofstream out_file {"output.vtk", std::ios::out};

  // check file opened safely before writing
  if (!out_file){
    std::cerr << "Error creating file" << std::endl;
    return 1;
  }

  // begin by building .vtk header
  out_file << "# vtk DataFile Version 3.0\nFirst example in 2D\nASCII\nDATASET UNSTRUCTURED_GRID\nPOINTS " << grid*grid << " double\n";


  // Store x,y,z coords of cell centres
  for (int i=0; i<grid; ++i){
    for (int j=0; j<grid; ++j){
      out_file << x[i][j] << " " << y[i][j] << " 0.5\n";
    }
  }


  // define data lookup table
  out_file << "POINT_DATA " << grid*grid << "\nSCALARS radius double\nLOOKUP_TABLE default\n";
  
  // assign data to each cell centre
  for (int i=0; i<grid; ++i){
    for (int j=0; j<grid; ++j){
      r[i][j] = f(x[i][j], y[i][j], 0.3, 0.5, 0.5);
      out_file << r[i][j] << std::endl;
    }
  }


  // for each cell
  for (int i=0; i<grid; ++i){
    for (int j=0; j<grid; ++j){
      std::cout << "\n\n\n\nCell at: " << x[i][j] << "  " << y[i][j] << " with value:  " << r[i][j] << std::endl;
      // store the 4 cell corner locations
      corners[0][0] = xv[i][j];
      corners[0][1] = yv[i][j];
      corners[1][0] = xv[i+1][j];
      corners[1][1] = yv[i+1][j];
      corners[2][0] = xv[i][j+1];
      corners[2][1] = yv[i][j+1];
      corners[3][0] = xv[i+1][j+1];
      corners[3][1] = yv[i+1][j+1];

      // for (int x=0; x<4; ++x){
      // 	std::cout << "Corner " << x+1 << "  at: " << corners[x][0] << "   " << corners[x][1] << std::endl;
      // }

      // corner 1
      cell_vals[0][0] = r[i-1][j-1];
      cell_vals[0][1] = r[i][j-1];
      cell_vals[0][2] = r[i-1][j];
      cell_vals[0][3] = r[i][j];

      // corner 2
      cell_vals[1][0] = r[i][j-1];
      cell_vals[1][1] = r[i+1][j-1];
      cell_vals[1][2] = r[i][j];
      cell_vals[1][3] = r[i+1][j];

      // corner 3
      cell_vals[2][0] = r[i-1][j];
      cell_vals[2][1] = r[i][j];
      cell_vals[2][2] = r[i-1][j+1];
      cell_vals[2][3] = r[i][j+1];

      // corner 4
      cell_vals[3][0] = r[i][j];
      cell_vals[3][1] = r[i+1][j];
      cell_vals[3][2] = r[i][j+1];
      cell_vals[3][3] = r[i+1][j+1];
      
      // for each of the corners
      for (int co=0; co<4; ++co){
	total = 0.0;
	std::cout << "\nCorner " << co+1 << std::endl;
	// for each cell around corner
	for (int ce=0; ce<4; ++ce){
	  total += cell_vals[co][ce];
	}
	corner_vals[co] = total / 4;
	std::cout << "Corner 1 value: " << corner_vals[co] << std::endl;
      }


      // check the four edges for intersections
      for (int co=0; co<4; ++co){
	if (corner_vals[co] >= iso_val){
	  for(int inco=0; inco<4; ++inco){
	    if (inco != co){
	      if (corner_vals[inco] <= iso_val){
		// found intersection
	      }
	    }
	  }
	}
	else {
	  // no isosurface here
	}
      }
    }
  }
   
  return 0;
}