1+
2+
3+
4+
5+ function count_saturated_N0f8 (img)
6+ sum ((s-> s == N0f8 (1 )). (img))
7+ end
8+
9+ function percentile_saturated_RGB (img)
10+ perc = 0.0
11+ perc += count_saturated_N0f8 ((s-> s. r). (img))
12+ perc += count_saturated_N0f8 ((s-> s. g). (img))
13+ perc += count_saturated_N0f8 ((s-> s. b). (img))
14+ perc /= 3
15+ perc /= * (size (img)... )
16+ return perc
17+ end
18+
19+
20+
21+ function estimateShutterTimes_RGB (
22+ imgs:: AbstractVector{<:AbstractMatrix} ,
23+ nominal:: Real = 1 / 100 ,
24+ )
25+ avgEk = []
26+ scalek = []
27+ for img in imgs
28+ avgE_ = (s-> (s. r + s. g + s. b)/ 3 )(mean (img))
29+ push! (avgEk, avgE_)
30+ # many saturated pixels mean stronger scaling on nominal / avgE
31+ push! (scalek, exp (percentile_saturated_RGB (img)) )
32+ end
33+ # avg energy in scene
34+ avgE = mean (avgEk)
35+
36+ ((nominal / avgE ) .* avgEk) .* scalek
37+ end
38+
39+ function percentile_saturated_gray (img)
40+ perc = 0.0
41+ perc += count_saturated_N0f8 ((s-> s. val). (img))
42+ perc /= * (size (img)... )
43+ return perc
44+ end
45+
46+ function estimateShutterTimes_gray (
47+ imgs:: AbstractVector{<:AbstractMatrix} ,
48+ nominal:: Real = 1 / 100 ,
49+ )
50+ avgEk = []
51+ scalek = []
52+ for img in imgs
53+ avgE_ = (s-> s. val)(mean (img))
54+ push! (avgEk, avgE_)
55+ # many saturated pixels mean stronger scaling on nominal / avgE
56+ push! (scalek, exp (percentile_saturated_gray (img)) )
57+ end
58+ # avg energy in scene
59+ avgE = mean (avgEk)
60+
61+ ((nominal / avgE ) .* avgEk) .* scalek
62+ end
63+
64+
65+ """
66+ _pixvalWindow
67+
68+ A utility weighting function for HDR smoothness regularization.
69+
70+ See: solveDetectorResponse for more details.
71+ """
72+ function _pixvalWindow (
73+ z:: Int ;
74+ zmin:: Int = 0 ,
75+ zmax:: Int = 2 ^ 8 - 1 ,
76+ )
77+ if z <= (zmin + zmax)/ 2
78+ return z - zmin
79+ else
80+ return zmax - z
81+ end
82+ end
83+
84+
85+ """
86+ solveDetectorResponse
87+
88+ Solve for imaging system response function, g[z] which is the log exposure
89+ corresponding to pixel value z:
90+ g := ln(f^-1) : pixelval --> log irradiance + log exposure.
91+
92+ This procedure uses various pixel values Z[j][i] from a set of physical scene
93+ objects i, as repeatedly observed with different accumulated energies j.
94+ logΔTs[j] is a vector of the log delta t / shutter speeds / ambient radiance.
95+ Note, log exposure[i] is the log film irradiance at pixel location i
96+
97+ The core assumption is that the different pixel values from
98+ different images can be attributed to the same underlying physical
99+ scene object with known received energy (exposure time). This allows for an
100+ optimization to be performed to solve for the imaging system
101+ response function gcurve. Various physical scene locations i intend to
102+ make the whole pixel value range of gcurve to be sufficiently observed.
103+
104+ Unscaled radiance is reconstructed using (gcurve[zij] - logΔts[j]) or better,
105+ - see `recoverRadianceWeighted`, `_pixvalWindow`
106+
107+ How different pixels of the same physical scene object are observed is
108+ less important, barring image warping/perspective limitations. This procedure
109+ is not capable of solving both object albedo and environmental changes simultaneously.
110+
111+ Returns a tuple (gcurve, logExposure)
112+
113+ Parameters:
114+ - λ::Real is the regularization constant rewarding camera response curve smoothness.
115+ - window: z --> weight, relative importance of different mid-range pixel values for the optimization.
116+ - see _pixvalWindow(z) for a common weighting function that de-emphasizes near saturation pixel values.
117+
118+
119+ Usage example:
120+ ```julia
121+ # stack images of the same scene with different exposures
122+ imgs = [img1, img2, img3]
123+ # extract average exposure as proxy for shutter speeds
124+ est_lΔT = log.(CameraModels.estimateShutterTimes_RGB(imgs))
125+ # Z[j][i] is the pixel values of pixel location number i in image j
126+ _getpixel(mg, ij) = mg[ij...]
127+ _collectpixs(mg) = _getpixel.(Ref(mg), uv)
128+ _Z = _collectpixs.(imgs)
129+
130+ # Convert FixedPoint N0f8 values to UInt8 (todo all channels)
131+ Zr = (s->((p->reinterpret(UInt8,p.r)).(s))).(_Z)
132+ gr, lEr = CameraModels.solveDetectorResponse(Zr, est_lΔT; λ)
133+ # ...
134+ normalized_image = CameraModels.recoverRadianceWeighted_RGB(imgs, est_lΔT, [gr, gg, gb])
135+ ```
136+
137+ See also: `randomPixels`, `normalizeRadianceMap`, `recoverRadianceWeighted_RGB`, `_pixvalWindow`
138+ """
139+ function solveDetectorResponse (
140+ Z:: AbstractVector ,
141+ logΔTs:: AbstractVector ;
142+ n:: Int = 2 ^ 8 , # common for 3*8 bit RGB;
143+ window:: AbstractVector{<:Real} = _pixvalWindow .(0 : (n- 1 )),
144+ λ:: Real = 200.0 ,
145+ diff_kernel:: AbstractVector{<:Real} = SA[1 , - 2 , 1 ],
146+ )
147+ # System size
148+ nimgs = length (Z)
149+ nlocs = length (Z[1 ])
150+
151+ # Prepare a linear system.
152+ # The total number of equations (neqs), which includes:
153+ # - nlocs*nimgs equations from pixel value observations
154+ # - 1 equation from fixing gcurve middle value to 0 for scale
155+ # - n equations from smoothness regularization
156+ # A is size (neqs, Zmax-Zmin+pixel_locations)
157+ # A = [A1 A2; A3 0], where
158+ # A{1,2} is has rows that each contain two column entries such that, gcurve[zij] = logExposure_ij + logΔt_j
159+ # - a1: the weight at column z (observed pixel value) and next available column
160+ # - a2: the negative weight at column n+i (corresponding to list_i of pixel locations)
161+ # A2 is also off-diagonal
162+ # A3 is the weighted regularization as off diagnonal entries for smoothness of gcurve, and
163+ # b has length (neqs) = [weighted logΔTs; zeros(n)]
164+ # x = [gcurve; = A * b
165+ # logExposure]
166+ #
167+ neqs = nlocs* nimgs + 1 + n
168+ A = zeros (neqs,n+ nlocs)
169+ b = zeros (neqs)
170+
171+ # Fill in the equations from pixel value observations in A and b,
172+ # Note the number of equations is nlocs*nimgs, where each pixel value observation contributes one equation.
173+ k = 1 ;
174+ for i in 1 : nlocs, j in 1 : nimgs
175+ # A1 columns correspond to gcurve possible pixel z values,
176+ # A2 columns correspond to logExposure of observed pixels from list of locations i
177+ pixz1 = Z[j][i] + 1
178+ a_12 = view (A, k, SA[pixz1, n+ i])
179+ a_12 .= window[pixz1] .* SA[1 , - 1 ]
180+ b[k] = window[pixz1] * logΔTs[j]
181+ k += 1
182+ end
183+
184+ # Add the equation to fix the curve by setting its middle value to 0
185+ # Assumption, middle value between (Zmax-min)/2 + 1, e.g. [0..255]
186+ A[k,Int (n// 2 + 1 )] = 1 ;
187+ k += 1
188+
189+ # Add n equations for smoothness regularization, which encourages the
190+ # response curve to be smooth by penalizing second derivatives.
191+ for i in 2 : n− 1
192+ # A3 off-diagonal entries numerically approximate the second derivative of gcurve at possible pixel values
193+ _A = view (A, k, (i- 1 ): (i+ 1 ))
194+ _A .= λ* window[i] .* diff_kernel
195+ k += 1
196+ end
197+
198+ # Solve linear system via Cholesky, QR, similar (internal Julia solver)
199+ x = A\ b
200+
201+ # Extract gcurve and log exposure values from solution vector
202+ gcurve = x[1 : n]
203+ logExposure = x[(n+ 1 ): end ]
204+ return gcurve, logExposure
205+ end
206+
207+
208+ function normalizeRadianceMap (
209+ rm,
210+ nper:: Real = 0.02 ;
211+ upper = maximum (rm),
212+ )
213+ lower = percentile (rm[:], nper)
214+ nrm = (rm .- lower) ./ (upper - lower)
215+ nrm[nrm .< 0 ] .= 0
216+ return nrm, upper
217+ end
218+
219+
220+ function recoverRadianceWeighted (
221+ imgs,
222+ logΔts,
223+ gcurve;
224+ window:: AbstractVector{<:Real} = _pixvalWindow .(0 : (2 ^ 8 - 1 )),
225+ normalize_percentile:: Real = 0.02 , # set <= 0 to disable
226+ )
227+ resim = zeros (size (imgs[1 ])... )
228+ wmap = zeros (size (imgs[1 ])... )
229+
230+ for (k,img) in enumerate (imgs)
231+ for i in axes (img,1 )
232+ for j in axes (img,2 )
233+ pixel = img[i,j]
234+ zij = Int16 (reinterpret (UInt8,pixel))+ 1
235+ wij = window[zij]
236+ resim[i,j] += wij* (gcurve[zij] - logΔts[k])
237+ wmap[i,j] += wij
238+ end
239+ end
240+ end
241+
242+ rim = resim ./ (wmap .+ 1 )
243+ if 0 < normalize_percentile
244+ return normalizeRadianceMap (rim, normalize_percentile)[1 ]
245+ else
246+ return rim
247+ end
248+ end
249+
250+
251+ function recoverRadianceWeighted_RGB (
252+ imgs,
253+ logΔts,
254+ gcurve_rgb;
255+ window:: AbstractVector{<:Real} = _pixvalWindow .(0 : (2 ^ 8 - 1 )),
256+ normalize_percentile:: Real = 0.02 , # set <= 0 to disable
257+ )
258+ imgs_r = (_img-> ((s-> s. r). (_img))). (imgs)
259+ imgs_g = (_img-> ((s-> s. g). (_img))). (imgs)
260+ imgs_b = (_img-> ((s-> s. b). (_img))). (imgs)
261+
262+ rim_r = recoverRadianceWeighted (imgs_r, logΔts, gcurve_rgb[1 ]; window, normalize_percentile) # = 0);
263+ rim_g = recoverRadianceWeighted (imgs_g, logΔts, gcurve_rgb[2 ]; window, normalize_percentile) # = 0);
264+ rim_b = recoverRadianceWeighted (imgs_b, logΔts, gcurve_rgb[3 ]; window, normalize_percentile) # = 0);
265+
266+ # upper = maximum(vcat(rim_r[:], rim_g[:], rim_b[:]))
267+ # rim_r, = normalizeRadianceMap(rim_r, 0.02; upper)
268+ # rim_g, = normalizeRadianceMap(rim_g, 0.02; upper)
269+ # rim_b, = normalizeRadianceMap(rim_b, 0.02; upper)
270+
271+ return colorview (RGB, rim_r, rim_g, rim_b)
272+ end
273+
274+ """
275+ randomPixels
276+
277+ Utility function to sample N random pixel locations from an image, which can be used for HDR curve solving.
278+
279+ See also: `solveDetectorResponse`, `recoverRadianceWeighted`
280+ """
281+ function randomPixels (
282+ img:: AbstractMatrix ,
283+ N:: Int = 1
284+ )
285+ uu = rand (1 : size (img,1 ), N)
286+ vv = rand (1 : size (img,2 ), N)
287+ zip (uu,vv) |> collect
288+ end
0 commit comments