Actually build an HDR-image and do something with it

.gitignore

@@ -1,3 +1,4 @@
 haader
 .*.swp
 *.o
+*.ppm

Makefile

@@ -1,15 +1,15 @@
 CXXFLAGS= -std=c++98 -pedantic -Wall -O3
 LIBS= -ljpeg
-OBJS= haader.o image.o
-
-run: haader
-	./haader
+OBJS= haader.o image.o exif.o csv_reader.o
 
 haader: main.cpp $(OBJS)
 	$(CXX) main.cpp -o haader $(OBJS) $(CXXFLAGS) $(LIBS)
 
+run: haader
+	./haader
+
 %.o: %.c %.h
-	gcc $< -c -o $@ $(CFLAGS)
+	$(CXX) $< -c -o $@ $(CFLAGS)
 
 clean:
 	rm -f haader $(OBJS)

csv_reader.cpp

@@ -0,0 +1,382 @@
+#include "csv_reader.h"
+
+#include <errno.h>
+#include <limits.h>
+#include <stdlib.h>
+
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+
+CSVField::CSVField(int start, int length, bool hasQuotes): start(start),
+                                                           length(length),
+                                                           hasQuotes(hasQuotes)
+{
+}
+
+CSVRow::CSVRow(string *row,
+               const std::string * const filePath,
+               long lineNumber,
+               char delimiter,
+               char quote): m_row(row),
+                            m_filePath(filePath),
+                            m_lineNumber(lineNumber),
+                            m_quote(quote)
+{
+    parse(delimiter, quote);
+}
+
+string CSVRow::toString() const {
+    stringstream ss;
+
+    ss << "CSVRow (" << endl;
+
+    if (m_lineNumber > 0) {
+        ss << "  line #: \"" << m_lineNumber << "\"" << endl;
+    }
+
+    if (m_filePath) {
+        ss << "  path: \"" << *m_filePath << "\"" << endl;
+    }
+
+    ss << "  value: \"" << *m_row << "\"" << endl;
+
+    for (unsigned int i = 0; i < m_fields.size(); i++) {
+        string row;
+        getFieldAsString(i, &row);
+        ss << "  field " << i << ": '" << row << "'" << endl;
+    }
+
+    ss << ")";
+    return ss.str();
+}
+
+static bool string_to_long(const string &s, long *out) {
+    const char *cstr = s.c_str();
+    char *end = NULL;
+
+    errno = 0;
+
+    long x = strtol(cstr, &end, 10);
+
+    if (errno != 0) {
+        return false;
+    } else if (end == &cstr[s.size()]) {
+        *out = x;
+        return true;
+    } else {
+        return false;
+    }
+}
+
+static bool cstring_to_long(const char *start, const char *end, long *out) {
+    char *parse_end = NULL;
+
+    errno = 0;
+
+    long x = strtol(start, &parse_end, 10);
+
+    if (errno != 0) {
+        return false;
+    } else if (parse_end == end) {
+        *out = x;
+        return true;
+    } else {
+        return false;
+    }
+}
+
+static bool string_to_double(const string &s, double *out) {
+    const char *cstr = s.c_str();
+    char *end = NULL;
+
+    errno = 0;
+
+    double x = strtod(cstr, &end);
+
+    if (errno != 0) {
+        return false;
+    } else if (end == &cstr[s.size()]) {
+        *out = x;
+        return true;
+    } else {
+        return false;
+    }
+}
+
+static bool cstring_to_double(const char *start, const char *end, double *out) {
+    char *parse_end = NULL;
+
+    errno = 0;
+
+    double x = strtod(start, &parse_end);
+
+    if (errno != 0) {
+        return false;
+    } else if (parse_end == end) {
+        *out = x;
+        return true;
+    } else {
+        return false;
+    }
+}
+
+bool CSVRow::getFieldAsString(unsigned int fieldIndex, string *out) const {
+    if (fieldIndex < m_fields.size()) {
+        CSVField field = m_fields[fieldIndex];
+
+        if (field.length == 0) {
+            *out = "";
+        } else {
+            string s;
+
+            if (field.hasQuotes) {
+                s = m_row->substr(field.start, field.length);
+                s = unquote(s, m_quote);
+            } else {
+                s = m_row->substr(field.start, field.length);
+            }
+
+            // remove newlines
+            s.erase(std::remove(s.begin(), s.end(), '\n'), s.end());
+            s.erase(std::remove(s.begin(), s.end(), '\r'), s.end());
+
+            *out = s;
+        }
+
+        return true;
+    } else {
+        return false;
+    }
+}
+
+bool CSVRow::getFieldAsLong(unsigned int fieldIndex, long *out) const {
+    if (fieldIndex >= m_fields.size()) {
+        return false;
+    }
+
+    if (m_fields[fieldIndex].hasQuotes) {
+        string field;
+        if (!getFieldAsString(fieldIndex, &field)) {
+            return false;
+        }
+
+        return string_to_long(field, out);
+    } else {
+        const CSVField *field = &m_fields[fieldIndex];
+        int sep_index = field->start + field->length;
+        char sep = (*m_row)[sep_index];
+        (*m_row)[sep_index] = '\0';
+
+        bool ok = cstring_to_long(m_row->c_str() + field->start, m_row->c_str() + sep_index, out);
+
+        (*m_row)[sep_index] = sep;
+
+        return ok;
+    }
+}
+
+bool CSVRow::getFieldAsInt(unsigned int fieldIndex, int *out) const {
+    if (fieldIndex >= m_fields.size()) {
+        return false;
+    }
+
+    long x;
+    if (!getFieldAsLong(fieldIndex, &x)) {
+        return false;
+    }
+
+    if (x < INT_MIN || x > INT_MAX) {
+        // underflow or overflow
+        return false;
+    }
+
+    *out = (int)x;
+
+    return true;
+}
+
+bool CSVRow::getFieldAsDouble(unsigned int fieldIndex, double *out) const {
+    if (fieldIndex >= m_fields.size()) {
+        return false;
+    }
+
+    if (m_fields[fieldIndex].hasQuotes) {
+        string field;
+        if (!getFieldAsString(fieldIndex, &field)) {
+            return false;
+        }
+
+        return string_to_double(field, out);
+    } else {
+        const CSVField *field = &m_fields[fieldIndex];
+        int sep_index = field->start + field->length;
+        char sep = (*m_row)[sep_index];
+        (*m_row)[sep_index] = '\0';
+
+        bool ok = cstring_to_double(m_row->c_str() + field->start, m_row->c_str() + sep_index, out);
+
+        (*m_row)[sep_index] = sep;
+
+        return ok;
+    }
+}
+
+string CSVRow::getFilePath() const {
+    if (m_filePath) {
+        return string(*m_filePath);
+    } else {
+        return "";
+    }
+}
+
+long CSVRow::getLineNumber() const {
+    return m_lineNumber;
+}
+
+unsigned int CSVRow::getNumberOfFields() const {
+    return m_fields.size();
+}
+
+void CSVRow::parse(char delimiter, char quote) {
+    m_fields.clear();
+
+    //TODO: properly handle unicode strings (could use QString::fromUtf8 and iterate over that)
+    int start = 0;
+    int length = -1;
+    bool in_quotes = false;
+    bool has_quotes = false;
+
+    unsigned int size = m_row->size();
+
+    // handle trailing \n and \r chars
+    while (size > 0 && ((*m_row)[size - 1] == '\n' || (*m_row)[size - 1] == '\r')) {
+        size--;
+    }
+
+    for (unsigned int i = 0; i < size; i++) {
+        length++;
+        char c = (*m_row)[i];
+
+        if (c == delimiter && !in_quotes) {
+            m_fields.push_back(CSVField(start, length, has_quotes));
+            start = i + 1;
+            length = -1;
+            has_quotes = false;
+        }
+
+        if (c == quote) {
+            has_quotes = true;
+            if (in_quotes) {
+                in_quotes = false;
+            } else {
+                in_quotes = true;
+            }
+        }
+    }
+    if (!(start == 0 && length == -1)) {
+        length++;
+        m_fields.push_back(CSVField(start, length, has_quotes));
+    }
+}
+
+string CSVRow::unquote(const string &fieldStr, char quoteChar) {
+    string ret;
+    bool in_quotes = false;
+    bool last_char_was_quote = false;
+
+    //TODO: properly handle unicode strings here too
+    for (unsigned int i = 0; i < fieldStr.size(); i++) {
+        char c = fieldStr[i];
+
+        if (in_quotes) {
+            if (c == quoteChar) {
+                in_quotes = false;
+                last_char_was_quote = true;
+            } else {
+                ret.push_back(c);
+                last_char_was_quote = false;
+            }
+        } else {
+            if (c == quoteChar) {
+                in_quotes = true;
+                if (last_char_was_quote) {
+                    ret.push_back(quoteChar);
+                }
+            }
+            last_char_was_quote = false;
+        }
+    }
+
+    return ret;
+}
+
+string CSVRow::quote(const string &fieldStr, char quoteChar) {
+    string ret;
+
+    //TODO: properly handle unicode strings here too
+    //TODO: optimize using string::find
+    for (unsigned int i = 0; i < fieldStr.size(); i++) {
+        char c = fieldStr[i];
+
+        if (c == quoteChar) {
+            ret.append(2, quoteChar);
+        } else {
+            ret.push_back(c);
+        }
+    }
+
+    return ret;
+}
+
+bool CSVReader::readFromFile(const std::string filePath,
+                             CSVReader::RowCallback callback,
+                             void *userData,
+                             char delimiter,
+                             char quote)
+{
+    ifstream stream;
+    stream.open(filePath.c_str(), ios::in);
+    if (stream) {
+        bool ok = CSVReader::readFromStream(stream,
+                                            callback,
+                                            userData,
+                                            &filePath,
+                                            delimiter,
+                                            quote);
+        stream.close();
+        return ok;
+    } else {
+        cerr << "Error: Cannot read file \"" << filePath << "\"" << endl;
+        return false;
+    }
+
+    return true;
+}
+
+
+bool CSVReader::readFromStream(std::istream &stream,
+                               CSVReader::RowCallback callback,
+                               void *userData,
+                               const string * const filePath,
+                               char delimiter,
+                               char quote)
+{
+    string line = "";
+    long nr = 0;
+    while (getline(stream, line, '\n')) {
+        nr++;
+        if (line.size() > 0) {
+            CSVRow row(&line, filePath, nr, delimiter, quote);
+            if (!callback(row, userData)) {
+                return false;
+            }
+        }
+    }
+
+    return true;
+}

