volmetrics: cfeecd093730 src/volume.cc

volmetrics

view src/volume.cc @ 31:cfeecd093730

minor change

author	Eleni Maria Stea <elene.mst@gmail.com>
date	Thu, 01 May 2014 16:28:12 +0300
parents	4e120dcd55ec
children	eb87d4a12bd3

line source

1 #include <GL/glew.h>

3 #include <assert.h>

5 #include <ctype.h>

6 #include <errno.h>

7 #include <stdio.h>

8 #include <string.h>

10 #include <string>

11 #include <vector>

13 #include <metasurf.h>

15 #include "mesh.h"

16 #include "volume.h"

18 static char *strip_whitespaces(char *buf);

20 Volume::Volume()

21 {

22 width = height = 0;

23 zaspect = 1;

24 x_rot = 0;

25 vol_tex_valid = false;

27 glGenTextures(1, &vol_tex);

28 }

30 Volume::~Volume()

31 {

32 glDeleteTextures(1, &vol_tex);

33 }

35 void Volume::create_vol_tex() const

36 {

37 if(slices.empty())

38 return;

40 glBindTexture(GL_TEXTURE_3D, vol_tex);

41 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

42 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

43 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

44 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

45 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);

47 assert(glGetError() == GL_NO_ERROR);

48 glTexImage3D(GL_TEXTURE_3D, 0, GL_LUMINANCE, width, height,

49 slices.size(), 0, GL_LUMINANCE, GL_FLOAT, 0);

50 assert(glGetError() == GL_NO_ERROR);

52 for(size_t i=0; i<slices.size(); i++) {

53 glTexSubImage3D(GL_TEXTURE_3D, 0, 0, 0, i, width, height,

54 1, GL_LUMINANCE, GL_FLOAT, slices[i].get_pixels());

55 }

56 }

58 bool Volume::load(const char *fname)

59 {

60 char up = 'y';

62 FILE *fp = fopen(fname, "r");

63 if(!fp) {

64 fprintf(stderr, "Failed to open file: %s: %s\n", fname, strerror(errno));

65 return false;

66 }

68 char buf[512];

69 if(!fgets(buf, sizeof buf, fp)) {

70 fprintf(stderr, "Empty file: %s.\n", fname);

71 return false;

72 }

73 if(strstr(buf, "VOLUME") != buf) {

74 fprintf(stderr, "Invalid volume file format: %s\n", fname);

75 return false;

76 }

78 bool reading_slices = false;

79 std::vector<std::string> img_fnames;

81 while(fgets(buf, sizeof buf, fp)) {

82 char *line = strip_whitespaces(buf);

83 if(!line || *line == 0 || *line == '#')

84 continue;

86 if(reading_slices == false) {

87 if(strcmp(line, "SLICES") == 0)

88 reading_slices = true;

90 sscanf(line, "zaspect = %f", &zaspect);

91 sscanf(line, "up = %c", &up);

92 }

93 else {

94 img_fnames.push_back(line);

95 }

96 }

98 fclose(fp);

100 for(size_t i=0; i<img_fnames.size(); i++) {

101 Image img;

102 if(img.load(img_fnames[i].c_str())) {

103 push_slice(std::move(img));

104 }

105 else {

106 fprintf(stderr, "Failed to load slice: %s\n", img_fnames[i].c_str());

107 }

108 }

109

110 if (up == 'z')

111 x_rot = 90;

112

113 return true;

114 }

115

116 bool Volume::push_slice(Image &&slice)

117 {

118 if(slices.empty()) {

119 width = slice.get_width();

120 height = slice.get_height();

121 }

122 else {

123 if(width != slice.get_width() || height != slice.get_height()) {

124 fprintf(stderr, "failed to load slice no: %d\n", (int)slices.size() + 1);

125 return false;

126 }

127 }

128

129 slices.push_back(slice);

130 return true;

131 }

132

133 const Image *Volume::get_slice(int num_slice) const

134 {

135 if(num_slice >= (int)slices.size() || num_slice < 0)

136 return 0;

137

138 return &slices[num_slice];

139 }

140

