processimage.cpp

/* -LICENSE-START-
 ** Copyright (c) 2018 Blackmagic Design
 **
 ** Permission is hereby granted, free of charge, to any person or organization
 ** obtaining a copy of the software and accompanying documentation covered by
 ** this license (the "Software") to use, reproduce, display, distribute,
 ** execute, and transmit the Software, and to prepare derivative works of the
 ** Software, and to permit third-parties to whom the Software is furnished to
 ** do so, all subject to the following:
 **
 ** The copyright notices in the Software and this entire statement, including
 ** the above license grant, this restriction and the following disclaimer,
 ** must be included in all copies of the Software, in whole or in part, and
 ** all derivative works of the Software, unless such copies or derivative
 ** works are solely in the form of machine-executable object code generated by
 ** a source language processor.
 **
 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 ** IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 ** FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 ** SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 ** FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 ** ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 ** DEALINGS IN THE SOFTWARE.
 ** -LICENSE-END-
 */

#include "ProcessImage.h"
#include "CintelRawImage.h"


ProcessImage::ProcessImage(const CintelRawImage& cri) :
		m_cri(cri),
		m_width(m_cri.width()),
		m_height(m_cri.height()),
		m_data(m_cri.getRawFrame())
{
	if (m_data.size() != m_width * m_height)
		throw std::logic_error("Can't process without valid image");

	if (m_width % 2 != 0 || m_height % 2 != 0)
		throw std::logic_error("Processing requires even width and height");
}

void ProcessImage::debayer()
{
	auto&	r = m_rgb[0];
	auto&	g = m_rgb[1];
	auto&	b = m_rgb[2];

	for (unsigned y = 0; y < m_height; y += 2)
	{
		for (unsigned x = 0; x < m_width; x += 2)
		{
			m_data[(y + 0) * m_width + x + 0];
		}
	}

	r.resize(m_data.size());
	g.resize(m_data.size());
	b.resize(m_data.size());

	// Using a single pass to extract the 3 color channels to minimise memory access time
	for (unsigned y = 0; y < m_height; y++)
	{
		for (unsigned x = 0; x < m_width; x++)
		{
			// Bayer pattern is:
			//    0  1  0  1
			//   +-----------
			// 0 |G0 R  G0 R
			// 1 |B  G1 B  G1
			// 0 |G0 R  G0 R
			// 1 |B  G1 B  G1
			const unsigned	bayerX = x % 2;
			const unsigned	bayerY = y % 2;

			if (bayerX == 0 && bayerY == 0)							// G0 sample
			{
				g[idx(x, y)] = m_data[idx(x, y)];						// straight copy for pure green sample
				r[idx(x, y)] = mean2(x, y, -1, 0, 1, 0);				// interpolate from 2 horizontal red neighbors
				b[idx(x, y)] = mean2(x, y, 0, -1, 0, 1);				// interpolate from 2 vertical blue neighbors
			}
			else if (bayerX == 1 && bayerY == 0)					// R sample
			{
				g[idx(x, y)] = mean4(x, y, {0,-1, -1,0, 1,0, 0,1});		// interpolate from 4 green neighbors
				r[idx(x, y)] = m_data[idx(x, y)];						// straight copy for pure red sample
				b[idx(x, y)] = mean4(x, y, {-1,-1, 1,-1, -1,1, 1,1});	// interpolate from 4 blue neighbors
			}
			else if (bayerX == 0 && bayerY == 1)					// B sample
			{
				g[idx(x, y)] = mean4(x, y, {0,-1, -1,0, 1,0, 0,1});		// interpolate from 4 green neighbors
				r[idx(x, y)] = mean4(x, y, {-1,-1, 1,-1, -1,1, 1,1});	// interpolate from 4 red neighbors
				b[idx(x, y)] = m_data[idx(x, y)];						// straight copy for pure blue sample
			}
			else													// G1 sample
			{
				g[idx(x, y)] = m_data[idx(x, y)];						// straight copy for pure green sample
				r[idx(x, y)] = mean2(x, y, 0, -1, 0, 1);				// interpolate from 2 vertical red neighbors
				b[idx(x, y)] = mean2(x, y, -1, 0, 1, 0);				// interpolate from 2 horizontal blue neighbors
			}
		}
	}
}