csv_reader.h

@@ -0,0 +1,115 @@
+#ifndef CSVREADER_H
+#define CSVREADER_H
+
+#include <string>
+#include <vector>
+#include <istream>
+
+
+// A Field in a CSVRow
+class CSVField {
+public:
+    CSVField(int start, int length, bool hasQuotes);
+
+    // byte offset
+    int start;
+
+    // length in bytes
+    int length;
+
+    // contains quotes (use CSVRow::unquote)
+    bool hasQuotes;
+};
+
+// A Row in a CSV Document
+class CSVRow {
+public:
+    CSVRow(std::string *row,
+           const std::string * const filePath=0,
+           long lineNumber=-1,
+           char delimiter=',',
+           char quote='"');
+
+    // Return a human readable representation of a CSVRow.
+    std::string toString() const;
+
+    // Get the value of a field.
+    // fieldIndex starts at 0.
+    // Returns false if fieldIndex is out of bounds or
+    // field could not be converted to the specified type.
+    bool getFieldAsString(unsigned int fieldIndex, std::string *out) const;
+    bool getFieldAsLong(unsigned int fieldIndex, long *out) const;
+    bool getFieldAsInt(unsigned int fieldIndex, int *out) const;
+    bool getFieldAsDouble(unsigned int fieldIndex, double *out) const;
+
+    std::string getFilePath() const;
+
+    // Get the line number
+    // starting at 1 for the first line.
+    // Return a value < 0 if there is no line number specified.
+    long getLineNumber() const;
+
+    unsigned int getNumberOfFields() const;
+
+    // Get the original string that was quoted using some of the rules from RFC 4180
+    // Examples: unquote("'abc'", "'") == "abc"
+    //           unquote("'ab''c'", "'") == "ab'c"
+    // Fields containing line breaks are not supported.
+    static std::string unquote(const std::string &fieldStr, char quoteChar);
+
+    // Quotes a string using some of the rules from RFC 4180
+    // Example: quote("ab'c", "'") == "ab''c"
+    // Fields containing line breaks are not supported.
+    static std::string quote(const std::string &fieldStr, char quoteChar);
+
+private:
+    // fill m_fields
+    void parse(char delimiter, char quote);
+
+    // raw bytes from CSV document
+    std::string *m_row;
+
+    // each field
+    std::vector<CSVField> m_fields;
+
+    // file path of the CSV document or NULL
+    const std::string * const m_filePath;
+
+    long m_lineNumber;
+
+    // quote character
+    char m_quote;
+};
+
+
+// static functions to parse CSV data
+class CSVReader {
+public:
+
+    // callback function signature
+    typedef bool (*RowCallback) (const CSVRow &row, void *userData);
+
+    // Parse a CSV document given by filePath.
+    // Invokes callback function for each row with the supplied userData.
+    // Stop parsing if callback returns false.
+    // filePath: pointer to a file path or NULL if there is no file specified
+    // Returns false if file cannot be opened.
+    static bool readFromFile(const std::string filePath,
+                             RowCallback callback,
+                             void *userData,
+                             char delimiter=',',
+                             char quote='"');
+
+    // Parse CSV data from stream.
+    // Could be a std::stringstream or a std::fstream, etc.
+    // Invokes callback function for each row with the supplied userData.
+    // Stop parsing and return false if callback returns false.
+    static bool readFromStream(std::istream &stream,
+                               RowCallback callback,
+                               void *userData,
+                               const std::string * const filePath=0,
+                               char delimiter=',',
+                               char quote='"');
+};
+
+#endif //CSVREADER_H

