- Seperate export and run into different function in the Blender addon

- Seed the Random Number Generator based on Time rather then fixed constant
- Calculate T60 from Impulse Response, Sabine and Norris-Eyring
- Calculate mesh statistics for use in Sabine and Norris-Eyring formulas
- Increase the maximum number of ray bounces

blender/render_EAR/__init__.py

@@ -348,7 +348,9 @@ def is_animated(ob):
                 num_recorders += 1
                 recorder_names.append(ob.name)
     fs.close()
-
+
+def export_and_run(filename):
+    export(filename)
     # Now start rendering the file
     if sys.platform.startswith('win'):
         os.system('start /LOW EAR render "%s"'%filename)
@@ -364,7 +366,7 @@ class EAR_OP_Export(bpy.types.Operator):
     filter_glob = StringProperty(default="*.ear",options={'HIDDEN'})
 
     def execute(self, context): 
-        export(self.filepath)
+        export_and_run(self.filepath)
         return {'FINISHED'}
 
     def invoke(self, context, event):

src/EAR.cpp

@@ -19,7 +19,7 @@
 
 #define BANNER "  ______                      _____  \n |  ____|         /\\         |  __ \\ \n | |__           /  \\        | |__) | \n |  __|         / /\\ \\       |  _  / \n | |____       / ____ \\      | | \\ \\ \n |______| (_) /_/    \\_\\ (_) |_|  \\_\\"
 #define PRODUCT "Evaluation of Acoustics using Ray-tracing"
-#define VERSION "0.1.3b"
+#define VERSION "0.1.4b"
 #define HEADER BANNER "\n " PRODUCT "\n version " VERSION
 
 #include <iostream>
@@ -52,10 +52,10 @@
 
 using boost::thread;
 