// Return the mean of the 2 neighboring pixels at the specified offsets
float ProcessImage::mean2(int xOrigin, int yOrigin, int xOffset1, int yOffset1, int xOffset2, int yOffset2)
{
	unsigned sum = 0;

	int x = xOrigin + xOffset1;
	int y = yOrigin + yOffset1;

	// If the first neighbor is out of bounds, the second will be in bounds
	if (x < 0 || x >= m_width || y < 0 || y >= m_height)
		return m_data[idx(xOrigin + xOffset2, yOrigin + yOffset2)];

	sum += m_data[idx(x, y)];

	x = xOrigin + xOffset2;
	y = yOrigin + yOffset2;

	// If the second neighbor is out of bounds, return the first neighbor
	if (x < 0 || x >= m_width || y < 0 || y >= m_height)
		return sum;

	sum += m_data[idx(x, y)];

	return float(sum / 2.0f);
}

// Return the mean of the 4 neighboring pixels at the specified offsets
float ProcessImage::mean4(int xOrigin, int yOrigin, std::initializer_list<int> offsets)
{
	unsigned sum = 0;
	unsigned count = 0;

	for (auto it = offsets.begin(); it != offsets.end(); )
	{
		const int x = xOrigin + *it++;
		const int y = yOrigin + *it++;

		if (x < 0 || x >= m_width || y < 0 || y >= m_height)
			continue;			// skip pixels outside image boundary

		sum += m_data[idx(x, y)];
		count++;
	}

	float mean = float(sum) / count;
	return mean;
}

unsigned ProcessImage::idx(unsigned x, unsigned y)
{
	return y * m_width + x;
}

void ProcessImage::linearMask()
{
	if (! (m_rgb[0].size() == m_rgb[1].size() && m_rgb[0].size() == m_rgb[2].size()))
		throw std::logic_error("RGB channels must be equal size");

	const auto linearMask = m_cri.linearMask();

	// Shortcut if linearMask is unity
	if (linearMask == std::array<float, 9> {{ 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }})
		return;

	for (size_t i = 0; i < m_rgb[0].size(); i++)
	{
		const auto r = m_rgb[0][i];
		const auto g = m_rgb[1][i];
		const auto b = m_rgb[2][i];
		m_rgb[0][i] = r * linearMask[0] + g * linearMask[3] + b * linearMask[6];
		m_rgb[1][i] = r * linearMask[1] + g * linearMask[4] + b * linearMask[7];
		m_rgb[2][i] = r * linearMask[2] + g * linearMask[5] + b * linearMask[8];
	}
}

namespace
{
	float removeOffset(float pixelDN)
	{
		// The Cintel scanner has a black offset of 20 that must be removed before applying the log function
		const float		kSensorFullRange = 4095.0f;
		const float		kCintelWhiteOffset = 4095.0f / kSensorFullRange;
		const float		kCintelBlackOffset = 20.0f / kSensorFullRange;

		return (pixelDN / kSensorFullRange - kCintelBlackOffset) / (kCintelWhiteOffset - kCintelBlackOffset);
	}

	float cintelLog(float pixel, bool usePrintLog)
	{
		// Not technically a 'gamma' in the normal sense but it's used the same way - this is from the Cineon
		// standard of 0.002 density units per stored DN (so at 10 bits it's 0.002 * 1023 = 2.046).
		// Note that earlier standards incorrectly state this as 2.048.
		const float		kCineonGamma	= 2.046f;
		const float		kFilmGammaNeg	= 1.0f / kCineonGamma;

		const float		kDisplayGamma	= 2.200f;
		const float		kFilmGammaPrint	= 1.0f / kDisplayGamma;

		const float		gamma = usePrintLog ? kFilmGammaPrint : kFilmGammaNeg;

		float a = 0.0f;
		float b = 1.0f;

		if (usePrintLog)
		{
			a = 1.0f / (pow(10.0f, 1.0f / gamma) - 1.0f);
			b = 1.0f -gamma * log10f(1.0f  + a);
		}

		// One sided clamp to avoid negative or zero values before calling log function
		pixel = std::max(pixel, std::numeric_limits<float>::min());

		// Apply log curve
		float logValue = gamma * log10f(pixel + a) + b;

		// Clamp to a range a little larger than pixel range to provide headroom for other processing i.e. [-0.2, 1.2]
		return std::max(-0.2f, std::min(logValue, 1.2f));
	}
}

