nuclear@1: #include nuclear@1: #include nuclear@1: #include nuclear@2: #include nuclear@1: #include nuclear@1: #include nuclear@1: #include "volume.h" nuclear@1: nuclear@1: static char *strip_space(char *s); nuclear@1: nuclear@1: nuclear@1: Volume::~Volume() nuclear@1: { nuclear@1: } nuclear@1: nuclear@1: int Volume::num_samples(int dim) const nuclear@1: { nuclear@1: return 0; nuclear@1: } nuclear@1: nuclear@1: void Volume::normalf(float *norm, float x, float y, float z, float delta) nuclear@1: { nuclear@1: float dx, dy, dz; nuclear@1: dx = dy = dz = delta; nuclear@1: nuclear@1: if(num_samples(0) > 0) { nuclear@1: // discrete volume nuclear@1: dx /= (float)num_samples(0); nuclear@1: dy /= (float)num_samples(1); nuclear@1: dz /= (float)num_samples(2); nuclear@1: } nuclear@1: nuclear@1: norm[0] = valuef(x + dx, y, z) - valuef(x - dx, y, z); nuclear@1: norm[1] = valuef(x, y + dy, z) - valuef(x, y - dy, z); nuclear@1: norm[2] = valuef(x, y, z + dz) - valuef(x, y, z - dz); nuclear@1: nuclear@1: float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]); nuclear@1: if(len != 0.0) { nuclear@1: norm[0] /= len; nuclear@1: norm[1] /= len; nuclear@1: norm[2] /= len; nuclear@1: } nuclear@1: } nuclear@1: nuclear@1: void Volume::normali(float *norm, int x, int y, int z) nuclear@1: { nuclear@1: int sz[3]; nuclear@1: nuclear@1: if((sz[0] = num_samples(0)) <= 0) { nuclear@1: // if it's a continuous volume, just call normalf nuclear@1: normalf(norm, x, y, z); nuclear@1: return; nuclear@1: } nuclear@1: sz[1] = num_samples(1); nuclear@1: sz[2] = num_samples(2); nuclear@1: nuclear@1: int prevx = x <= 0 ? 0 : x; nuclear@1: int nextx = x >= sz[0] - 1 ? sz[0] - 1 : x; nuclear@1: int prevy = y <= 0 ? 0 : y; nuclear@1: int nexty = y >= sz[1] - 1 ? sz[1] - 1 : y; nuclear@1: int prevz = z <= 0 ? 0 : z; nuclear@1: int nextz = z >= sz[2] - 1 ? sz[2] - 1 : z; nuclear@1: nuclear@1: norm[0] = valuei(nextx, y, z) - valuei(prevx, y, z); nuclear@1: norm[1] = valuei(x, nexty, z) - valuei(x, prevy, z); nuclear@1: norm[2] = valuei(x, y, nextz) - valuei(x, y, prevz); nuclear@1: nuclear@1: float len = sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]); nuclear@1: if(len != 0.0) { nuclear@1: norm[0] /= len; nuclear@1: norm[1] /= len; nuclear@1: norm[2] /= len; nuclear@1: } nuclear@1: } nuclear@1: nuclear@1: // ---- VoxelVolume (discrete) implementation ---- nuclear@1: VoxelVolume::VoxelVolume() nuclear@1: { nuclear@1: size[0] = size[1] = size[2] = 0; nuclear@1: } nuclear@1: nuclear@1: bool VoxelVolume::load(const char *fname) nuclear@1: { nuclear@2: char *prefix = (char*)alloca(strlen(fname) + 1); nuclear@2: strcpy(prefix, fname); nuclear@2: char *slash = strrchr(prefix, '/'); nuclear@2: if(slash) { nuclear@2: *slash = 0; nuclear@2: } else { nuclear@2: prefix = 0; nuclear@2: } nuclear@2: nuclear@1: FILE *fp = fopen(fname, "r"); nuclear@1: if(!fp) { nuclear@1: fprintf(stderr, "failed to open file: %s: %s\n", fname, strerror(errno)); nuclear@1: return false; nuclear@1: } nuclear@1: nuclear@2: char buf[256], path[300]; nuclear@1: while(fgets(buf, sizeof buf, fp)) { nuclear@1: char *line = strip_space(buf); nuclear@2: sprintf(path, "%s/%s", prefix, line); nuclear@1: nuclear@1: Image img; nuclear@2: if(!img.load(path)) { nuclear@1: slices.clear(); nuclear@1: return false; nuclear@1: } nuclear@1: nuclear@1: if(slices.empty()) { nuclear@1: size[0] = img.width; nuclear@1: size[1] = img.height; nuclear@1: } else { nuclear@1: if(img.width != size[0] || img.height != size[1]) { nuclear@1: fprintf(stderr, "slice %d \"%s\" size mismatch (%dx%d, %dx%d expected)\n", nuclear@1: (int)slices.size(), line, img.width, img.height, size[0], size[1]); nuclear@1: slices.clear(); nuclear@1: return false; nuclear@1: } nuclear@1: } nuclear@1: nuclear@1: slices.push_back(img); nuclear@1: img.pixels = 0; // otherwise the destructor will free it nuclear@1: } nuclear@1: nuclear@1: size[2] = slices.size(); nuclear@1: fclose(fp); nuclear@1: return true; nuclear@1: } nuclear@1: nuclear@1: int VoxelVolume::num_samples(int dim) const nuclear@1: { nuclear@1: return size[dim]; nuclear@1: } nuclear@1: nuclear@1: static inline int clamp(int x, int low, int high) nuclear@1: { nuclear@1: return x < low ? low : (x > high ? high : x); nuclear@1: } nuclear@1: nuclear@1: static inline float lookup(int x0, int y0, int x1, int y1, float tx, float ty, nuclear@1: float *pixels, int width, int height) nuclear@1: { nuclear@1: float v00, v01, v10, v11; nuclear@1: nuclear@1: v00 = pixels[y0 * width + x0]; nuclear@1: v01 = pixels[y1 * width + x0]; nuclear@1: v10 = pixels[y0 * width + x1]; nuclear@1: v11 = pixels[y1 * width + x1]; nuclear@1: nuclear@1: float v0 = v00 + (v01 - v00) * ty; nuclear@1: float v1 = v10 + (v11 - v10) * ty; nuclear@1: nuclear@1: return v0 + (v1 - v0) * tx; nuclear@1: } nuclear@1: nuclear@1: float VoxelVolume::valuef(float x, float y, float z) const nuclear@1: { nuclear@1: if(slices.empty()) { nuclear@1: return 0.0f; nuclear@1: } nuclear@1: nuclear@1: float floor_x = floor(x); nuclear@1: float ceil_x = ceil(x); nuclear@1: float tx = x - floor_x; nuclear@1: nuclear@1: float floor_y = floor(y); nuclear@1: float ceil_y = ceil(y); nuclear@1: float ty = y - floor_y; nuclear@1: nuclear@1: float floor_z = floor(z); nuclear@1: float ceil_z = ceil(z); nuclear@1: float tz = z - floor_z; nuclear@1: nuclear@1: int x0 = clamp(floor_x, 0, size[0] - 1); nuclear@1: int x1 = clamp(ceil_x, 0, size[0] - 1); nuclear@1: nuclear@1: int y0 = clamp(floor_y, 0, size[1] - 1); nuclear@1: int y1 = clamp(ceil_y, 0, size[1] - 1); nuclear@1: nuclear@1: int s0 = clamp(floor_z, 0, size[2] - 1); nuclear@1: int s1 = clamp(ceil_z, 0, size[2] - 1); nuclear@1: nuclear@1: float val_s0 = lookup(x0, y0, x1, y1, tx, ty, slices[s0].pixels, size[0], size[1]); nuclear@1: float val_s1 = lookup(x0, y0, x1, y1, tx, ty, slices[s1].pixels, size[0], size[1]); nuclear@1: nuclear@1: return val_s0 + (val_s1 - val_s0) * tz; nuclear@1: } nuclear@1: nuclear@1: float VoxelVolume::valuei(int x, int y, int z) const nuclear@1: { nuclear@1: x = clamp(x, 0, size[0] - 1); nuclear@1: y = clamp(y, 0, size[1] - 1); nuclear@1: z = clamp(z, 0, size[2] - 1); nuclear@1: nuclear@1: return slices[z].pixels[y * size[0] + x]; nuclear@1: } nuclear@1: nuclear@1: static char *strip_space(char *s) nuclear@1: { nuclear@1: while(*s && isspace(*s)) s++; nuclear@1: if(!*s) return 0; nuclear@1: nuclear@1: char *end = s + strlen(s) - 1; nuclear@1: while(end > s && isspace(*end)) end--; nuclear@1: end[1] = 0; nuclear@1: nuclear@1: return s; nuclear@1: }