implemented merging

.gitignore

@@ -1 +1,6 @@
 CMakeFiles
+Makefile
+CMakeCache.txt
+hdrmerge
+cmake_install.cmake
+rawspeed/librawspeed.a

CMakeLists.txt

@@ -1,6 +1,7 @@
 cmake_minimum_required(VERSION 2.6)
 project(hdrmerge)
 
+set(CMAKE_BUILD_TYPE Release)
 set(CMAKE_C_COMPILER /usr/local/bin/gcc)
 set(CMAKE_C_COMPILER /usr/local/bin/gcc)
 
@@ -21,7 +22,7 @@ include_directories(${OPENEXR_INCLUDE_PATHS})
 include_directories(${EXIV2_INCLUDE_DIR})
 
 # Compile with C++11 features
-set (CMAKE_CXX_FLAGS "-std=c++11 ${CMAKE_CXX_FLAGS}")
+set (CMAKE_CXX_FLAGS "-std=c++11 -g -fno-omit-frame-pointer ${CMAKE_CXX_FLAGS}")
 
 add_subdirectory(rawspeed)
 add_executable(hdrmerge input.cpp util.cpp main.cpp)

hdrmerge.h

@@ -20,7 +20,7 @@ struct Exposure {
 	std::string filename;
 	float exposure;
 	float shown_exposure;
-	float *image;
+	uint16_t *image;
 
 	inline Exposure(const std::string &filename)
 	 : filename(filename), exposure(-1), image(NULL) { }
@@ -45,7 +45,15 @@ struct ExposureSeries {
 	std::vector<Exposure> exposures;
 	StringMap metadata;
 	size_t width, height;
-	float saturation;
+	uint16_t blacklevel, saturation, whitepoint;
+	float *image;
+
+	inline ExposureSeries() : image(NULL) { }
+
+	~ExposureSeries() {
+		if (image)
+			delete image;
+	}
 
 	/**
 	 * Add a file to the exposure series (or, optionally, a sequence
@@ -72,16 +80,16 @@ struct ExposureSeries {
 	void load();
 
 	/// Evaluate a given pixel of an exposure series
-	inline float eval(size_t img, size_t x, size_t y) const {
+	inline uint16_t &eval(size_t img, size_t x, size_t y) const {
 		return exposures[img].image[x+width*y];
 	}
 
 	/// Evaluate a given pixel of an exposure series (Bayer grid)
-	inline float eval(size_t img, size_t x, size_t y, int ch) const {
+/*	inline float &eval(size_t img, size_t x, size_t y, int ch) const {
 		x = x*2 + ch&1;
 		y = y*2 + (ch&2) >> 1;
 		return exposures[img].image[x+width*y];
-	}
+	}*/
 
 	/// Return the number of exposures
 	inline size_t size() const {

input.cpp

@@ -59,6 +59,15 @@ void ExposureSeries::add(const char *fmt) {
 	}
 }
 
+float exposureTime(float shutterSpeedValue) {
+	/* lifted from libexiv2 */
+	double tmp = std::exp(std::log(2.0) * shutterSpeedValue);
+	if (tmp > 1)
+		return 1.0f / ((int) (tmp + 0.5));
+	else
+		return (int) (1/tmp + 0.5);
+}
+
 void ExposureSeries::check() {
 	float isoSpeed = -1, aperture = -1;
 
@@ -91,7 +100,9 @@ void ExposureSeries::check() {
 		it = exifData.findKey(Exiv2::ExifKey("Exif.Photo.ShutterSpeedValue"));
 		if (it == exifData.end())
 			throw std::runtime_error("\"" + exp.filename + "\": could not extract the exposure time!");
+
 		exp.exposure = std::pow(2, -it->toFloat());
+		exp.exposure = exposureTime(it->toFloat());
 
 		it = Exiv2::exposureTime(exifData);
 		if (it == exifData.end())
@@ -160,8 +171,8 @@ void ExposureSeries::check() {
 			return a.exposure == b.exposure;
 	});
 
-	if (it != exposures.end())
-		throw std::runtime_error((boost::format("Duplicate exposure time: %1%") % it->toString()).str());
+//	if (it != exposures.end())
+//		throw std::runtime_error((boost::format("Duplicate exposure time: %1%") % it->toString()).str());
 
 	cout << "Collected " << metadata.size() << " metadata entries." << endl;
 }
@@ -207,21 +218,15 @@ void ExposureSeries::load() {
 		if (i == 0) {
 			this->width = width;
 			this->height = height;
+			this->blacklevel = raw->blackLevel;
+			this->whitepoint = raw->whitePoint;
 		}
 
 		uint16_t *data = (uint16_t *) raw->getData(0, 0);
+		uint16_t *image = new uint16_t[width*height];
 
-		uint16_t offset = raw->blackLevel;
-		float factor = 1.0f / (raw->whitePoint - offset);
-
-		float *image = new float[width*height];
-		for (int y=0; y<height; ++y) {
-			uint16_t *src = data + y*pitch;
-			float *dst = image + y*width;
-
-			for (int x=0; x<width; ++x)
-				*dst++ = (*src++ - offset) * factor;
-		}
+		for (int y=0; y<height; ++y)
+			memcpy(image+y*width, data+y*pitch, sizeof(uint16_t)*width);
 
 		exposures[i].image = image;
 
@@ -233,17 +238,18 @@ void ExposureSeries::load() {
 	}
 
 	cout << " done (" << width << "x" << height << ", using "
-		 << (width*height*sizeof(float) * exposures.size()) / (float) (1024*1024)
+		 << (width*height*sizeof(uint16_t) * exposures.size()) / (float) (1024*1024)
 		 << " MiB of memory)" << endl;
 
 	/* Determine the value of a pixel considered to be overexposured */
 	size_t npix = width*height;
-	float *temp = new float[npix];
-	memcpy(temp, exposures[size()-1].image, npix*sizeof(float));
+	uint16_t *temp = new uint16_t[npix];
+	memcpy(temp, exposures[size()-1].image, npix*sizeof(uint16_t));
 	size_t percentile = (size_t) (npix*0.999);
 	std::nth_element(temp, temp+percentile, temp+npix);
 	saturation = *(temp+percentile);
 	delete[] temp;
 
-	cout << "Saturation detected to be around " << saturation << "." << endl;
+	cout << "Saturation determined to be around "
+		<< (saturation-blacklevel) / (float) (whitepoint-blacklevel) * 100 << "\% of the dynamic range." << endl;
 }

main.cpp

@@ -1,43 +1,100 @@
 #include "hdrmerge.h"
 
-void merge(ExposureSeries &es) {
-	int x = (int) (randf() * es.width);
-	int y = (int) (randf() * es.height);
+float compute_weight(uint16_t value, uint16_t blacklevel, uint16_t saturation) {
+	const float alpha = -1.0f / 10.0f;
+	const float beta = 1.0f / std::exp(4.0f*alpha);
 
-	cout << "x=" << x << ", y=" << y << endl;
+	float scaled = (value - blacklevel) / (float) (saturation - blacklevel);
 
-	cout << "[";
-	for (int i=0; i<es.size(); ++i) {
-		cout << es.eval(i, x, y);
-		if (i+1 < es.size())
-			cout << ", ";
-	}
-	cout << "]" << endl;
+	if (scaled <= 0 || scaled >= 1)
+		return 0;
+
+	return beta * std::exp(alpha * (1/scaled + 1/(1-scaled)));
+}
+
+float compute_weight(float scaled) {
+	const float alpha = -1.0f / 10.0f;
+	const float beta = 1.0f / std::exp(4.0f*alpha);
+
+	if (scaled <= 0 || scaled >= 1)
+		return 0;
 
-	cout << "[";
-	for (int i=0; i<es.size(); ++i) {
-		cout << es.exposures[i].exposure;
-		if (i+1 < es.size())
-			cout << ", ";
+	return beta * std::exp(alpha * (1/scaled + 1/(1-scaled)));
+}
+
+void merge(ExposureSeries &es) {
+	cout << "Merging " << es.size() << " exposures .." << endl;
+	es.image = new float[es.width * es.height];
+
+	/* Precompute some tables for weights and normalized pixel values */
+	float weight_tbl[0xFFFF], value_tbl[0xFFFF];
+	for (int i=0; i<0xFFFF; ++i) {
+		weight_tbl[i] = compute_weight((uint16_t) i, es.blacklevel, es.saturation);
+		value_tbl[i] = (float) (i - es.blacklevel) / (float) (es.whitepoint - es.blacklevel);
 	}
-	cout << "]" << endl;
 
+
+//	#pragma omp parallel for schedule(dynamic, 1)
+	for (int y=0; y<es.height; ++y) {
+		uint32_t offset = y * es.width;
+		for (int x=0; x<es.width; ++x) {
+			float value = 0, total_exposure = 0;
+			offset = 11002202;
+
+			cout << "Offset=" << offset << endl;
+			for (int img=0; img<es.size(); ++img) {
+				uint16_t pxvalue = es.exposures[img].image[offset];
+				float weight = weight_tbl[pxvalue];
+				value += value_tbl[pxvalue] * weight;
+				total_exposure += es.exposures[img].exposure * weight;
+
+				cout << img << ": " << pxvalue << " => weight=" << weight << ", value=" << value << endl;
+			}
+
+			if (total_exposure > 0)
+				value /= total_exposure;
+
+			cout << " ====>> reference = " << value << endl;
+			float reference = value;
+			value = total_exposure = 0;
+
+			float blacklevel = es.blacklevel, scale = es.whitepoint - es.blacklevel;
+
+			for (int img=0; img<es.size(); ++img) {
+				float predicted = reference * es.exposures[img].exposure * scale + blacklevel;
+				if (predicted <= 0 || predicted >= 65535.0f)
+					continue;
+				uint16_t predicted_pxvalue = (uint16_t) (predicted + 0.5f);
+				uint16_t pxvalue = es.exposures[img].image[offset];
+				float weight = weight_tbl[predicted_pxvalue];
+
+				value += value_tbl[pxvalue] * weight;
+				total_exposure += es.exposures[img].exposure * weight;
+//				cout << img << ": " << value_tbl[pxvalue] << " (pred=" << reference * es.exposures[img].exposure << ")  => weight=" << weight << ", value=" << value << endl;
+				cout << pxvalue << ", " << predicted_pxvalue << "; ";
+			}
+			cout << endl;
+
+			if (total_exposure > 0)
+				value /= total_exposure;
+			cout << " ====>> reference = " << value << endl;
+
+			es.image[offset++] = value;
+			exit(-1);
+		}
+	}
 }
 
 int main(int argc, char **argv) {
 	ExposureSeries es;
 //	es.add("test-%02i.cr2");
-	es.add("/mnt/raid0/layered2/meas-%05i-00000.cr2");
+	es.add("/mnt/raid0/layered2/meas2-%05i-00000.cr2");
 	es.check();
 	if (es.size() == 0)
 		throw std::runtime_error("No input found / list of exposures to merge is empty!");
 	es.load();
+	merge(es);
 
-	for (int i=0; i<10; ++i) {
-		merge(es);
-	}
-
-
-	writeOpenEXR("test.exr", es.width, es.height, 1, es.exposures[2].image, es.metadata, false);
+	writeOpenEXR("test.exr", es.width, es.height, 1, es.image, es.metadata, false);
 	return 0;
 }

rawspeed/cmake_install.cmake

@@ -12,7 +12,7 @@ IF(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME)
     STRING(REGEX REPLACE "^[^A-Za-z0-9_]+" ""
            CMAKE_INSTALL_CONFIG_NAME "${BUILD_TYPE}")
   ELSE(BUILD_TYPE)
-    SET(CMAKE_INSTALL_CONFIG_NAME "")
+    SET(CMAKE_INSTALL_CONFIG_NAME "Release")
   ENDIF(BUILD_TYPE)
   MESSAGE(STATUS "Install configuration: \"${CMAKE_INSTALL_CONFIG_NAME}\"")
 ENDIF(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME)

util.cpp

@@ -13,6 +13,7 @@ void writeOpenEXR(const std::string &filename, size_t w, size_t h, int nChannels
 
 	Imf::ChannelList &channels = header.channels();
 
+	cout << "Writing " << filename << " .. " << endl;
 	if (nChannels == 3) {
 		if (writeHalf) {
 			channels.insert("R", Imf::Channel(Imf::HALF));