volmetrics: eb87d4a12bd3 src/volume.cc

volmetrics

view src/volume.cc @ 34:eb87d4a12bd3

histogram

author	Eleni Maria Stea <elene.mst@gmail.com>
date	Mon, 18 Aug 2014 10:44:11 +0300
parents	e548d95e0667
children	df4a277adb82

line source

1 #include <GL/glew.h>

3 #include <assert.h>

5 #include <ctype.h>

6 #include <errno.h>

7 #include <float.h>

8 #include <stdio.h>

9 #include <string.h>

11 #include <string>

12 #include <vector>

14 #include <metasurf.h>

16 #include "mesh.h"

17 #include "volume.h"

19 static char *strip_whitespaces(char *buf);

21 Volume::Volume()

22 {

23 width = height = 0;

24 zaspect = 1;

25 x_rot = 0;

26 vol_tex_valid = false;

28 memset(histogram, 0, HIST_SIZE * sizeof(int));

29 max_histogram_value = 0;

31 glGenTextures(1, &vol_tex);

32 }

34 Volume::~Volume()

35 {

36 glDeleteTextures(1, &vol_tex);

37 }

39 void Volume::create_vol_tex() const

40 {

41 if(slices.empty())

42 return;

44 glBindTexture(GL_TEXTURE_3D, vol_tex);

45 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

46 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);

47 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

48 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

49 glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE);

51 assert(glGetError() == GL_NO_ERROR);

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

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

54 assert(glGetError() == GL_NO_ERROR);

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

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

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

59 }

60 }

62 void Volume::create_histogram()

63 {

64 if (slices.empty())

65 return;

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

68 {

69 float *slice_pixels = (float*)slices[i].get_pixels();

70 for (int j=0; j < width * height; j++)

71 {

72 float val = slice_pixels[j];

74 int ival = (int)(val * 256);

75 histogram[ival]++;

76 }

77 }

79 max_histogram_value = 0;

80 for (int i=0; i<HIST_SIZE; i++) {

81 if (i < 10)

82 continue;

83 max_histogram_value = std::max(max_histogram_value, histogram[i]);

84 }

85 }

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

88 {

89 char up = 'y';

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

92 if(!fp) {

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

94 return false;

95 }

97 char buf[512];

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

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

100 return false;

101 }

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

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

104 return false;

105 }

106

107 bool reading_slices = false;

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

109

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

111 char *line = strip_whitespaces(buf);

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

113 continue;

114

115 if(reading_slices == false) {

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

117 reading_slices = true;

118

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

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

121 }

122 else {

123 img_fnames.push_back(line);

124 }

125 }

126

127 fclose(fp);

128

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

130 Image img;

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

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

133 }

134 else {

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

136 }

137 }

138

139 if (up == 'z')

140 x_rot = 90;

141

142 create_histogram();

143 return true;

144 }

145

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

147 {

148 if(slices.empty()) {

149 width = slice.get_width();

150 height = slice.get_height();

151 }

152 else {

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

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

155 return false;

156 }

157 }

158

159 slices.push_back(slice);

160 return true;

161 }

162

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

164 {

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

166 return 0;

167

168 return &slices[num_slice];

169 }

170

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

172 {

173 int idx = get_slice_idx_by_z(z);

174 return get_slice(idx);

175 }

176

177 int Volume::get_slice_count() const

178 {

179 return (int)slices.size();

180 }

181

182 int Volume::get_slice_idx_by_z(float z) const

183 {

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

185 int idx = z * last_slice_idx;

186 if(idx < 0) {

187 idx = 0;

188 }

189 else if(idx > last_slice_idx) {

190 idx = last_slice_idx;

191 }

192

193 return idx;

194 }

195

196 unsigned int Volume::get_texture() const

197 {

198 if(!vol_tex_valid) {

199 create_vol_tex();

200 vol_tex_valid = true;

201 }

202

203 return vol_tex;

204 }

205

206 float Volume::get_volume_rotation() const

207 {

208 return x_rot;

209 }

210

211 static Mesh *cur_mesh;

212 static Volume *cur_vol;

213 static float low_thres, high_thres;

214

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

216 {

217 if(x < low) return 0.0;

218 if(x >= high) return 1.0;

219

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

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

222 }

223

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

225 {

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

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

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

229 }

230

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

232 {

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

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

235

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

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

238

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

240 return transfer_function(val, low_thres, high_thres);

241 }

242

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

244 {

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

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

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

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

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

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

251

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

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

254 }

255

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

257 {

258 if (tmin > tmax) {

259 float t = tmax;

260 tmax = tmin;

261 tmin = t;

262 }

263 low_thres = tmin; high_thres = tmax;

264

265 if(xres <= 0)

266 xres = width;

267 if(yres <= 0)

268 yres = height;

269 if(zres <= 0)

270 zres = slices.size();

271

272 cur_mesh = mesh;

273 cur_vol = this;

274

275 metasurface *ms = msurf_create();

276 msurf_threshold(ms, 0.5);

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

278 msurf_vertex_func(ms, cb_vertex);

279 msurf_eval_func(ms, cb_eval);

280 msurf_polygonize(ms);

281 msurf_free(ms);

282 }

283

284 void Volume::draw() const

285 {

286 }

287

288 static char *strip_whitespaces(char *buf)

289 {

290 while(*buf && isspace(*buf))

291 buf++;

292

293 if(!*buf)

294 return 0;

295

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

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

298 end--;

299 end[1] = 0;

300 return buf;

301 }