void ProcessImage::logMask()
{
	if (! (m_rgb[0].size() == m_rgb[1].size() && m_rgb[0].size() == m_rgb[2].size()))
		throw std::logic_error("RGB channels must be equal size");

	bool invert;
	bool printLog;
	switch (m_cri.filmType())
	{
		case CintelRawImage::FilmType::Positive:
			invert = true;
			printLog = true;
			break;
		case CintelRawImage::FilmType::Negative:
			invert = false;
			printLog = false;
			break;
		case CintelRawImage::FilmType::InterPositive:
			invert = true;
			printLog = false;
			break;
		case CintelRawImage::FilmType::InterNegative:
			invert = false;
			printLog = false;
			break;
		default:
			throw std::logic_error("Unsupported FilmType " + std::to_string(unsigned(m_cri.filmType())));
	}

	const auto logMask = m_cri.logMask();

	for (size_t i = 0; i < m_rgb[0].size(); i++)
	{
		// Remove the black offset and scale to full pixel range
		const float r = removeOffset(m_rgb[0][i]);
		const float g = removeOffset(m_rgb[1][i]);
		const float b = removeOffset(m_rgb[2][i]);

		// Apply log curve and invert (log mask is always applied on inverted data)
		float rInv = 1.0f - cintelLog(r, printLog);
		float gInv = 1.0f - cintelLog(g, printLog);
		float bInv = 1.0f - cintelLog(b, printLog);

		float rMasked = rInv * logMask[0] + gInv * logMask[3] + bInv * logMask[6];
		float gMasked = rInv * logMask[1] + gInv * logMask[4] + bInv * logMask[7];
		float bMasked = rInv * logMask[2] + gInv * logMask[5] + bInv * logMask[8];

		if (invert)
		{
			// Undo inversion for prints
			m_rgb[0][i] = 1.0f - rMasked;
			m_rgb[1][i] = 1.0f - gMasked;
			m_rgb[2][i] = 1.0f - bMasked;
		}
		else
		{
			// Leave negative inverted
			m_rgb[0][i] = rMasked;
			m_rgb[1][i] = gMasked;
			m_rgb[2][i] = bMasked;
		}
	}
}

void ProcessImage::gains()
{
	const auto gains = m_cri.gains();

	// Shortcut if gains are all unity
	if (gains == std::array<float, 3> {{ 1.0f, 1.0f, 1.0f }})
		return;

	for (size_t i = 0; i < m_rgb[0].size(); i++)
	{
		m_rgb[0][i] = m_rgb[0][i] * gains[0];
		m_rgb[1][i] = m_rgb[1][i] * gains[1];
		m_rgb[2][i] = m_rgb[2][i] * gains[2];
	}
}

void ProcessImage::lifts()
{
	const auto lifts = m_cri.lifts();

	// Shortcut if lifts are all zero
	if (lifts == std::array<float, 3> {{ 0.0f, 0.0f, 0.0f }})
		return;

	for (size_t i = 0; i < m_rgb[0].size(); i++)
	{
		m_rgb[0][i] = m_rgb[0][i] + lifts[0];
		m_rgb[1][i] = m_rgb[1][i] + lifts[1];
		m_rgb[2][i] = m_rgb[2][i] + lifts[2];
	}
}

void ProcessImage::applyStabilityOffsetsAndFlips()
{
	const auto flipHorizontal = m_cri.flipHorizontal();
	const auto flipVertical = m_cri.flipVertical();
	const auto offsetDetectedHorizontal = m_cri.offsetDetectedHorizontal();
	const auto offsetDetectedVertical = m_cri.offsetDetectedVertical();

	if (! flipHorizontal && ! flipVertical)
		return;

	// For each color
	for (int i = 0; i < m_rgb.size(); i++)
	{
		const auto&		color = m_rgb[i];

		// To apply an offset or flip we use a temporary copy of the image data for this color
		std::vector<float>	flipped;
		flipped.reserve(color.size());

		for (int y = 0; y < m_height; y++)
		{
			for (int x = 0; x < m_width; x++)
			{
				int row = flipVertical ? (m_height - 1 - y) : y;
				int col = flipHorizontal ? (m_width - 1 - x) : x;

				row += offsetDetectedVertical;
				col += offsetDetectedHorizontal;

				if (row < 0 || row >= m_height || col < 0 || col >= m_width)
					flipped.emplace_back(0.0f);									// set out-of-bounds pixels to black
				else
					flipped.emplace_back(color[row * m_width + col]);
			}
		}

		// Copy flipped image back
		m_rgb[i] = flipped;
	}
}