added raster, untested

CMakeLists.txt

@@ -58,7 +58,7 @@ polyfem_python_download_data()
 ################################################################################
 # Subdirectories
 ################################################################################
-add_library(polyfempy MODULE src/binding.cpp)
+add_library(polyfempy MODULE src/binding.cpp src/raster.cpp)
 target_link_libraries(polyfempy PRIVATE polyfem)
 target_link_libraries(polyfempy PRIVATE pybind11::module)
 set_target_properties(polyfempy PROPERTIES PREFIX "${PYTHON_MODULE_PREFIX}" SUFFIX "${PYTHON_MODULE_EXTENSION}")

src/binding.cpp

@@ -5,6 +5,8 @@
 #include <polyfem/GenericProblem.hpp>
 #include <polyfem/StringUtils.hpp>
 
+#include "raster.hpp"
+
 #include <geogram/basic/command_line.h>
 #include <geogram/basic/command_line_args.h>
 
@@ -659,7 +661,157 @@ PYBIND11_MODULE(polyfempy, m)
 								   throw "Updating BC works only for Tensor GenericProblems";
 							   }
 						   },
-						   "updates pressure boundary", py::arg("id"), py::arg("val"), py::arg("interp") = std::string(""));
+						   "updates pressure boundary", py::arg("id"), py::arg("val"), py::arg("interp") = std::string(""))
+					   .def(
+						   "render", [](polyfem::State &self, int width, int height, const Eigen::Vector3d &camera_position, const double camera_fov, const double camera_fl, const bool is_perspective, const Eigen::Vector3d &left, const Eigen::Vector3d &up, const Eigen::Vector3d &ambient_light, const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>> &lights, const Eigen::Vector3d &diffuse_color, const Eigen::Vector3d &specular_color, const double specular_exponent)
+						   {
+							   using namespace renderer;
+							   using namespace Eigen;
+
+							   Material material;
+							   material.diffuse_color = diffuse_color;
+							   material.specular_color = specular_color;
+							   material.specular_exponent = specular_exponent;
+
+							   Eigen::Matrix<FrameBufferAttributes, Eigen::Dynamic, Eigen::Dynamic> frameBuffer(width, height);
+							   UniformAttributes uniform;
+
+							   const Vector3d w = -left.normalized();
+							   const Vector3d u = up.cross(w).normalized();
+							   const Vector3d v = w.cross(u);
+
+							   Matrix4d M_cam_inv;
+							   M_cam_inv << u(0), v(0), w(0), camera_position(0),
+								   u(1), v(1), w(1), camera_position(1),
+								   u(2), v(2), w(2), camera_position(2),
+								   0, 0, 0, 1;
+
+							   uniform.M_cam = M_cam_inv.inverse();
+
+							   {
+								   const double camera_ar = double(width) / height;
+								   double t = tan(camera_fov / 2) * camera_fl;
+								   double b = -t;
+								   double r = t * camera_ar;
+								   double l = -r;
+								   double n = -camera_fl;
+								   double f = 5 * n;
+
+								   uniform.M_orth << 2 / (r - l), 0, 0, -(r + l) / (r - l),
+									   0, 2 / (t - b), 0, -(t + b) / (t - b),
+									   0, 0, 2 / (n - f), -(n + f) / (n - f),
+									   0, 0, 0, 1;
+								   Matrix4d P;
+								   if (is_perspective)
+								   {
+									   P << n, 0, 0, 0,
+										   0, n, 0, 0,
+										   0, 0, n + f, -f * n,
+										   0, 0, 1, 0;
+								   }
+								   else
+								   {
+									   P << 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1;
+								   }
+
+								   uniform.M = uniform.M_orth * P * uniform.M_cam;
+							   }
+
+							   Program program;
+							   program.VertexShader = [&](const VertexAttributes &va, const UniformAttributes &uniform)
+							   {
+								   VertexAttributes out;
+								   out.position = uniform.M * va.position;
+								   Vector3d color = ambient_light;
+
+								   Vector3d hit(va.position(0), va.position(1), va.position(2));
+								   for (const auto &l : lights)
+								   {
+									   Vector3d Li = (l.first - hit).normalized();
+									   Vector3d N = va.normal;
+									   Vector3d diffuse = va.material.diffuse_color * std::max(Li.dot(N), 0.0);
+									   Vector3d H;
+									   if (is_perspective)
+									   {
+										   H = (Li + (camera_position - hit).normalized()).normalized();
+									   }
+									   else
+									   {
+										   H = (Li - left.normalized()).normalized();
+									   }
+									   const Vector3d specular = va.material.specular_color * std::pow(std::max(N.dot(H), 0.0), va.material.specular_exponent);
+									   const Vector3d D = l.first - hit;
+									   color += (diffuse + specular).cwiseProduct(l.second) / D.squaredNorm();
+								   }
+								   out.color = color;
+								   return out;
+							   };
+
+							   program.FragmentShader = [](const VertexAttributes &va, const UniformAttributes &uniform)
+							   {
+								   FragmentAttributes out(va.color(0), va.color(1), va.color(2));
+								   out.depth = -va.position(2);
+								   return out;
+							   };
+
+							   program.BlendingShader = [](const FragmentAttributes &fa, const FrameBufferAttributes &previous)
+							   {
+								   if (fa.depth < previous.depth)
+								   {
+									   FrameBufferAttributes out(fa.color[0] * 255, fa.color[1] * 255, fa.color[2] * 255, fa.color[3] * 255);
+									   out.depth = fa.depth;
+									   return out;
+								   }
+								   else
+								   {
+									   return previous;
+								   }
+							   };
+
+							   const Eigen::MatrixXd tmp = self.boundary_nodes_pos + self.sol;
+							   Eigen::MatrixXd vnormals(tmp.rows(), 3);
+							   //    Eigen::MatrixXd areas(tmp.rows(), 1);
+							   vnormals.setZero();
+							   //    areas.setZero();
+							   Eigen::MatrixXd fnormals(self.boundary_triangles.rows(), 3);
+
+							   for (int i = 0; i < self.boundary_triangles.rows(); ++i)
+							   {
+								   const Vector3d l1 = tmp.row(self.boundary_triangles(i, 1)) - tmp.row(self.boundary_triangles(i, 0));
+								   const Vector3d l2 = tmp.row(self.boundary_triangles(i, 2)) - tmp.row(self.boundary_triangles(i, 0));
+								   const auto nn = l1.cross(l2);
+								   const double area = nn.norm();
+								   fnormals.row(i) = nn / area;
+
+								   for (int j = 0; j < 3; j++)
+								   {
+									   int vid = self.boundary_triangles(i, j);
+									   vnormals.row(vid) += nn;
+									   //    areas(vid) += area;
+								   }
+							   }
+
+							   std::vector<VertexAttributes> vertices;
+							   for (int i = 0; i < self.boundary_triangles.rows(); ++i)
+							   {
+								   for (int j = 0; j < 3; j++)
+								   {
+									   int vid = self.boundary_triangles(i, j);
+									   VertexAttributes va(tmp(vid, 0), tmp(vid, 1), tmp(vid, 2));
+									   va.material = material;
+									   va.normal = vnormals.row(vid).normalized();
+									   vertices.push_back(va);
+								   }
+							   }
+
+							   rasterize_triangles(program, uniform, vertices, frameBuffer);
+
+							   std::vector<uint8_t> image;
+							   framebuffer_to_uint8(frameBuffer, image);
+
+							   return image;
+						   },
+						   "renders scene");
 
 	solver.doc() = "Polyfem solver";
 