exif.cpp

@@ -0,0 +1,342 @@
+#include "exif.h"
+
+#include <iostream>
+#include <fstream>
+#include <stdint.h>
+#include <math.h>
+
+using namespace std;
+
+
+urational_t::urational_t(): p(0), q(0)
+{}
+
+urational_t::urational_t(uint32_t p, uint32_t q): p(p), q(q)
+{}
+
+srational_t::srational_t(): p(0), q(0)
+{}
+
+srational_t::srational_t(int32_t p, int32_t q): p(p), q(q)
+{}
+
+
+ExifTags::ExifTags():
+    ifd0_parsed_ok(false),
+    exififd_parsed_ok(false),
+    exififd_offset(0)
+{
+}
+
+
+static uint16_t read_be_uint16(istream &stream) {
+    unsigned char buf[2];
+    stream.read((char*)buf, 2);
+    return ((uint16_t)buf[0]) * 256 + (uint16_t)buf[1];
+}
+
+static uint16_t read_le_uint16(istream &stream) {
+    unsigned char buf[2];
+    stream.read((char*)buf, 2);
+    return (uint16_t)buf[0] + ((uint16_t)buf[1]) * 256;
+}
+
+static uint16_t read_uint16(istream &stream, bool is_little_endian) {
+    if (is_little_endian) {
+        return read_le_uint16(stream);
+    } else {
+        return read_be_uint16(stream);
+    }
+}
+
+static uint32_t read_be_uint32(istream &stream) {
+    unsigned char buf[4];
+    stream.read((char*)buf, 4);
+    return ((uint32_t)buf[0]) * 16777216 +
+           ((uint32_t)buf[1]) * 65536 +
+           ((uint32_t)buf[2]) * 256 +
+           (uint32_t)buf[3];
+}
+
+static uint32_t read_le_uint32(istream &stream) {
+    unsigned char buf[4];
+    stream.read((char*)buf, 4);
+    return ((uint32_t)buf[3]) * 16777216 +
+           ((uint32_t)buf[2]) * 65536 +
+           ((uint32_t)buf[1]) * 256 +
+           (uint32_t)buf[0];
+}
+
+static uint32_t read_uint32(istream &stream, bool is_little_endian) {
+    if (is_little_endian) {
+        return read_le_uint32(stream);
+    } else {
+        return read_be_uint32(stream);
+    }
+}
+
+static int32_t read_be_int32(istream &stream) {
+    unsigned char buf[4];
+    stream.read((char*)buf, 4);
+    return (int32_t)(((uint32_t)buf[0]) * 16777216 +
+                     ((uint32_t)buf[1]) * 65536 +
+                     ((uint32_t)buf[2]) * 256 +
+                     (uint32_t)buf[3]);
+}
+
+static int32_t read_le_int32(istream &stream) {
+    unsigned char buf[4];
+    stream.read((char*)buf, 4);
+    return (int32_t)(((uint32_t)buf[3]) * 16777216 +
+                     ((uint32_t)buf[2]) * 65536 +
+                     ((uint32_t)buf[1]) * 256 +
+                     (uint32_t)buf[0]);
+}
+
+/*
+static int32_t read_int32(istream &stream, bool is_little_endian) {
+    if (is_little_endian) {
+        return read_le_int32(stream);
+    } else {
+        return read_be_int32(stream);
+    }
+}
+*/
+
+static urational_t read_urational(istream &stream, bool is_little_endian) {
+    if (is_little_endian) {
+        uint32_t p = read_le_uint32(stream);
+        uint32_t q = read_le_uint32(stream);
+        return urational_t(p, q);
+    } else {
+        uint32_t p = read_be_uint32(stream);
+        uint32_t q = read_be_uint32(stream);
+        return urational_t(p, q);
+    }
+}
+
+static srational_t read_srational(istream &stream, bool is_little_endian) {
+    if (is_little_endian) {
+        int32_t p = read_le_int32(stream);
+        int32_t q = read_le_int32(stream);
+        return srational_t(p, q);
+    } else {
+        int32_t p = read_be_int32(stream);
+        int32_t q = read_be_int32(stream);
+        return srational_t(p, q);
+    }
+}
+
+
+static bool read_exif_tags(istream &stream,
+                           unsigned int num_entries,
+                           ExifTags &tags,
+                           uint64_t tiff_start_pos,
+                           bool is_little_endian)
+{
+    for (unsigned int i = 0; i < num_entries; i++) {
+        uint16_t tag_num = read_uint16(stream, is_little_endian);
+        uint16_t data_format = read_uint16(stream, is_little_endian);
+        uint32_t num_components = read_uint32(stream, is_little_endian);
+        uint32_t tag_data = read_uint32(stream, is_little_endian);
+
+        if (tag_num == 34665) {
+            if (data_format == 4 && num_components == 1) {
+                if (tag_data > 65535) {
+                    cerr << "Exif IFD offset too big (" << tag_data << ")." << endl;
+                    return false;
+                }
+                tags.exififd_offset = tag_data;
+            } else {
+                cerr << "Exif IFD offset has strange format." << endl;
+                return false;
+            }
+
+        }
+
+        if (tag_num == 0x829a) {
+            uint64_t cur = stream.tellg();
+            stream.seekg(tiff_start_pos + tag_data, std::ios_base::beg);
+
+            tags.exposure_time = read_urational(stream, is_little_endian);
+
+            stream.seekg(cur, std::ios_base::beg);
+        }
+
+        if (tag_num == 0x9201) {
+            uint64_t cur = stream.tellg();
+            stream.seekg(tiff_start_pos + tag_data, std::ios_base::beg);
+
+            tags.shutter_speed_value = read_srational(stream, is_little_endian);
+
+            stream.seekg(cur, std::ios_base::beg);
+        }
+    }
+
+    return true;
+}
+
+static bool read_exif_ifd(istream &stream,
+                          uint64_t tiff_start_pos,
+                          uint16_t ifd_offset,
+                          ExifTags &tags,
+                          bool is_little_endian)
+{
+    if (ifd_offset == 0) {
+        return true;
+    }
+
+    {
+        uint64_t cur = stream.tellg();
+        if (cur < tiff_start_pos) {
+            cerr << "Read backwards?" << endl;
+            return false;
+        }
+
+        if (ifd_offset < (cur - tiff_start_pos)) {
+            cerr << "IFD offset too small." << endl;
+            return false;
+        }
+    }
+
+    stream.seekg(tiff_start_pos + ifd_offset, std::ios_base::beg);
+
+    uint16_t num_entries = read_le_uint16(stream);
+
+    return read_exif_tags(stream, num_entries, tags, tiff_start_pos, is_little_endian);
+}
+
+static bool parse_exif(istream &stream, uint16_t length, ExifTags &output) {
+    if (length < 8) {
+        return false;
+    }
+
+    {
+        char buf[6];
+        const char *compare = "Exif\0\0";
+        stream.read(buf, 6);
+        for (unsigned int i = 0; i < 6; i++) {
+            if (buf[i] != compare[i]) {
+                return false;
+            }
+        }
+    }
+
+    uint64_t tiff_start_pos = stream.tellg();
+
+    bool is_little_endian = false;
+    {
+
+        char buf[2];
+        stream.read(buf, 2);
+
+        if (buf[0] == 'I' && buf[1] == 'I') {
+            is_little_endian = true;
+        } else if (buf[0] == 'M' && buf[1] == 'M') {
+            is_little_endian = false;
+        } else {
+            cerr << "Endian fail" << endl;
+            return false;
+        }
+    }
+
+    {
+        uint16_t magic = read_uint16(stream, is_little_endian);
+
+        if (magic != 0x002a) {
+            cerr << "TIFF fail" << endl;
+            return false;
+        }
+    }
+
+    uint32_t ifd0_offset = read_uint32(stream, is_little_endian);
+
+    // Read IFD0
+    output.ifd0_parsed_ok = read_exif_ifd(stream, tiff_start_pos, ifd0_offset, output, is_little_endian);
+
+    if (output.exififd_offset != 0) {
+        // Read Exif IFD
+        output.exififd_parsed_ok = read_exif_ifd(stream, tiff_start_pos, output.exififd_offset, output, is_little_endian);
+    }
+
+    if (stream.fail()) {
+        cerr << "stream failed" << endl;
+        return false;
+    } else {
+        return true;
+    }
+}
+
+bool ExifTags::read_from_jpeg(istream &stream, ExifTags &output) {
+    unsigned char buf[2];
+    stream.read((char*)buf, 2);
+    if (!(buf[0] == 0xff && buf[1] == 0xd8)) {
+        return false;
+    }
+
+    while (stream.good()) {
+        stream.read((char*)buf, 2);
+        if (buf[0] != 0xff) {
+            cerr << "not a tag" << endl;
+            return false;
+        }
+        if (buf[1] == 0xff || buf[1] == 0x00) {
+            cerr << "not a marker" << endl;
+            return false;
+        }
+
+        if (buf[1] == 0xd9) {
+            break;
+        }
+
+        if (buf[1] == 0xda || (buf[1] >= 0xd0 && buf[1] <= 0xd7)) {
+            // SOS marker or segment without length
+            unsigned char c = 0;
+            while(stream.good()) {
+                c = stream.get();
+                if (c == 0xff) {
+                    c = stream.get();
+                    if (c != 0x00) {
+                        // end of section
+                        stream.seekg(-2, std::ios_base::cur);
+                        break;
+                    }
+                }
+            }
+        } else {
+            // read size, seek to next marker
+
+            uint16_t s = read_be_uint16(stream);
+
+            if (s < 2) {
+                cerr << "invalid size" << endl;
+                return false;
+            }
+
+            if (buf[1] == 0xe1) {
+                return parse_exif(stream, s, output);
+            }
+
+            stream.seekg(s - 2, std::ios_base::cur);
+        }
+    }
+
+    if (stream.fail()) {
+        cerr << "stream failed" << endl;
+        return false;
+    } else {
+        return true;
+    }
+}
+
+bool ExifTags::get_exposure_time(double &output) {
+    if (exposure_time.p != 0 or exposure_time.q != 0) {
+        output = (double)exposure_time.p / (double)exposure_time.q;
+        return true;
+    } else if (shutter_speed_value.p != 0 or shutter_speed_value.q != 0) {
+        output = 1.0 / pow(2.0, (double)shutter_speed_value.p / (double)shutter_speed_value.q);
+        return true;
+    } else {
+        return false;
+    }
+}

