1 #pragma once
2 /*
3 The EdgePreservingDecomposition files contain standard C++ (standard except the
4 first line) code for creating and, to a limited extent (create your own uses!),
5 messing with multi scale edge preserving decompositions of a 32 bit single
6 channel image. As a byproduct it contains a lot of linear algebra which can be
7 useful for optimization problems that you want to solve in rectangles on
8 rectangular grids.
9 
10 Anyway. Basically, this is an implementation of what's presented in the
11 following papers:
12 - Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation
13 - An Iterative Solution Method for Linear Systems of Which the Coefficient
14   Matrix is a Symmetric M-Matrix
15 - Color correction for tone mapping
16 - Wikipedia, the free encyclopedia
17 
18 First one is most of what matters, next two are details, last everything else.
19 I did a few things differently, especially:
20 - Reformulated the minimization with finite elements instead of finite
21   differences. This results in better conditioning, slightly better accuracy
22   (less artifacts), the possibility of a better picked edge stopping function,
23   but more memory consumption.
24 - A single rotationally invariant edge stopping function is used instead of two
25   non-invariant ones.
26 - Incomplete Cholseky factorization instead of Szeliski's LAHBF. Slower, but
27   not subject to any patents.
28 - For tone mapping, original images are decomposed instead of their logarithms,
29   and just one decomposition is made;
30 - I find that this way works plenty good (theirs isn't better or worse... just
31   different) and is simpler.
32 
33 Written by ben_pcc in Portland, Oregon, USA. Some history:
34 - Late April 2010, I develop interest in this stuff because photos of my
35   ceramics lack local contrast.
36 - Mid 2010, it works but is too slow to be useful.
37 - Fall 2010, various unsuccessful attempts at speeding up are tried.
38 - Early December 2010, I get off the path of least resistance and write a
39   matrix storage class with incomplete Cholesky decomposition.
40 - 31 December 2010, the FEM reformulation works very well.
41 - 1 January 2011, I'm cleaning up this file and readying it for initial release.
42 - 12 - 14 November 2011, further cleanup, improvements, bug fixes, integration
43   into Raw Therapee.
44 
45 It's likely that I'll take apart and rerelease contents of this file (in the
46 distant future) as most of it isn't edge preserving decomposition and rather
47 supporting material. SparseConjugateGradient alone is a workhorse I and a few
48 others have been exploiting for a few years.
49 
50 EdgePreservingDecomposition.h and EdgePreservingDecomposition.cpp are released
51 under the following licence:
52 - It's free.
53 - You may not incorporate this code as part of proprietary or commercial
54   software, but via freeware you may use its output for profit.
55 - You may modify and redistribute, but keep this big comment block intact and
56   not for profit in any way unless I give specific permission.
57 - If you're unsure about anything else, treat as public domain.
58 - Don't be a dick.
59 
60 My email address is my screen name followed by @yahoo.com. I'm also known as
61 ben_s or nonbasketless. Enjoy!
62 */
63 
64 
65 
66 #include <cmath>
67 #include <cstdio>
68 #include <cstdlib>
69 #include <cstring>
70 
71 #include "opthelper.h"
72 #include "noncopyable.h"
73 
74 //This is for solving big symmetric positive definite linear problems.
75 float *SparseConjugateGradient(void Ax(float *Product, float *x, void *Pass), float *b, int n, bool OkToModify_b = true, float *x = nullptr, float RMSResidual = 0.0f, void *Pass = nullptr, int MaximumIterates = 0, void Preconditioner(float *Product, float *x, void *Pass) = nullptr);
76 
77 //Storage and use class for symmetric matrices, the nonzero contents of which are confined to diagonals.
78 class MultiDiagonalSymmetricMatrix :
79     public rtengine::NonCopyable
80 {
81 public:
82     MultiDiagonalSymmetricMatrix(int Dimension, int NumberOfDiagonalsInLowerTriangle);
83     ~MultiDiagonalSymmetricMatrix();
84 
85     /* Storage of matrix data, and a function to create memory for Diagonals[index].
86     Here's the storage scheme, designed to be simple both on paper and in C++:
87 
88     Let's say you have some diagonal. The StartRows is the row on which, at the left edge of the matrix, the diagonal "starts",
89     and StartRows must strictly increase with its index. The main diagonal for example has start row 0, its subdiagonal has 1, etc.
90     Then, Diagonal[j] is the matrix entry on the diagonal at column j. For efficiency, you're expected to learn this and fill in
91     public Diagonals manually. Symmetric matrices are represented by this class, and all symmetry is handled internally, you
92     only every worry or think about the lower triangular (including main diagonal) part of the matrix.
93     */
94     float **Diagonals;
95     char *buffer;
96     float *DiagBuffer;
97     int *StartRows;
98     bool CreateDiagonal(int index, int StartRow);
99     int n, m;   //The matrix is n x n, with m diagonals on the lower triangle. Don't change these. They should be private but aren't for convenience.
DiagonalLength(int StartRow)100     inline int DiagonalLength(int StartRow)     //Gives number of elements in a diagonal.
101     {
102         return n - StartRow;
103     };
104 
105     //Not efficient, but you can use it if you're lazy, or for early tests. Returns false if the row + column falls on no loaded diagonal, true otherwise.
106     bool LazySetEntry(float value, int row, int column);
107 
108     //Calculates the matrix-vector product of the matrix represented by this class onto the vector x.
109     void VectorProduct(float *Product, float *x);
110 
111     //Given the start row, attempts to find the corresponding index, or -1 if the StartRow doesn't exist.
112     inline int FindIndex(int StartRow) __attribute__((always_inline));
113 
114     //This is the same as above, but designed to take this class as a pass through variable. By this way you can feed
115     //the meat of this class into an independent function, such as SparseConjugateGradient.
PassThroughVectorProduct(float * Product,float * x,void * Pass)116     static void PassThroughVectorProduct(float *Product, float *x, void *Pass)
117     {
118         (static_cast<MultiDiagonalSymmetricMatrix *>(Pass))->VectorProduct(Product, x);
119     };
120 
121     /* CreateIncompleteCholeskyFactorization creates another matrix which is an incomplete (or complete if MaxFillAbove is big enough)
122     LDLt factorization of this matrix. Storage is like this: the first diagonal is the diagonal matrix D and the remaining diagonals
123     describe all of L except its main diagonal, which is a bunch of ones. Read up on the LDLt Cholesky factorization for what all this means.
124     Note that VectorProduct is nonsense. More useful to you is CholeskyBackSolve which fills x, where LDLt x = b. */
125     bool CreateIncompleteCholeskyFactorization(int MaxFillAbove = 0);
126     void KillIncompleteCholeskyFactorization(void);
127     void CholeskyBackSolve(float *x, float *b);
128     MultiDiagonalSymmetricMatrix *IncompleteCholeskyFactorization;
129 
PassThroughCholeskyBackSolve(float * Product,float * x,void * Pass)130     static void PassThroughCholeskyBackSolve(float *Product, float *x, void *Pass)
131     {
132         (static_cast<MultiDiagonalSymmetricMatrix *>(Pass))->CholeskyBackSolve(Product, x);
133     };
134 
135 };
136 
137 class EdgePreservingDecomposition :
138     public rtengine::NonCopyable
139 {
140 public:
141     EdgePreservingDecomposition(int width, int height);
142     ~EdgePreservingDecomposition();
143 
144     //Create an edge preserving blur of Source. Will create and return, or fill into Blur if not NULL. In place not ok.
145     //If UseBlurForEdgeStop is true, supplied not NULL Blur is used to calculate the edge stopping function instead of Source.
146     float *CreateBlur(float *Source, float Scale, float EdgeStopping, int Iterates, float *Blur = nullptr, bool UseBlurForEdgeStop = false);
147 
148     //Iterates CreateBlur such that the smoothness term approaches a specific norm via iteratively reweighted least squares. In place not ok.
149     float *CreateIteratedBlur(float *Source, float Scale, float EdgeStopping, int Iterates, int Reweightings, float *Blur = nullptr);
150 
151     /*Lowers global contrast while preserving or boosting local contrast. Can fill into Compressed. The smaller Compression
152     the more compression is applied, with Compression = 1 giving no effect and above 1 the opposite effect. You can totally
153     use Compression = 1 and play with DetailBoost for some really sweet unsharp masking. If working on luma/grey, consider giving it a logarithm.
154     In place calculation to save memory (Source == Compressed) is totally ok. Reweightings > 0 invokes CreateIteratedBlur instead of CreateBlur. */
155     void CompressDynamicRange(float *Source, float Scale = 1.0f, float EdgeStopping = 1.4f, float CompressionExponent = 0.8f, float DetailBoost = 0.1f, int Iterates = 20, int Reweightings = 0);
156 
157 private:
158     MultiDiagonalSymmetricMatrix *A;    //The equations are simple enough to not mandate a matrix class, but fast solution NEEDS a complicated preconditioner.
159     int w, h, n;
160 
161     //Convenient access to the data in A.
162     float * RESTRICT a0, * RESTRICT a_1, * RESTRICT a_w, * RESTRICT a_w_1, * RESTRICT a_w1;
163 };
164 
165