Their intention is to solve:
$$ \arg \min_{ J \left( x, y \right) } \sum_{x, y} \left\| J \left( x, y \right) - 255 \right\|_{2}^{2} + \lambda \sum_{x, y} \left\| \nabla J \left( x, y \right) - \nabla I \left( x, y \right) \right\|_{2}^{2} $$
Where $ I \left( x, y \right) $ is the input image and $ J \left( x, y \right) $ is the output image.
They state that usually the solution is constant locally and can be approximated by the map of Gain and Offset.
They do not say the method they use to calculate it but you can use few ideas:
- Solve the above optimization problem for a downsampled version of $ I \left( x, y \right) $. From the downsampled solution of $ J \left( x, y \right) $ infer the downsampled map of the Gain and Offset. Then upsample the maps and apply them on the full size $ I \left( x, y \right) $.
- Just like the Guided Filter approximate a local solution of the Weighted Least Squares filter, reformulate it for solution of the above problem (You'll see they are very similar).
Solving the Optimization Problem
As written in Dropbox's Blog Post - Fast Document Rectification and Enhancement the above is a Poisson Equation which has many dedicated solvers.
To give you a glimpse of straight forward way solve it we can vectorize the image into a vector. Let's say the input image is $ y $ and the output image is $ x $ (Vectors) then the problem can be rewritten:
$$ \hat{x} = \arg \min_{x} \frac{1}{2} {\left\| x - 255 \cdot \boldsymbol{1} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{h} \left( x - y \right) \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{v} \left( x - y \right) \right\|}_{2}^{2} $$
Where $ {D}_{h} $ is the horizontal Derivative Operator, $ {D}_{v} $ is the vertical Derivative Operator and $ \boldsymbol{1} $ is vector of ones.
Then the solution is given by:
$$ \hat{x} = { \left( I + \lambda {D}_{h}^{T} {D}_{h} + \lambda {D}_{v}^{T} {D}_{v} \right) }^{-1} \left( \lambda {D}_{h}^{T} {D}_{h} y + \lambda {D}_{v}^{T} {D}_{v} y + 255 \cdot \boldsymbol{1} \right) $$
Pay attention that if the image $ I \left( x, y \right) \in \mathbb{R}^{m \times n} $ then $ D \in \mathbb{R}^{m n \times m n} $.
Hence from memory and computation this is hard to do for large images unless the solver takes advantage of the matrix $ D $ structure (Sparse) and that's what Poisson Solvers are all about (Be more accurate $ {D}^{T} D $ has a structure to take advantage of without ever formally creating it).
For very small images you can play with this.
For large you can do follwoing:
- Use Direct Solvers
For example just using the \
operator in MATLAB.
- Use Iterative Solvers
For example the Conjugate Gradient. Nice property of the Conjugate Gradient is it doesn't require the matrices themselves but the products of the matrices with the image which can be calculated by convolution.
So basically the above is about solving Linear Equation with special structure. You may read about solving those kind of equations at the course The Poisson Equation in Image & Shape Processing.
MATLAB Code
I added a MATLAB Code to solve the problem of whitening the image as described in the Blog Post.
Here is the result for different values of $ \lambda $:

As expected, one could see that the higher the regularization the better the contrast.
The full MATLAB code which implements the solution is at my Signal Processing StackExchange Question 50329 - GitHub Repository.
The code uses direct solver. In practice one should use Preconditioned Conjugate Gradient solver which only requires the result of the operator on the image which can be done using Convolution in much faster manner than Matrix Multiplication (At least for "Denser" operators).
As reference image I used Certificate of Arrival for Berta Werner.