exif.h

@@ -0,0 +1,36 @@
+#include <stdint.h>
+#include <iostream>
+
+// unsigned rational number: p/q
+struct urational_t {
+    urational_t();
+    urational_t(uint32_t p, uint32_t q);
+
+    uint32_t p; // numerator
+    uint32_t q; // denominator
+};
+
+// signed rational number: p/q
+struct srational_t {
+    srational_t();
+    srational_t(int32_t p, int32_t q);
+
+    int32_t p; // numerator
+    int32_t q; // denominator
+};
+
+struct ExifTags {
+    ExifTags();
+
+    static bool read_from_jpeg(std::istream &stream, ExifTags &output);
+
+    // Returns true => output reference is set to exposure time or shutter speed in seconds.
+    // Returns false => do not set output because there is no valid value.
+    bool get_exposure_time(double &output);
+
+    bool ifd0_parsed_ok;
+    bool exififd_parsed_ok;
+    uint16_t exififd_offset;
+    urational_t exposure_time;
+    srational_t shutter_speed_value;
+};

haader.cpp

@@ -1,13 +1,247 @@
 #include <iostream>
 #include <fstream>
+#include <algorithm>
 #include <stdlib.h>
 #include <time.h>
+#include <math.h>
 
 #include "haader.h"
 
 using namespace std;
 using namespace haader;
 