-int Render(std::string filename, int id) {
+int Render(std::string filename, float* calc_T60=0, float* T60_Sabine=0, float* T60_Eyring=0) {
 
 	// Init RNG, scene, and file input
-	gmtl::Math::seedRandom(id + 1000);
+	gmtl::Math::seedRandom((unsigned int) time(0));
 	Scene* scene = new Scene();
 	bool valid_file = Datatype::SetInput(filename);
 	if ( ! valid_file ) {
@@ -128,6 +128,9 @@ int Render(std::string filename, int id) {
 		int sf_id = 0;
 		for ( std::vector<AbstractSoundFile*>::const_iterator it = scene->sources.begin(); it != scene->sources.end(); ++ it ) {
 			for ( int band_id = 0; band_id < 3; ++ band_id ) {
+				// If we are only here to calculate the T60 reverberation time
+				// we are only going to render the mid frequency range.
+				if ( calc_T60 && band_id != 1 ) continue;
 				SoundFile* band = (*it)->Band(band_id);
 				WaveFile w;
 				w.FromFloat(band->data,band->sample_length);
@@ -136,6 +139,9 @@ int Render(std::string filename, int id) {
 				w.Save(ss.str().c_str());
 				delete band;
 			}
+			// If we are only here to calculate the T60 reverberation time
+			// we are only going to render the first sound file encountered.
+			if ( calc_T60 ) break;
 			sf_id ++;
 		}
 	}
@@ -153,13 +159,24 @@ int Render(std::string filename, int id) {
 	std::vector<SceneContext> scs;
 	for( unsigned int sound_id = 0; sound_id < scene->sources.size(); sound_id ++ ) {
 		AbstractSoundFile* sf = scene->sources[sound_id];
+		// This for loop iterates over all keyframes. If the scene contains a static
+		// configuration and no keyframes are present, keyframe_id is assigned -1.
 		for( int keyframe_id = keys?0:-1; keyframe_id < (int)(keys?keys->keys.size():0); keyframe_id ++ ) {
 			for( int band_id = 0; band_id < 3; band_id ++ ) {
+				// If we are only here to calculate the T60 reverberation time
+				// we are only going to render the mid frequency range.
+				if ( calc_T60 && band_id != 1 ) continue;
 				const float absorption_factor = 1.0f-absorption[band_id];
 				SceneContext s(scene,band_id,sound_id,num_samples,absorption_factor,dry_level,keyframe_id);
 				scs.push_back(s);
 			}
+			// If we are only here to calculate the T60 reverberation time
+			// we are only going to render the first keyframe.
+			if ( calc_T60 ) break;
 		}
+		// If we are only here to calculate the T60 reverberation time
+		// we are only going to render the first sound file encountered.
+		if ( calc_T60 ) break;
 	}
 
 	if ( max_threads > 0 )
@@ -180,7 +197,7 @@ int Render(std::string filename, int id) {
 
 	// Calculate max response
 	float max = 0.0f;
-	for( std::vector<SceneContext>::iterator it = scs.begin(); it != scs.end(); ++it ) {
+	for( std::vector<SceneContext>::const_iterator it = scs.begin(); it != scs.end(); ++it ) {
 		for ( std::vector<Recorder*>::const_iterator rit = it->recorders.begin(); rit != it->recorders.end(); ++ rit ) {
 			Recorder* r1 = *rit;
 			r1->Power(0.335f);
@@ -211,8 +228,49 @@ int Render(std::string filename, int id) {
 	}	
 
 	const bool noprocess = Settings::IsSet("noprocessing") && Settings::GetBool("noprocessing");
-	if ( noprocess ) {
+	if ( noprocess || calc_T60 ) {
 		std::cout << std::endl << "Not processing data" << std::endl;
+
+		if ( calc_T60 ) {
+
+			// If we are only here to calculate the T60 reverberation time the rendered result
+			// does not need to be convoluted. Instead, the T60 is determined based on the
+			// rendered impulse response, as well as by the two well-known formulas Sabine
+			// and Norris-Eyring. These deal with the prediction of reverberation time on a
+			// statistical level. For a 'conventional' setup, the T60 that is calculated from
+			// the impulse response should not deviate too much from the statistical prediction.
+
+			const SceneContext sc = *scs.begin();
+			const Recorder* rec = *sc.recorders.begin();
+			const RecorderTrack* track = *rec->tracks.begin();
+			*calc_T60 = track->T60();
+
+			const Mesh* mesh = scene->meshes[0];
+			const float V = mesh->Volume();
+			const float A = mesh->TotalAbsorption();
+			const float S = mesh->Area();
+			const float a = mesh->AverageAbsorption();
+			const float m = absorption[1];
+
+			// Sabine:
+			//     0.1611 V
+			// T = -------
+			//        A
+			//
+			// Norris-Eyring:
+			//     -0.1611 V 
+			// T = ---------
+			//     S ln(1-a)
+
+			if ( T60_Sabine )
+				*T60_Sabine = 0.1611f*V/A;
+			if ( T60_Eyring )
+				*T60_Eyring = -0.1611f*V/(S*log(1.0f-a));
+		}
+
+		Datatype::Dispose();
+		Keyframes::Dispose();
+		delete scene;
 		return 0;
 	}
 
@@ -306,24 +364,37 @@ int Render(std::string filename, int id) {
 int main(int argc, char** argv) {
 	std::cout << HEADER << std::endl << std::endl << std::endl;
 	std::cout << std::setprecision(3) << std::fixed;
-	int id = -1;
 	for ( int i = 0; i < argc; i ++ ) {
-		std::string s(argv[i]);
-		if ( s == "render" && ((i+1)<argc) ) {
-			std::string filename(argv[i+1]);
+		const std::string cmd(argv[i]);
+		const std::string arg1 = ((i+1)<argc) ? std::string(argv[i+1]) : "";
+		const std::string arg2 = ((i+2)<argc) ? std::string(argv[i+2]) : "";
+		if ( cmd == "render" && !arg1.empty() ) {
 			int ret_value = 1;
 			try {
-				ret_value = Render(filename,id);
+				ret_value = Render(arg1);
 			} catch ( std::exception& e ) {
 				std::cout << std::endl << "Error: " << e.what() << std::endl << std::endl;
 			}
 			std::cout << "Press a key to exit..." << std::endl;
 			std::cin.get();
 			return ret_value;
+		} else if ( cmd == "calc" && arg1 == "T60" && !arg2.empty() ) {
+			float T60_value, T60_sabine, T60_eyring;
+			int ret_value = 1;
+			try {
+				ret_value = Render(arg2,&T60_value,&T60_sabine,&T60_eyring);
+				std::cout << "T60_ear   : " << std::setprecision(9) << std::fixed << T60_value << "s" << std::endl;
+				std::cout << "T60_sabine: " << std::setprecision(9) << std::fixed << T60_sabine << "s" << std::endl;
+				std::cout << "T60_eyring: " << std::setprecision(9) << std::fixed << T60_eyring << "s" << std::endl;
+			} catch ( std::exception& e ) {
+				std::cout << std::endl << "Error: " << e.what() << std::endl << std::endl;				
+			}
+			return ret_value;
 		}
 	}
 	std::cout << "Usage:" << std::endl 
-		<< " EAR render <filename>" << std::endl;
+		<< " EAR render <filename>" << std::endl
+		<< " EAR calc T60 <filename>" << std::endl;
 }
 
 boost::mutex MonoRecorder::mutex;

src/Mesh.cpp

@@ -72,16 +72,19 @@ Mesh* Mesh::Empty() {
 
 Mesh::Mesh(bool from_file) {
 	total_area = 0;
+	total_weighted_area = 0;
 	if ( ! from_file ) return;
 	Read(false);
 	assertid("MESH");		
 	std::string m = ReadString();
 	material = materials[m];
+	const float absorption = 1.0f - material->absorption_coefficient[1];
 	while ( Datatype::PeakId() == "tri " ) {
 		Triangle* tri = new Triangle();
 		tri->m = material;
 		tris.push_back(tri);
 		total_area += tri->area;
+		total_weighted_area += tri->area * absorption;
 	}
 	has_boundingbox = false;
 	BoundingBox();
@@ -91,6 +94,8 @@ Mesh::Mesh(bool from_file) {
 	ss << Datatype::prefix << "Mesh \r\n" << indent << " +- faces: " << tris.size() << std::endl << indent << " +- material: '" << material->name << "'" << std::endl;
 	ss << indent << " +- bounds: (" << xmin << ", " << ymin << ", " << zmin << ") - (" << xmax << ", " << ymax << ", " << zmax << ")" << std::endl;
 	ss << indent << " +- surface area: " << total_area << std::endl;
+	ss << indent << " +- total absorption: " << total_weighted_area << std::endl;
+	ss << indent << " +- volume: " << Volume() << std::endl;
 	if ( indent.size() ) {
 		Datatype::stringblock = ss.str();
 	} else {
@@ -144,4 +149,26 @@ void Mesh::SamplePoint(gmtl::Point3f& p, gmtl::Vec3f& n) {
 	}
 }
 
+float Mesh::Area() const {
+	return total_area;
+}
+
+float Mesh::TotalAbsorption() const {
+	return total_weighted_area;
+}
+
+float Mesh::Volume() const {
+	// http://stackoverflow.com/questions/1406029/how-to-calculate-the-volume-of-a-3d-mesh-object-the-surface-of-which-is-made-up
+	float volume = 0.0f;
+	std::vector<Triangle*>::const_iterator it;
+	for( it = tris.begin(); it != tris.end(); it ++ ) {
+		volume += (*it)->SignedVolume();
+	}
+	return volume;
+}
+
+float Mesh::AverageAbsorption() const {
+	return total_weighted_area / total_area;
+}
+
 std::map<std::string,Material*> Mesh::materials;

src/Mesh.h

@@ -42,6 +42,7 @@ class Mesh : public Datatype {
 private:
 	bool has_boundingbox;
 	float total_area;
+	float total_weighted_area;
 public:
 	std::vector<Triangle*> tris;
 	Material* material;
@@ -58,6 +59,22 @@ class Mesh : public Datatype {
 	static Mesh* Empty();
 	void Combine(Mesh* m);
 	void BoundingBox();
+	/// Returns the surface area of the mesh, useful for example to determine the T60
+	/// reverberation time using Sabine, Eyring or Millington-Sette.
+	float Area() const;
+	/// Returns the surface area of the mesh times the average absorption of the
+	/// surfaces, commonly called the Total Absorption measured in Sabins. Useful 
+	/// for example to determine the T60 reverberation time using Sabine, Eyring 
+	/// or Millington-Sette.
+	float TotalAbsorption() const;
+	/// Calculates the internal volume of the mesh. In case of a non-manifold or open
+	/// mesh, this function returns wrong results. For example to determine the T60
+	/// reverberation time using Sabine, Eyring or Millington-Sette.
+	float Volume() const;
+	/// Calculates the internal volume of the mesh. In case of a non-manifold or open
+	/// mesh, this function returns wrong results. For example to determine the T60
+	/// reverberation time using Sabine, Eyring or Millington-Sette.
+	float AverageAbsorption() const;
 };
 
 #endif

src/Recorder.h

@@ -200,6 +200,49 @@ class RecorderTrack : public FloatBuffer {
 			_this[i] += _other[i];
 		}
 	}
+	/// Returns the T60 reverberation time for the samples stored in this recorder track.
+	/// From http://en.wikipedia.org/wiki/Reverberation
+	/// T60 is the time required for reflections of a direct sound to decay by 60 dB below 
+	/// the level of the direct sound.
+	float T60() const {
+		const float attenuation_db = 60.0f;
+		const float attenuation_gain = powf(10.0f,attenuation_db/20.0f);
+
+		float direct_intensity = 0.0f;
+		float min_gain = 0.0f;
+		int last_significant_offset;
+		int direct_sound_offset;
+
+		const RecorderTrack& _this = *this;	
+
+		float previous_sample = -1.0f;
+		bool inside_indirect_lobe = false;
+		for ( unsigned int j = first_sample; j < real_length; ++ j ) {
+			const float sample = _this[j];
+			if ( inside_indirect_lobe ) {
+				if ( sample > min_gain ) {
+					last_significant_offset = j;
+				}
+			} else if ( sample < previous_sample ) {
+				// This is a rather silly way to determine the end of the direct sound
+				// field, for it may not even been present in this track and it is 
+				// explicitely calculated seperately from the reflections anyway in the
+				// rendering function, but this information is no longer available at
+				// this stage.
+				inside_indirect_lobe = true;
+				direct_intensity = previous_sample;
+				min_gain = direct_intensity / attenuation_gain;
+				direct_sound_offset = j;
+			}
+			previous_sample = sample;
+		}
+
+		const int reverberation_length = last_significant_offset - direct_sound_offset;
+		return (float)reverberation_length/44100.0f;
+	}
+	/// TODO: It would be great to implement other statistics as well. For example EDT
+	/// (Early Decay Time) & STI (Speach Transmission Index). The StereoRecorder class
+	/// could for example implement the IACC (Inter Aural Cross Correlation).
 };
 
 /// This class is the abstract base class for all classes of listeners. It defines

src/Scene.cpp

@@ -124,7 +124,7 @@ void Scene::Render(int band, int sound, float absorbtion_factor, int num_samples
 		Material* mat = 0;
 		BounceType bt;		
 
-		for( int num_bounces = 0; num_bounces < 200; num_bounces ++ ) {
+		for( int num_bounces = 0; num_bounces < 500; num_bounces ++ ) {
 			surface_normal = 0;
 			mat = 0;
 			if ( ! sound_ray ) {
@@ -212,6 +212,8 @@ void Scene::Render(int band, int sound, float absorbtion_factor, int num_samples
 
 			}
 
+			// Arbitrary constant, ideally this would be determined based
+			// on some heuristics or previously collected samples.
 			if ( sample_intensity < 0.00000001 ) break;
 
 			previous_ray_direction = gmtl::makeNormal(sound_ray->mDir);
@@ -230,6 +232,9 @@ void Scene::Render(int band, int sound, float absorbtion_factor, int num_samples
 		rec->Multiply(1.0f / amount);
 
 		// The direct sound lobe is added...
+		// Ideally this lobe would be stored separately from the rest of the samples, it could then
+		// be subject of dopler-effect calculations for example and would ease the calculation of
+		// some of the statistical properties of the rendered impulse response.
 		if ( !currentSound->isMeshSource() ) {
 		const gmtl::Point3f listener_location = rec->getLocation(keyframeID);
 		gmtl::LineSegf* ls = Connect(&listener_location,sfloc);
@@ -254,5 +259,7 @@ Scene::~Scene() {
 	for ( std::vector<AbstractSoundFile*>::const_iterator it = sources.begin(); it != sources.end(); ++ it ) {
 		delete *it;
 	}
-	delete meshes[0];
+	for ( std::vector<Mesh*>::const_iterator it = meshes.begin(); it != meshes.end(); ++ it ) {
+		delete *it;
+	}
 }

src/Triangle.cpp

@@ -50,4 +50,18 @@ void Triangle::SamplePoint(gmtl::Point3f& P) {
 	const gmtl::Point3f& B = (*this)[1];
 	const gmtl::Point3f& C = (*this)[2];
 	P = (1.0f - sr1) * A + (sr1 * (1.0f - r2)) * B + (sr1 * r2) * C;
+}
+
+float Triangle::SignedVolume() const {
+	// http://stackoverflow.com/questions/1406029/how-to-calculate-the-volume-of-a-3d-mesh-object-the-surface-of-which-is-made-up
+	const gmtl::Point3f& p1 = mVerts[0];
+	const gmtl::Point3f& p2 = mVerts[1];
+	const gmtl::Point3f& p3 = mVerts[2];
+	const float v321 = p3[0]*p2[1]*p1[2];
+    const float v231 = p2[0]*p3[1]*p1[2];
+    const float v312 = p3[0]*p1[1]*p2[2];
+    const float v132 = p1[0]*p3[1]*p2[2];
+    const float v213 = p2[0]*p1[1]*p3[2];
+    const float v123 = p1[0]*p2[1]*p3[2];
+    return (1.0f/6.0f)*(-v321 + v231 + v312 - v132 - v213 + v123);
 }

src/Triangle.h

@@ -39,6 +39,7 @@ class Triangle : public gmtl::Trif, public Datatype {
 	Triangle(const gmtl::Point3f& a,const gmtl::Point3f& b,const gmtl::Point3f& c);
 	Triangle();
 	void SamplePoint(gmtl::Point3f& p);
+	float SignedVolume() const;
 };
 
 #endif