src/raster.cpp

@@ -0,0 +1,147 @@
+#include "raster.hpp"
+#include <iostream>
+
+namespace renderer
+{
+	void rasterize_triangle(const Program &program, const UniformAttributes &uniform, const VertexAttributes &v1, const VertexAttributes &v2, const VertexAttributes &v3, Eigen::Matrix<FrameBufferAttributes, Eigen::Dynamic, Eigen::Dynamic> &frameBuffer)
+	{
+		Eigen::Matrix<double, 3, 4> p;
+		p.row(0) = v1.position.array() / v1.position[3];
+		p.row(1) = v2.position.array() / v2.position[3];
+		p.row(2) = v3.position.array() / v3.position[3];
+
+		p.col(0) = ((p.col(0).array() + 1.0) / 2.0) * frameBuffer.rows();
+		p.col(1) = ((p.col(1).array() + 1.0) / 2.0) * frameBuffer.cols();
+
+		int lx = std::floor(p.col(0).minCoeff());
+		int ly = std::floor(p.col(1).minCoeff());
+		int ux = std::ceil(p.col(0).maxCoeff());
+		int uy = std::ceil(p.col(1).maxCoeff());
+
+		lx = std::min(std::max(lx, int(0)), int(frameBuffer.rows() - 1));
+		ly = std::min(std::max(ly, int(0)), int(frameBuffer.cols() - 1));
+		ux = std::max(std::min(ux, int(frameBuffer.rows() - 1)), int(0));
+		uy = std::max(std::min(uy, int(frameBuffer.cols() - 1)), int(0));
+
+		Eigen::Matrix3d A;
+		A.col(0) = p.row(0).segment(0, 3);
+		A.col(1) = p.row(1).segment(0, 3);
+		A.col(2) = p.row(2).segment(0, 3);
+		A.row(2) << 1.0, 1.0, 1.0;
+
+		Eigen::Matrix3d Ai = A.inverse();
+
+		for (unsigned i = lx; i <= ux; i++)
+		{
+			for (unsigned j = ly; j <= uy; j++)
+			{
+
+				Eigen::Vector3d pixel(i + 0.5, j + 0.5, 1);
+				Eigen::Vector3d b = Ai * pixel;
+				if (b.minCoeff() >= 0)
+				{
+					VertexAttributes va = VertexAttributes::interpolate(v1, v2, v3, b[0], b[1], b[2]);
+
+					if (va.position[2] >= -1 && va.position[2] <= 1)
+					{
+						FragmentAttributes frag = program.FragmentShader(va, uniform);
+						frameBuffer(i, j) = program.BlendingShader(frag, frameBuffer(i, j));
+					}
+				}
+			}
+		}
+	}
+
+	void rasterize_triangles(const Program &program, const UniformAttributes &uniform, const std::vector<VertexAttributes> &vertices, Eigen::Matrix<FrameBufferAttributes, Eigen::Dynamic, Eigen::Dynamic> &frameBuffer)
+	{
+		std::vector<VertexAttributes> v(vertices.size());
+		for (unsigned i = 0; i < vertices.size(); i++)
+			v[i] = program.VertexShader(vertices[i], uniform);
+
+		for (unsigned i = 0; i < vertices.size() / 3; i++)
+			rasterize_triangle(program, uniform, v[i * 3 + 0], v[i * 3 + 1], v[i * 3 + 2], frameBuffer);
+	}
+
+	void rasterize_line(const Program &program, const UniformAttributes &uniform, const VertexAttributes &v1, const VertexAttributes &v2, double line_thickness, Eigen::Matrix<FrameBufferAttributes, Eigen::Dynamic, Eigen::Dynamic> &frameBuffer)
+	{
+		Eigen::Matrix<double, 2, 4> p;
+		p.row(0) = v1.position.array() / v1.position[3];
+		p.row(1) = v2.position.array() / v2.position[3];
+
+		p.col(0) = ((p.col(0).array() + 1.0) / 2.0) * frameBuffer.rows();
+		p.col(1) = ((p.col(1).array() + 1.0) / 2.0) * frameBuffer.cols();
+
+		int lx = std::floor(p.col(0).minCoeff() - line_thickness);
+		int ly = std::floor(p.col(1).minCoeff() - line_thickness);
+		int ux = std::ceil(p.col(0).maxCoeff() + line_thickness);
+		int uy = std::ceil(p.col(1).maxCoeff() + line_thickness);
+
+		lx = std::min(std::max(lx, int(0)), int(frameBuffer.rows() - 1));
+		ly = std::min(std::max(ly, int(0)), int(frameBuffer.cols() - 1));
+		ux = std::max(std::min(ux, int(frameBuffer.rows() - 1)), int(0));
+		uy = std::max(std::min(uy, int(frameBuffer.cols() - 1)), int(0));
+
+		Eigen::Vector2f l1(p(0, 0), p(0, 1));
+		Eigen::Vector2f l2(p(1, 0), p(1, 1));
+
+		double t = -1;
+		double ll = (l1 - l2).squaredNorm();
+
+		for (unsigned i = lx; i <= ux; i++)
+		{
+			for (unsigned j = ly; j <= uy; j++)
+			{
+
+				Eigen::Vector2f pixel(i + 0.5, j + 0.5);
+
+				if (ll == 0.0)
+					t = 0;
+				else
+				{
+					t = (pixel - l1).dot(l2 - l1) / ll;
+					t = std::fmax(0, std::fmin(1, t));
+				}
+
+				Eigen::Vector2f pixel_p = l1 + t * (l2 - l1);
+
+				if ((pixel - pixel_p).squaredNorm() < (line_thickness * line_thickness))
+				{
+					VertexAttributes va = VertexAttributes::interpolate(v1, v2, v1, 1 - t, t, 0);
+					FragmentAttributes frag = program.FragmentShader(va, uniform);
+					frameBuffer(i, j) = program.BlendingShader(frag, frameBuffer(i, j));
+				}
+			}
+		}
+	}
+
+	void rasterize_lines(const Program &program, const UniformAttributes &uniform, const std::vector<VertexAttributes> &vertices, double line_thickness, Eigen::Matrix<FrameBufferAttributes, Eigen::Dynamic, Eigen::Dynamic> &frameBuffer)
+	{
+		std::vector<VertexAttributes> v(vertices.size());
+		for (unsigned i = 0; i < vertices.size(); i++)
+			v[i] = program.VertexShader(vertices[i], uniform);
+
+		for (unsigned i = 0; i < vertices.size() / 2; i++)
+			rasterize_line(program, uniform, v[i * 2 + 0], v[i * 2 + 1], line_thickness, frameBuffer);
+	}
+
+	void framebuffer_to_uint8(const Eigen::Matrix<FrameBufferAttributes, Eigen::Dynamic, Eigen::Dynamic> &frameBuffer, std::vector<uint8_t> &image)
+	{
+		const int w = frameBuffer.rows();
+		const int h = frameBuffer.cols();
+		const int comp = 4;
+		const int stride_in_bytes = w * comp;
+		image.resize(w * h * comp, 0);
+
+		for (unsigned wi = 0; wi < w; ++wi)
+		{
+			for (unsigned hi = 0; hi < h; ++hi)
+			{
+				unsigned hif = h - 1 - hi;
+				image[(hi * w * 4) + (wi * 4) + 0] = frameBuffer(wi, hif).color[0];
+				image[(hi * w * 4) + (wi * 4) + 1] = frameBuffer(wi, hif).color[1];
+				image[(hi * w * 4) + (wi * 4) + 2] = frameBuffer(wi, hif).color[2];
+				image[(hi * w * 4) + (wi * 4) + 3] = frameBuffer(wi, hif).color[3];
+			}
+		}
+	}
+}