+
+int ResponseFunction::lookup(double value) const {
+    // + 0.5 for rounding to nearest integer
+    int bin = (double)m_lookup_table.size() * (value - m_min_value) * m_recip_range + 0.5;
+
+    if (bin < 0) {
+        return 0;
+    } else if ((unsigned int)bin >= m_lookup_table.size()) {
+        return 255;
+    } else {
+        return m_lookup_table[bin];
+    }
+}
+
+double InverseResponseFunction::lookup(int index) const {
+    return m_lookup_table[index];
+}
+
+int InverseResponseFunction::inverse_lookup(double value) const {
+    if (value <= m_lookup_table[0]) {
+        return 0;
+    }
+    if (value >= m_lookup_table[255]) {
+        return 255;
+    }
+
+    for (unsigned int i = 1; i < 256; i++) {
+        if (value > m_lookup_table[i - 1] && value <= m_lookup_table[i]) {
+            double diff_a = value - m_lookup_table[i - 1];
+            double diff_b = m_lookup_table[i] - value;
+            if (diff_a < diff_b) {
+                return i - 1;
+            } else {
+                return i;
+            }
+        }
+    }
+
+    // should not happen
+    return 0;
+}
+
+ResponseFunction InverseResponseFunction::to_response_function(int bins) {
+    if (!is_monotonic()) {
+        cerr << "InverseResponseFunction::to_response_function: not monotonic." << endl;
+        return ResponseFunction();
+    }
+
+    double min_value = m_lookup_table[0];
+    double max_value = m_lookup_table[m_lookup_table.size() - 1];
+
+    ResponseFunction rf;
+    rf.m_min_value = min_value;
+    rf.m_max_value = max_value;
+    rf.m_recip_range = 1.0 / (max_value - min_value);
+    rf.m_lookup_table.resize(bins);
+
+
+    //TODO reduce quadratic complexity to linear complexity
+    for (int i = 0; i < bins; i++) {
+        double y = min_value + ((max_value - min_value) * i) / (double)(bins - 1.0);
+        rf.m_lookup_table[i] = inverse_lookup(y);
+    }
+
+    return rf;
+}
+
+bool InverseResponseFunction::repair() {
+
+    // Fill missing values
+    fill_zeros();
+
+    if (is_monotonic()) {
+        return true;
+    }
+
+    for (unsigned int i = 0; i < 128; i++) {
+        cout << "low pass filter " << i + 1 << endl;
+        low_pass_filter();
+        if (is_monotonic()) {
+            // repair success
+
+            // offset values so that m_lookup_table[128] is 0.0 again
+            double offset = -m_lookup_table[128];
+            for (unsigned int i = 0; i < 256; i++) {
+                m_lookup_table[i] += offset;
+            }
+
+            return true;
+        }
+    }
+
+    // repair failed
+    return false;
+}
+
+void InverseResponseFunction::fill_zeros() {
+    //TODO test correctness
+
+    // fill at front
+    if (m_lookup_table[0] == 0) {
+        for (unsigned int i = 0; i < 128; i++) {
+            if (m_lookup_table[i] != 0.0) {
+                for (unsigned int k = 0; k <= i; k++) {
+                    m_lookup_table[k] = m_lookup_table[i];
+                }
+                break;
+            }
+        }
+    }
+
+    // fill at back
+    if (m_lookup_table[255] == 0) {
+        for (unsigned int i = 255; i > 128; i--) {
+            if (m_lookup_table[i] != 0.0) {
+                for (unsigned int k = 255; k >= i; k--) {
+                    m_lookup_table[k] = m_lookup_table[i];
+                }
+                break;
+            }
+        }
+    }
+
+    {
+        bool inside_zeros = false;
+        int start_index = 0;
+        for (unsigned int i = 1; i < 256; i++) {
+            if (m_lookup_table[i] == 0.0 && !inside_zeros && i != 128) {
+                start_index = i;
+                inside_zeros = true;
+            } else if ((m_lookup_table[i] != 0.0 || i == 128) && inside_zeros) {
+                inside_zeros = false;
+                double a_val = m_lookup_table[start_index - 1];
+                double b_val = m_lookup_table[i];
+                for (unsigned int k = start_index; k < i; k++) {
+                    double t = (double)(k - start_index + 1) / (double)(i - start_index + 1);
+                    m_lookup_table[k] = (1.0 - t) * a_val + t * b_val;
+                }
+            }
+        }
+    }
+}
+
+void InverseResponseFunction::low_pass_filter() {
+    vector<double> newvals;
+    newvals.resize(256);
+
+    newvals[0] = m_lookup_table[0];
+    newvals[255] = m_lookup_table[255];
+
+    for (unsigned int i = 1; i < 255; i++) {
+        newvals[i] = 0.25 * m_lookup_table[i - 1] +  0.5 * m_lookup_table[i] + 0.25 * m_lookup_table[i + 1];
+    }
+
+    for (unsigned int i = 0; i < 256; i++) {
+        m_lookup_table[i] = newvals[i];
+    }
+}
+
+bool InverseResponseFunction::is_monotonic() {
+    for (unsigned int i = 1; i < 256; i++) {
+        if (m_lookup_table[i - 1] > m_lookup_table[i]) {
+            return false;
+        }
+    }
+    return true;
+}
+
+
+HdrImage::HdrImage():
+    m_width(0),
+    m_height(0),
+    m_components(0)
+{
+}
+
+Image HdrImage::get_log_image() {
+    if (m_image_data.size() == 0) {
+        return Image();
+    }
+
+    Image output;
+    output.m_image_data.resize(m_image_data.size());
+    output.m_width = m_width;
+    output.m_height = m_height;
+    output.m_components = m_components;
+
+    double min_value = m_image_data[0];
+    double max_value = m_image_data[0];
+
+    unsigned int size = m_image_data.size();
+    for (unsigned int i = 1; i < size; i++) {
+        double val = m_image_data[i];
+        if (val < min_value) {
+            min_value = val;
+        }
+        if (val > max_value) {
+            max_value = val;
+        }
+    }
+
+    double scale = 256.0 / (max_value - min_value);
+
+    for (unsigned int i = 0; i < size; i++) {
+        output.m_image_data[i] = (int)((m_image_data[i] - min_value) * scale);
+    }
+
+    return output;
+}
+
+Image HdrImage::expose(double exposure_time, const ResponseFunction &rf, double compression) {
+    if (m_image_data.size() == 0) {
+        return Image();
+    }
+
+    Image output;
+    output.m_image_data.resize(m_image_data.size());
+    output.m_width = m_width;
+    output.m_height = m_height;
+    output.m_components = m_components;
+
+    double log_exposure = log(exposure_time);
+
+    unsigned int size = m_image_data.size();
+    for (unsigned int i = 0; i < size; i++) {
+        output.m_image_data[i] = rf.lookup((m_image_data[i] + log_exposure) * compression);
+    }
+
+    return output;
+}
+
+
 HdrImageStack::HdrImageStack() {
 }
 
