1.1 --- a/src/volume.cc	Wed Dec 31 07:57:42 2014 +0200
     1.2 +++ b/src/volume.cc	Thu Jan 01 06:36:45 2015 +0200
     1.3 @@ -94,6 +94,10 @@
     1.4  VoxelVolume::VoxelVolume()
     1.5  {
     1.6  	size[0] = size[1] = size[2] = 0;
     1.7 +
     1.8 +	hist_valid = false;
     1.9 +	hist = 0;
    1.10 +	num_hist_samples = 0;
    1.11  }
    1.12  
    1.13  VoxelVolume::~VoxelVolume()
    1.14 @@ -151,6 +155,7 @@
    1.15  	}
    1.16  
    1.17  	size[2] = slices.size();
    1.18 +	hist_valid = false;
    1.19  	fclose(fp);
    1.20  	return true;
    1.21  }
    1.22 @@ -223,6 +228,39 @@
    1.23  	return slices[z].pixels[y * size[0] + x];
    1.24  }
    1.25  
    1.26 +float *VoxelVolume::calc_histogram(int hist_samples)
    1.27 +{
    1.28 +	if(hist && hist_samples == num_hist_samples && hist_valid) {
    1.29 +		return hist;
    1.30 +	}
    1.31 +
    1.32 +	int num_pixels = size[0] * size[1] * size[2];
    1.33 +	if(!num_pixels) {
    1.34 +		return 0;
    1.35 +	}
    1.36 +
    1.37 +	delete [] hist;
    1.38 +	hist = new float[hist_samples];
    1.39 +	memset(hist, 0, hist_samples * sizeof *hist);
    1.40 +
    1.41 +	for(int i=0; i<size[2]; i++) {
    1.42 +		float *pptr = slices[i].pixels;
    1.43 +		for(int j=0; j<size[0] * size[1]; j++) {
    1.44 +			int idx = (int)(*pptr++ * (float)hist_samples);
    1.45 +
    1.46 +			hist[idx] += 1.0;
    1.47 +		}
    1.48 +	}
    1.49 +
    1.50 +	for(int i=0; i<hist_samples; i++) {
    1.51 +		hist[i] /= (float)num_pixels;
    1.52 +	}
    1.53 +
    1.54 +	hist_valid = true;
    1.55 +	num_hist_samples = hist_samples;
    1.56 +	return hist;
    1.57 +}
    1.58 +
    1.59  static char *strip_space(char *s)
    1.60  {
    1.61  	while(*s && isspace(*s)) s++;