elene@6: #include <GL/glew.h>
elene@6: 
elene@6: #include <assert.h>
elene@6: 
elene@5: #include <ctype.h>
elene@5: #include <errno.h>
elene@4: #include <stdio.h>
elene@5: #include <string.h>
elene@5: 
elene@5: #include <string>
elene@5: #include <vector>
elene@4: 
elene@9: #include <metasurf.h>
elene@9: 
elene@9: #include "mesh.h"
elene@3: #include "volume.h"
elene@4: 
elene@5: static char *strip_whitespaces(char *buf);
elene@5: 
elene@4: Volume::Volume()
elene@4: {
elene@4: 	width = height = 0;
elene@5: 	zaspect = 1;
elene@30: 	x_rot = 0;
elene@6: 	vol_tex_valid = false;
elene@6: 
elene@6: 	glGenTextures(1, &vol_tex);
elene@4: }
elene@4: 
elene@6: Volume::~Volume()
elene@6: {
elene@6: 	glDeleteTextures(1, &vol_tex);
elene@6: }
elene@6: 
elene@6: void Volume::create_vol_tex() const
elene@6: {
elene@6: 	if(slices.empty())
elene@6: 		return;
elene@6: 
elene@6: 	glBindTexture(GL_TEXTURE_3D, vol_tex);
elene@6: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
elene@6: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
elene@6: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
elene@6: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
elene@6: 	glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);
elene@6: 
elene@6: 	assert(glGetError() == GL_NO_ERROR);
elene@6: 	glTexImage3D(GL_TEXTURE_3D, 0, GL_LUMINANCE, width, height,
elene@6: 			slices.size(), 0, GL_LUMINANCE, GL_FLOAT, 0);
elene@6: 	assert(glGetError() == GL_NO_ERROR);
elene@6: 
elene@6: 	for(size_t i=0; i<slices.size(); i++) {
elene@6: 		glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, width, height,
elene@6: 				1, GL_LUMINANCE, GL_FLOAT, slices[i].get_pixels());
elene@6: 	}
elene@6: }
elene@6: 
elene@6: bool Volume::load(const char *fname)
elene@4: {
elene@30: 	char up = 'y';
elene@30: 
elene@5: 	FILE *fp = fopen(fname, "r");
elene@5: 	if(!fp) {
elene@5: 		fprintf(stderr, "Failed to open file: %s: %s\n", fname, strerror(errno));
elene@5: 		return false;
elene@5: 	}
elene@5: 
elene@5: 	char buf[512];
elene@5: 	if(!fgets(buf, sizeof buf, fp)) {
elene@5: 		fprintf(stderr, "Empty file: %s.\n", fname);
elene@5: 		return false;
elene@5: 	}
elene@5: 	if(strstr(buf, "VOLUME") != buf) {
elene@5: 		fprintf(stderr, "Invalid volume file format: %s\n", fname);
elene@5: 		return false;
elene@5: 	}
elene@5: 
elene@5: 	bool reading_slices = false;
elene@5: 	std::vector<std::string> img_fnames;
elene@5: 
elene@5: 	while(fgets(buf, sizeof buf, fp)) {
elene@5: 		char *line = strip_whitespaces(buf);
elene@5: 		if(!line || *line == 0 || *line == '#')
elene@5: 			continue;
elene@5: 
elene@5: 		if(reading_slices == false) {
elene@5: 			if(strcmp(line, "SLICES") == 0)
elene@5: 				reading_slices = true;
elene@5: 
elene@5: 			sscanf(line, "zaspect = %f", &zaspect);
elene@30: 			sscanf(line, "up = %c", &up);
elene@5: 		}
elene@5: 		else {
elene@5: 			img_fnames.push_back(line);
elene@5: 		}
elene@5: 	}
elene@5: 
elene@5: 	fclose(fp);
elene@5: 
elene@5: 	for(size_t i=0; i<img_fnames.size(); i++) {
elene@5: 		Image img;
elene@5: 		if(img.load(img_fnames[i].c_str())) {
elene@5: 			push_slice(std::move(img));
elene@5: 		}
elene@5: 		else {
elene@5: 			fprintf(stderr, "Failed to load slice: %s\n", img_fnames[i].c_str());
elene@5: 		}
elene@5: 	}
elene@5: 
elene@30: 	if (up == 'z')
elene@30: 		x_rot = 90;
elene@30: 
elene@4: 	return true;
elene@4: }
elene@4: 
elene@4: bool Volume::push_slice(Image &&slice)
elene@4: {
elene@4: 	if(slices.empty()) {
elene@4: 		width = slice.get_width();
elene@4: 		height = slice.get_height();
elene@4: 	}
elene@4: 	else {
elene@4: 		if(width != slice.get_width() || height != slice.get_height()) {
elene@4: 			fprintf(stderr, "failed to load slice no: %d\n", (int)slices.size() + 1);
elene@4: 			return false;
elene@4: 		}
elene@4: 	}
elene@4: 
elene@4: 	slices.push_back(slice);
elene@4: 	return true;
elene@4: }
elene@5: 
elene@6: const Image *Volume::get_slice(int num_slice) const
elene@6: {
elene@6: 	if(num_slice >= (int)slices.size() || num_slice < 0)
elene@6: 		return 0;
elene@6: 
elene@6: 	return &slices[num_slice];
elene@6: }
elene@6: 
elene@6: const Image *Volume::get_slice_by_z(float z) const
elene@6: {
elene@6: 	int idx = get_slice_idx_by_z(z);
elene@6: 	return get_slice(idx);
elene@6: }
elene@6: 
elene@6: int Volume::get_slice_count() const
elene@6: {
elene@6: 	return (int)slices.size();
elene@6: }
elene@6: 
elene@6: int Volume::get_slice_idx_by_z(float z) const
elene@6: {
elene@6: 	int last_slice_idx = (int)slices.size() - 1;
elene@6: 	int idx = z * last_slice_idx;
elene@6: 	if(idx < 0) {
elene@6: 		idx = 0;
elene@6: 	}
elene@6: 	else if(idx > last_slice_idx) {
elene@6: 		idx = last_slice_idx;
elene@6: 	}
elene@6: 
elene@6: 	return idx;
elene@6: }
elene@6: 
elene@6: unsigned int Volume::get_texture() const
elene@6: {
elene@6: 	if(!vol_tex_valid) {
elene@6: 		create_vol_tex();
elene@6: 		vol_tex_valid = true;
elene@6: 	}
elene@6: 
elene@6: 	return vol_tex;
elene@6: }
elene@6: 
elene@30: float Volume::get_volume_rotation() const
elene@30: {
elene@30: 	return x_rot;
elene@30: }
elene@30: 
elene@9: static Mesh *cur_mesh;
elene@9: static Volume *cur_vol;
elene@12: static float low_thres, high_thres;
elene@12: 
elene@12: static float smoothstep(float x, float low, float high)
elene@12: {
elene@12: 	if(x < low) return 0.0;
elene@12: 	if(x >= high) return 1.0;
elene@12: 
elene@12: 	x = (x - low) / (high - low);
elene@12: 	return x * x * (3.0 - 2.0 * x);
elene@12: }
elene@16: 
elene@16: float transfer_function(float x, float low_thres, float high_thres)
elene@16: {
elene@16: 	float dt = 0.25 * (high_thres - low_thres);
elene@16: 	return smoothstep(x, low_thres - dt, low_thres + dt) *
elene@16: 		(1 - smoothstep(x, high_thres - dt, high_thres + dt));
elene@16: }
elene@16: 
elene@9: static float cb_eval(float x, float y, float z)
elene@9: {
elene@9: 	const Image *img = cur_vol->get_slice_by_z((z + 1) / 2.0);
elene@9: 	const float *pixels = img->get_pixels();
elene@9: 
elene@9: 	int px = (x + 1) / 2.0 * img->get_width();
elene@9: 	int py = (y + 1) / 2.0 * img->get_height();
elene@9: 
elene@12: 	float val = pixels[px + img->get_width() * py];
elene@16: 	return transfer_function(val, low_thres, high_thres);
elene@9: }
elene@9: 
elene@11: static void cb_vertex(float x, float y, float z)
elene@11: {
elene@11: 	float dx = 1.0 / cur_vol->get_slice(0)->get_width();
elene@11: 	float dy = 1.0 / cur_vol->get_slice(0)->get_height();
elene@11: 	float dz = 1.0 / cur_vol->get_slice_count();
elene@11: 	float dfdx = cb_eval(x - dx, y, z) - cb_eval(x + dx, y, z);
elene@11: 	float dfdy = cb_eval(x, y - dy, z) - cb_eval(x, y + dy, z);
elene@11: 	float dfdz = cb_eval(x, y, z - dz) - cb_eval(x, y, z + dz);
elene@11: 
elene@11: 	cur_mesh->add_normal(Vector3(dfdx, dfdy, dfdz));
elene@11: 	cur_mesh->add_vertex(Vector3(x, y, z));
elene@11: }
elene@11: 
elene@12: void Volume::create_mesh(Mesh *mesh, float tmin, float tmax, int xres, int yres, int zres)
elene@9: {
elene@12: 	if (tmin > tmax) {
elene@12: 		float t = tmax;
elene@12: 		tmax = tmin;
elene@12: 		tmin = t;
elene@12: 	}
elene@12: 	low_thres = tmin; high_thres = tmax;
elene@12: 
elene@12: 	if(xres <= 0)
elene@12: 		xres = width;
elene@12: 	if(yres <= 0)
elene@12: 		yres = height;
elene@12: 	if(zres <= 0)
elene@12: 		zres = slices.size();
elene@12: 
elene@9: 	cur_mesh = mesh;
elene@9: 	cur_vol = this;
elene@9: 
elene@9: 	metasurface *ms = msurf_create();
elene@12: 	msurf_threshold(ms, 0.5);
elene@12: 	msurf_resolution(ms, xres, yres, zres);
elene@11: 	msurf_vertex_func(ms, cb_vertex);
elene@9: 	msurf_eval_func(ms, cb_eval);
elene@9: 	msurf_polygonize(ms);
elene@9: 	msurf_free(ms);
elene@9: }
elene@9: 
elene@8: void Volume::draw() const
elene@5: {
elene@5: }
elene@5: 
elene@5: static char *strip_whitespaces(char *buf)
elene@5: {
elene@5: 	while(*buf && isspace(*buf))
elene@5: 		buf++;
elene@5: 
elene@5: 	if(!*buf)
elene@5: 		return 0;
elene@5: 
elene@5: 	char *end = buf + strlen(buf) - 1;
elene@5: 	while(end > buf && isspace(*end))
elene@5: 		end--;
elene@5: 	end[1] = 0;
elene@5: 	return buf;
elene@5: }