@@ -30,11 +264,20 @@ bool HdrImageStack::read_from_files(char * const* file_paths, int number_of_path
                      << endl;
                 return false;
             }
+            if (m_images[i].get_exposure_time() == 0.0) {
+                cerr << "HdrImageStack::read_from_files: Could not get exposure time (\""
+                     << file_paths[i]
+                     << "\")"
+                     << endl;
+                return false;
+            }
         } else {
             return false;
         }
     }
 
+    sort_by_exposure();
+
     return true;
 }
 
@@ -91,7 +334,7 @@ bool HdrImageStack::get_average_image(Image &output) {
         }
     }
 
-    // TODO test this with an image that is prime numbered dimensions
+    // TODO test this with an image that has prime numbered dimensions
     unsigned int rest = bytes_size % buf_size;
     for (unsigned int i = bytes_size - rest; i < bytes_size; i++) {
         int v = 0;
@@ -104,34 +347,86 @@ bool HdrImageStack::get_average_image(Image &output) {
     return true;
 }
 
-bool HdrImageStack::samples_to_csv(const char *file_path) {
-    if (m_images.size() == 0) {
-        return false;
-    }
-
+bool HdrImageStack::get_inverse_response_function(InverseResponseFunction &output, unsigned int num_samples) {
     srand(time(0));
-    ofstream output(file_path);
-
-    if (output.is_open()) {
-        unsigned int images_size = m_images.size();
-        unsigned int c = m_images[0].m_components;
-        unsigned int width = m_images[0].m_width;
-        unsigned int height = m_images[0].m_height;
-
-        for (unsigned int i = 0; i < 100; i++) {
-            unsigned int x = rand() % width;
-            unsigned int y = rand() % height;
-            unsigned int index = (y * width + x) * c;
-            for (unsigned int k = 0; k < images_size; k++) {
-                int val = m_images[k].m_image_data[index];
-                output << val;
-                if (k + 1 < images_size) {
-                    output << ',';
-                } else {
-                    output << '\n';
+
+    unsigned int images_size = m_images.size();
+    unsigned int c = m_images[0].m_components;
+    unsigned int width = m_images[0].m_width;
+    unsigned int height = m_images[0].m_height;
+
+    vector<vector<double> > resp;
+    resp.resize(256);
+
+    for (unsigned int i = 0; i < num_samples; i++) {
+        vector<int> values;
+        vector<double> exposure_times;
+        unsigned int x = rand() % width;
+        unsigned int y = rand() % height;
+        unsigned int index = (y * width + x) * c;
+        for (unsigned int k = 0; k < images_size; k++) {
+            int val = m_images[k].m_image_data[index];
+            if (val == 255) {
+                break;
+            }
+
+            values.push_back(val);
+            exposure_times.push_back(log(m_images[k].get_exposure_time()));
+        }
+
+        if (values.size() > 0 && values[0] < 128 && values[values.size() - 1] > 128) {
+            double offset = 0.0;
+            for (unsigned int k = 0; k < values.size(); k++) {
+                if (values[k] == 128) {
+                    offset = -exposure_times[k];
+                    break;
+                } else if (values[k] > 128) {
+                    double delta_val = values[k] - values[k -1];
+                    double t = (double)(128 - values[k - 1]) / delta_val;
+                    offset = -((1.0 - t) * exposure_times[k - 1] + t * exposure_times[k]);
+                    break;
                 }
             }
+
+            for (unsigned int k = 0; k < values.size(); k++) {
+                resp[values[k]].push_back(exposure_times[k] + offset);
+            }
+        }
+    }
+
+    output.m_lookup_table.resize(256);
+
+    for (unsigned int i = 0; i < 256; i++) {
+        double avg = 0.0;
+        for (unsigned int k = 0; k < resp[i].size(); k++) {
+            avg += resp[i][k];
+        }
+        if (resp[i].size() > 0) {
+            output.m_lookup_table[i] = avg / (double)resp[i].size();
+        } else {
+            output.m_lookup_table[i] = 0.0;
+        }
+    }
+
+    for (unsigned int i = 0; i < 256; i++) {
+        double avg = 0.0;
+        for (unsigned int k = 0; k < resp[i].size(); k++) {
+            avg += resp[i][k];
+        }
+        if (resp[i].size() > 0) {
+            output.m_lookup_table[i] = avg / (double)resp[i].size();
+        } else {
+            output.m_lookup_table[i] = 0.0;
         }
+    }
+
+    if (output.repair()) {
+        cout << "Inverse response function" << endl;
+        cout << endl;
+        for (unsigned int i = 0; i < 256; i++) {
+            cout << i << "," << output.m_lookup_table[i] << endl;
+        }
+        cout << endl;
 
         return true;
     } else {
@@ -139,8 +434,50 @@ bool HdrImageStack::samples_to_csv(const char *file_path) {
     }
 }
 
-bool HdrImageStack::recover_response_function(ResponseFunction &respFunc) {
-    //TODO implement
+HdrImage HdrImageStack::get_hdr_image(const InverseResponseFunction &invRespFunc) {
+    if (m_images.size() == 0) {
+        return HdrImage();
+    }
+
+    HdrImage output;
 
-    return false;
+    output.m_image_data.resize(m_images[0].m_image_data.size());
+    output.m_width = m_images[0].m_width;
+    output.m_height = m_images[0].m_height;
+    output.m_components = m_images[0].m_components;
+
+    unsigned int bytes_size = m_images[0].m_image_data.size();
+    unsigned int images_size = m_images.size();
+
+    // precompute log of exposure times
+    vector<double> log_exposures;
+    log_exposures.resize(images_size);
+    for (unsigned int k = 0; k < images_size; k++) {
+        log_exposures[k] = log(m_images[k].get_exposure_time());
+    }
+
+    for (unsigned int i = 0; i < bytes_size; i++) {
+
+        double scale_sum = 0.0;
+        double output_val = 0.0;
+
+        for (unsigned int k = 0; k < images_size; k++) {
+            int val = m_images[k].m_image_data[i];
+            double scale = 1.0 - fabs((double)(128.0 - val) * (1.0 / 128.0));
+            output_val += scale * (invRespFunc.lookup(val) - log_exposures[k]);
+            scale_sum += scale;
+        }
+
+        output.m_image_data[i] = output_val / scale_sum;
+    }
+
+    return output;
+}
+
+static bool compare_image_by_exposure(const Image &a, const Image &b) {
+    return a.get_exposure_time() < b.get_exposure_time();
+}
+
+void HdrImageStack::sort_by_exposure() {
+    std::sort(m_images.begin(), m_images.end(), compare_image_by_exposure);
 }

haader.h

@@ -7,28 +7,121 @@
 
 namespace haader {
     class ResponseFunction {
+        friend class InverseResponseFunction;
+        public:
+            // Returns a pixel value (0-255) given the logarithm of a luminance value.
+            int lookup(double value) const;
+
+        private:
+            // minimum log of luminance value
+            double m_min_value;
+
+            // maximum log of luminance value
+            double m_max_value;
+
+            // Precomputed value:
+            // 1.0 / (m_max_value - m_min_value)
+            double m_recip_range;
+
+            // Response function is implemented as a lookup table of with 256 values
+            std::vector<int> m_lookup_table;
+    };
+
+    class InverseResponseFunction {
+        friend class HdrImageStack;
+        public:
+
+            // Return logarithm of luminance given a pixel value (0-255)
+            double lookup(int index) const;
+
+            // The inverse of the lookup function:
+            //   y == lookup(x)  ->  x == inverse_lookup(y)
+            // O(n) complexity, proportional to size of lookup table.
+            int inverse_lookup(double value) const;
+
+            // Construct a ResponseFunction object that allows efficient O(1) inverse_lookup.
+            // bins: Size of the lookup table for the response function
+            ResponseFunction to_response_function(int bins=1024);
+
         private:
+
+            // Returns true if a monotonically increasing lookup table without missing values
+            // can be constructed.
+            bool repair();
+
+            // Remove missing values by linear interpolation of neighboring values
+            void fill_zeros();
+
+            // Apply binomial low pass filter
+            void low_pass_filter();
+
+            // Returns true if values are monotonically increasing
+            bool is_monotonic();
+
+            // Inverse response function is implemented as a lookup table with 256 values
             std::vector<double> m_lookup_table;
     };
 
-    class HdrImageStack {
-        friend class ResponseFunction;
+    // High Dynamic Range Image.
+    // Stores the logarithm of the scene luminance for each pixel/channel in a high resolution.
+    class HdrImage {
+        friend class HdrImageStack;
+
+        public:
+            HdrImage();
+
+            // Construct low dynamic range image with the logarithm of the scene luminance.
+            Image get_log_image();
+
+            // Construct low dynamic range image by a simulated exposure with a given response function.
+            // exposure_time: Exposure time in seconds
+            // compression:
+            //   A value > 1.0 will compress the luminace values and give more contrast.
+            //   A value < 1.0 will reduce the contrast but lower the risk of
+            //     under-/overexposure in parts of the image.
+            Image expose(double exposure_time, const ResponseFunction &rf, double compression=1.0);
+
+        private:
+            // Pixel values (logarithm of scene luminance). The channels are interleaved (RGBRGBRGB...).
+            std::vector<double> m_image_data;
+
+            // Width in pixels
+            unsigned int m_width;
+
+            // Height in pixels
+            unsigned int m_height;
+
+            // Number of channels (usually 3 for RGB images)
+            unsigned int m_components;
+    };
 
+    // A stack of images/photographs of the same scene, with the same resolution but different exposure times.
+    class HdrImageStack {
         public:
             HdrImageStack();
 
+            // Read images from the given array of file paths
             bool read_from_files(char * const* file_paths, int number_of_paths);
 
             bool get_average_image(Image &output);
             bool get_average_image_slow(Image &output);
 
-            bool samples_to_csv(const char *file_path);
+            // Approximate the inverse response function from samples of random pixel locations.
+            // num_samples: Number of samples
+            // Returns true if a well-formed inverse response function could be found.
+            bool get_inverse_response_function(InverseResponseFunction &output, unsigned int num_samples);
 
-            bool recover_response_function(ResponseFunction &respFunc);
+            // Return a HDR-image given an inverse response function
+            HdrImage get_hdr_image(const InverseResponseFunction &invRespFunc);
 
         private:
+            // Sort the images by exposure time
+            void sort_by_exposure();
+
+            // The images
             std::vector<Image> m_images;
     };
+
 }
 
 #endif // HAADER_ILUCF4Z0

image.cpp

@@ -14,6 +14,7 @@
 #include "stb_image.h"
 
 #include "image.h"
+#include "exif.h"
 
 using namespace std;
 using namespace haader;
@@ -22,11 +23,12 @@ using namespace haader;
 Image::Image():
     m_width(0),
     m_height(0),
-    m_components(0)
+    m_components(0),
+    m_exposure_time(0.0)
 {
 }
 
-void Image::to_string(string &s) {
+void Image::to_string(string &s) const {
     stringstream ss;
     ss << "Image(";
     ss << m_width << 'x' << m_height << 'x' << m_components << ")";
@@ -47,6 +49,19 @@ bool Image::read_from_file(const char *file_path) {
         size_t size = (size_t)x * (size_t)y * (size_t)components;
         m_image_data.assign(data, data + size);
 
+        // try to read EXIF data
+        {
+            ExifTags tags;
+            ifstream stream;
+            stream.open(file_path);
+
+            if (stream.is_open() and ExifTags::read_from_jpeg(stream, tags)) {
+                if (tags.get_exposure_time(m_exposure_time)) {
+                    cout << "Exposure time: " << m_exposure_time << " sec" << endl;
+                }
+            }
+        }
+
         return true;
     } else {
         cerr << "Can not read \"" << file_path << "\" as an image." << endl;
@@ -218,7 +233,7 @@ bool Image::read_from_jpeg_file(const char *file_path) {
     }
 }
 
-bool Image::save_as_ppm_file(const char *file_path, bool ascii) {
+bool Image::save_as_ppm_file(const char *file_path, bool ascii) const {
     if (ascii) {
         // ASCII
 
@@ -256,20 +271,24 @@ bool Image::save_as_ppm_file(const char *file_path, bool ascii) {
     }
 }
 
-unsigned int Image::get_width() {
+unsigned int Image::get_width() const {
     return m_width;
 }
 
-unsigned int Image::get_height() {
+unsigned int Image::get_height() const {
     return m_height;
 }
 
-unsigned int Image::get_components() {
+unsigned int Image::get_components() const {
     return m_components;
 }
 
-bool Image::has_equal_dimensions(const Image &other) {
+bool Image::has_equal_dimensions(const Image &other) const {
     return (m_width == other.m_width &&
             m_height == other.m_height &&
             m_components == other.m_components);
 }
+
+double Image::get_exposure_time() const {
+    return m_exposure_time;
+}

image.h

@@ -6,28 +6,36 @@
 
 namespace haader {
 
+    // A low dynamic range image/photograph with an associated exposure time.
     class Image {
         friend class HdrImageStack;
+        friend class HdrImage;
 
         public:
             Image();
 
-            void to_string(std::string &s);
+            void to_string(std::string &s) const;
 
             bool read_from_file(const char *file_path);
             bool read_from_jpeg_file(const char *file_path);
-            bool save_as_ppm_file(const char *file_path, bool ascii=false);
+            bool save_as_ppm_file(const char *file_path, bool ascii=false) const;
 
-            unsigned int get_width();
-            unsigned int get_height();
-            unsigned int get_components();
-            bool has_equal_dimensions(const Image &other);
+            unsigned int get_width() const;
+            unsigned int get_height() const;
+            unsigned int get_components() const;
+            bool has_equal_dimensions(const Image &other) const;
+            double get_exposure_time() const;
 
         private:
             std::vector<unsigned char> m_image_data;
             unsigned int m_width;
             unsigned int m_height;
+
+            // Number of channels (usually 3 for RGB images)
             unsigned int m_components;
+
+            // Exposure time in seconds
+            double m_exposure_time;
     };
 
 }

main.cpp

@@ -1,32 +1,84 @@
 #include "haader.h"
 
 #include <time.h>
+#include <math.h>
 
 #include <iostream>
+#include <sstream>
+#include <iomanip>
 
 using namespace std;
 
 
+void print_usage(const char *argv_0) {
+    cout << "Construct an HDR-image from a series of photographs\n";
+    cout << "with different exposure times and then synthesize\n";
+    cout << "new images from it.\n\n";
+    cout << "Usage:\n";
+    cout << argv_0 << " <image_path_1> <image_path_2> ...\n\n";
+}
+
 int main(int argc, char *argv[])
 {
     if (argc >= 2) {
         haader::HdrImageStack stack;
 
-        if (stack.read_from_files(argv + 1, argc - 1)) {
-            haader::Image img;
+        // read image files.
+        bool ok = stack.read_from_files(argv + 1, argc - 1);
+
+        if (ok) {
+            clock_t start = clock();
 
+            // construct and save average image
             {
-                clock_t start = clock();
+                haader::Image img;
                 stack.get_average_image(img);
-                clock_t end = clock();
-                double seconds = (end - start) / (double)CLOCKS_PER_SEC;
-                cout << "secs: " << seconds << endl;
+                img.save_as_ppm_file("average.ppm");
             }
 
-            img.save_as_ppm_file("average.ppm");
-            stack.samples_to_csv("samples.csv");
-        }
+            // approximate inverse response function
+            haader::InverseResponseFunction irf;
+            bool irf_ok = stack.get_inverse_response_function(irf, 2000);
+
+            if (irf_ok) {
+                haader::HdrImage hdr_img;
+                hdr_img = stack.get_hdr_image(irf);
+
+                // save logarithmic image of scene luminance values
+                {
+                    haader::Image x = hdr_img.get_log_image();
+                    x.save_as_ppm_file("log_image.ppm");
+                }
 
+                // convert irf to response function
+                haader::ResponseFunction rf;
+                rf = irf.to_response_function(1024);
+
+                // create images by a simulated exposure of the HDR image with the response function
+                cout << "Expose..." << endl;
+                for (int i = 0; i < 20; i++) {
+                    // exposure time
+                    double t = pow(2.0, (double)i * (-12.0 / 20.0));
+
+                    // create image
+                    haader::Image x = hdr_img.expose(t, rf, 1.0);
+
+                    // save image
+                    stringstream ss;
+                    ss << "expose_" << setw(5) << setfill('0') << i << ".ppm";
+                    string path = ss.str();
+                    x.save_as_ppm_file(path.c_str());
+
+                    cout << "\nexposure time: " << t << " sec" << endl;
+                    cout << "save as " << path << endl;
+                }
+            }
+            clock_t end = clock();
+            double seconds = (end - start) / (double)CLOCKS_PER_SEC;
+            cout << "\nsecs: " << seconds << endl;
+        }
+    } else {
+        print_usage(argv[0]);
     }
 
     return 0;