141 const Image *Volume::get_slice_by_z(float z) const

142 {

143 int idx = get_slice_idx_by_z(z);

144 return get_slice(idx);

145 }

146

147 int Volume::get_slice_count() const

148 {

149 return (int)slices.size();

150 }

151

152 int Volume::get_slice_idx_by_z(float z) const

153 {

154 int last_slice_idx = (int)slices.size() - 1;

155 int idx = z * last_slice_idx;

156 if(idx < 0) {

157 idx = 0;

158 }

159 else if(idx > last_slice_idx) {

160 idx = last_slice_idx;

161 }

162

163 return idx;

164 }

165

166 unsigned int Volume::get_texture() const

167 {

168 if(!vol_tex_valid) {

169 create_vol_tex();

170 vol_tex_valid = true;

171 }

172

173 return vol_tex;

174 }

175

176 float Volume::get_volume_rotation() const

177 {

178 return x_rot;

179 }

180

181 static Mesh *cur_mesh;

182 static Volume *cur_vol;

183 static float low_thres, high_thres;

184

185 static float smoothstep(float x, float low, float high)

186 {

187 if(x < low) return 0.0;

188 if(x >= high) return 1.0;

189

190 x = (x - low) / (high - low);

191 return x * x * (3.0 - 2.0 * x);

192 }

193

194 float transfer_function(float x, float low_thres, float high_thres)

195 {

196 float dt = 0.25 * (high_thres - low_thres);

197 return smoothstep(x, low_thres - dt, low_thres + dt) *

198 (1 - smoothstep(x, high_thres - dt, high_thres + dt));

199 }

200

201 static float cb_eval(float x, float y, float z)

202 {

203 const Image *img = cur_vol->get_slice_by_z((z + 1) / 2.0);

204 const float *pixels = img->get_pixels();

205

206 int px = (x + 1) / 2.0 * img->get_width();

207 int py = (y + 1) / 2.0 * img->get_height();

208

209 float val = pixels[px + img->get_width() * py];

210 return transfer_function(val, low_thres, high_thres);

211 }

212

213 static void cb_vertex(float x, float y, float z)

214 {

215 float dx = 1.0 / cur_vol->get_slice(0)->get_width();

216 float dy = 1.0 / cur_vol->get_slice(0)->get_height();

217 float dz = 1.0 / cur_vol->get_slice_count();

218 float dfdx = cb_eval(x - dx, y, z) - cb_eval(x + dx, y, z);

219 float dfdy = cb_eval(x, y - dy, z) - cb_eval(x, y + dy, z);

220 float dfdz = cb_eval(x, y, z - dz) - cb_eval(x, y, z + dz);

221

222 cur_mesh->add_normal(Vector3(dfdx, dfdy, dfdz));

223 cur_mesh->add_vertex(Vector3(x, y, z));

224 }

225

226 void Volume::create_mesh(Mesh *mesh, float tmin, float tmax, int xres, int yres, int zres)

227 {

228 if (tmin > tmax) {

229 float t = tmax;

230 tmax = tmin;

231 tmin = t;

232 }

233 low_thres = tmin; high_thres = tmax;

234

235 if(xres <= 0)

236 xres = width;

237 if(yres <= 0)

238 yres = height;

239 if(zres <= 0)

240 zres = slices.size();

241

242 cur_mesh = mesh;

243 cur_vol = this;

244

245 metasurface *ms = msurf_create();

246 msurf_threshold(ms, 0.5);

247 msurf_resolution(ms, xres, yres, zres);

248 msurf_vertex_func(ms, cb_vertex);

249 msurf_eval_func(ms, cb_eval);

250 msurf_polygonize(ms);

251 msurf_free(ms);

252 }

253

254 void Volume::draw() const

255 {

256 }

257

258 static char *strip_whitespaces(char *buf)

259 {

260 while(*buf && isspace(*buf))

261 buf++;

262

263 if(!*buf)

264 return 0;

265

266 char *end = buf + strlen(buf) - 1;

267 while(end > buf && isspace(*end))

268 end--;

269 end[1] = 0;

270 return buf;

271 }