FFmpeg coverage

Line	Branch	Exec	Source
1			/**
2			* Copyright (c) 2015-2021, Facebook, Inc.
3			* All rights reserved.
4			*
5			* This file is part of FFmpeg.
6			*
7			* FFmpeg is free software; you can redistribute it and/or
8			* modify it under the terms of the GNU Lesser General Public
9			* License as published by the Free Software Foundation; either
10			* version 2.1 of the License, or (at your option) any later version.
11			*
12			* FFmpeg is distributed in the hope that it will be useful,
13			* but WITHOUT ANY WARRANTY; without even the implied warranty of
14			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15			* Lesser General Public License for more details.
16			*
17			* You should have received a copy of the GNU Lesser General Public
18			* License along with FFmpeg; if not, write to the Free Software
19			* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20			*/
21
22			/* Computes the Structural Similarity Metric between two 360 video streams.
23			* original SSIM algorithm:
24			* Z. Wang, A. C. Bovik, H. R. Sheikh and E. P. Simoncelli,
25			* "Image quality assessment: From error visibility to structural similarity,"
26			* IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600-612, Apr. 2004.
27			*
28			* To improve speed, this implementation uses the standard approximation of
29			* overlapped 8x8 block sums, rather than the original gaussian weights.
30			*
31			* To address warping from 360 projections for videos with same
32			* projection and resolution, the 8x8 blocks sampled are weighted by
33			* their location in the image.
34			*
35			* To apply SSIM across projections and video sizes, we render the video on to
36			* a flat "tape" from which the 8x8 are selected and compared.
37			*/
38
39			/*
40			* @file
41			* Calculate the SSIM between two input 360 videos.
42			*/
43
44			#include <math.h>
45
46			#include "libavutil/avstring.h"
47			#include "libavutil/file_open.h"
48			#include "libavutil/mem.h"
49			#include "libavutil/opt.h"
50			#include "libavutil/pixdesc.h"
51
52			#include "avfilter.h"
53			#include "drawutils.h"
54			#include "filters.h"
55			#include "framesync.h"
56
57			#define RIGHT 0
58			#define LEFT 1
59			#define TOP 2
60			#define BOTTOM 3
61			#define FRONT 4
62			#define BACK 5
63
64			#define DEFAULT_HEATMAP_W 32
65			#define DEFAULT_HEATMAP_H 16
66
67			#define M_PI_F ((float)M_PI)
68			#define M_PI_2_F ((float)M_PI_2)
69			#define M_PI_4_F ((float)M_PI_4)
70			#define M_SQRT2_F ((float)M_SQRT2)
71
72			#define DEFAULT_EXPANSION_COEF 1.01f
73
74			#define BARREL_THETA_RANGE (DEFAULT_EXPANSION_COEF * 2.0f * M_PI_F)
75			#define BARREL_PHI_RANGE (DEFAULT_EXPANSION_COEF * M_PI_2_F)
76
77			// Use fixed-point with 16 bit precision for fast bilinear math
78			#define FIXED_POINT_PRECISION 16
79
80			// Use 1MB per channel for the histogram to get 5-digit precise SSIM value
81			#define SSIM360_HIST_SIZE 131072
82
83			// The last number is a marker < 0 to mark end of list
84			static const double PERCENTILE_LIST[] = {
85			1.0, 0.9, 0.8, 0.7, 0.6,
86			0.5, 0.4, 0.3, 0.2, 0.1, 0, -1
87			};
88
89			typedef enum StereoFormat {
90			STEREO_FORMAT_TB,
91			STEREO_FORMAT_LR,
92			STEREO_FORMAT_MONO,
93			STEREO_FORMAT_N
94			} StereoFormat;
95
96			typedef enum Projection {
97			PROJECTION_CUBEMAP32,
98			PROJECTION_CUBEMAP23,
99			PROJECTION_BARREL,
100			PROJECTION_BARREL_SPLIT,
101			PROJECTION_EQUIRECT,
102			PROJECTION_N
103			} Projection;
104
105			typedef struct Map2D {
106			int w, h;
107			double *value;
108			} Map2D;
109
110			typedef struct HeatmapList {
111			Map2D map;
112			struct HeatmapList *next;
113			} HeatmapList;
114
115			typedef struct SampleParams {
116			int sample_stride;
117			int planewidth;
118			int planeheight;
119			int x_image_offset;
120			int y_image_offset;
121			int x_image_range;
122			int y_image_range;
123			int projection;
124			float expand_coef;
125			} SampleParams;
126
127			typedef struct BilinearMap {
128			// Indices to the 4 samples to compute bilinear
129			int tli;
130			int tri;
131			int bli;
132			int bri;
133
134			// Fixed point factors with which the above 4 sample vector's
135			// dot product needs to be computed for the final bilinear value
136			int tlf;
137			int trf;
138			int blf;
139			int brf;
140			} BilinearMap;
141
142			typedef struct SSIM360Context {
143			const AVClass *class;
144
145			FFFrameSync fs;
146			// Stats file configuration
147			FILE *stats_file;
148			char *stats_file_str;
149
150			// Component properties
151			int nb_components;
152			double coefs[4];
153			char comps[4];
154			int max;
155
156			// Channel configuration & properties
157			int compute_chroma;
158
159			int is_rgb;
160			uint8_t rgba_map[4];
161
162			// Standard SSIM computation configuration & workspace
163			uint64_t frame_skip_ratio;
164
165			int *temp;
166			uint64_t nb_ssim_frames;
167			uint64_t nb_net_frames;
168			double ssim360[4], ssim360_total;
169			double *ssim360_hist[4];
170			double ssim360_hist_net[4];
171			double ssim360_percentile_sum[4][256];
172
173			// 360 projection configuration & workspace
174			int ref_projection;
175			int main_projection;
176			int ref_stereo_format;
177			int main_stereo_format;
178			float ref_pad;
179			float main_pad;
180			int use_tape;
181			char *heatmap_str;
182			int default_heatmap_w;
183			int default_heatmap_h;
184
185			Map2D density;
186			HeatmapList *heatmaps;
187			int ref_planewidth[4];
188			int ref_planeheight[4];
189			int main_planewidth[4];
190			int main_planeheight[4];
191			int tape_length[4];
192			BilinearMap *ref_tape_map[4][2];
193			BilinearMap *main_tape_map[4][2];
194			float angular_resolution[4][2];
195			double (*ssim360_plane)(
196			uint8_t *main, int main_stride,
197			uint8_t *ref, int ref_stride,
198			int width, int height, void *temp,
199			int max, Map2D density);
200			} SSIM360Context;
201
202			#define OFFSET(x) offsetof(SSIM360Context, x)
203			#define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
204
205			static const AVOption ssim360_options[] = {
206			{ "stats_file", "Set file where to store per-frame difference information",
207			OFFSET(stats_file_str), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, FLAGS },
208			{ "f", "Set file where to store per-frame difference information",
209			OFFSET(stats_file_str), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, FLAGS },
210
211			{ "compute_chroma",
212			"Specifies if non-luma channels must be computed",
213			OFFSET(compute_chroma), AV_OPT_TYPE_INT, {.i64 = 1},
214			0, 1, .flags = FLAGS },
215
216			{ "frame_skip_ratio",
217			"Specifies the number of frames to be skipped from evaluation, for every evaluated frame",
218			OFFSET(frame_skip_ratio), AV_OPT_TYPE_INT, {.i64 = 0},
219			0, 1000000, .flags = FLAGS },
220
221			{ "ref_projection", "projection of the reference video",
222			OFFSET(ref_projection), AV_OPT_TYPE_INT, {.i64 = PROJECTION_EQUIRECT},
223			0, PROJECTION_N - 1, .flags = FLAGS, .unit = "projection" },
224
225			{ "e", "equirectangular", 0, AV_OPT_TYPE_CONST, {.i64 = PROJECTION_EQUIRECT}, 0, 0, FLAGS, .unit = "projection" },
226			{ "equirect", "equirectangular", 0, AV_OPT_TYPE_CONST, {.i64 = PROJECTION_EQUIRECT}, 0, 0, FLAGS, .unit = "projection" },
227			{ "c3x2", "cubemap 3x2", 0, AV_OPT_TYPE_CONST, {.i64 = PROJECTION_CUBEMAP32}, 0, 0, FLAGS, .unit = "projection" },
228			{ "c2x3", "cubemap 2x3", 0, AV_OPT_TYPE_CONST, {.i64 = PROJECTION_CUBEMAP23}, 0, 0, FLAGS, .unit = "projection" },
229			{ "barrel", "barrel facebook's 360 format", 0, AV_OPT_TYPE_CONST, {.i64 = PROJECTION_BARREL}, 0, 0, FLAGS, .unit = "projection" },
230			{ "barrelsplit", "barrel split facebook's 360 format", 0, AV_OPT_TYPE_CONST, {.i64 = PROJECTION_BARREL_SPLIT}, 0, 0, FLAGS, .unit = "projection" },
231
232			{ "main_projection", "projection of the main video",
233			OFFSET(main_projection), AV_OPT_TYPE_INT, {.i64 = PROJECTION_N},
234			0, PROJECTION_N, .flags = FLAGS, .unit = "projection" },
235
236			{ "ref_stereo", "stereo format of the reference video",
237			OFFSET(ref_stereo_format), AV_OPT_TYPE_INT, {.i64 = STEREO_FORMAT_MONO},
238			0, STEREO_FORMAT_N - 1, .flags = FLAGS, .unit = "stereo_format" },
239
240			{ "mono", NULL, 0, AV_OPT_TYPE_CONST, {.i64 = STEREO_FORMAT_MONO }, 0, 0, FLAGS, .unit = "stereo_format" },
241			{ "tb", NULL, 0, AV_OPT_TYPE_CONST, {.i64 = STEREO_FORMAT_TB }, 0, 0, FLAGS, .unit = "stereo_format" },
242			{ "lr", NULL, 0, AV_OPT_TYPE_CONST, {.i64 = STEREO_FORMAT_LR }, 0, 0, FLAGS, .unit = "stereo_format" },
243
244			{ "main_stereo", "stereo format of main video",
245			OFFSET(main_stereo_format), AV_OPT_TYPE_INT, {.i64 = STEREO_FORMAT_N},
246			0, STEREO_FORMAT_N, .flags = FLAGS, .unit = "stereo_format" },
247
248			{ "ref_pad",
249			"Expansion (padding) coefficient for each cube face of the reference video",
250			OFFSET(ref_pad), AV_OPT_TYPE_FLOAT, {.dbl = .0f}, 0, 10, .flags = FLAGS },
251
252			{ "main_pad",
253			"Expansion (padding) coefficient for each cube face of the main video",
254			OFFSET(main_pad), AV_OPT_TYPE_FLOAT, {.dbl = .0f}, 0, 10, .flags = FLAGS },
255
256			{ "use_tape",
257			"Specifies if the tape based SSIM 360 algorithm must be used independent of the input video types",
258			OFFSET(use_tape), AV_OPT_TYPE_INT, {.i64 = 0},
259			0, 1, .flags = FLAGS },
260
261			{ "heatmap_str",
262			"Heatmap data for view-based evaluation. For heatmap file format, please refer to EntSphericalVideoHeatmapData.",
263			OFFSET(heatmap_str), AV_OPT_TYPE_STRING, {.str = NULL}, 0, 0, .flags = FLAGS },
264
265			{ "default_heatmap_width",
266			"Default heatmap dimension. Will be used when dimension is not specified in heatmap data.",
267			OFFSET(default_heatmap_w), AV_OPT_TYPE_INT, {.i64 = 32}, 1, 4096, .flags = FLAGS },
268
269			{ "default_heatmap_height",
270			"Default heatmap dimension. Will be used when dimension is not specified in heatmap data.",
271			OFFSET(default_heatmap_h), AV_OPT_TYPE_INT, {.i64 = 16}, 1, 4096, .flags = FLAGS },
272
273			{ NULL }
274			};
275
276		✗	FRAMESYNC_DEFINE_CLASS(ssim360, SSIM360Context, fs);
277
278		✗	static void set_meta(AVDictionary **metadata, const char *key, char comp, float d)
279			{
280			char value[128];
281		✗	snprintf(value, sizeof(value), "%0.2f", d);
282		✗	if (comp) {
283			char key2[128];
284		✗	snprintf(key2, sizeof(key2), "%s%c", key, comp);
285		✗	av_dict_set(metadata, key2, value, 0);
286			} else {
287		✗	av_dict_set(metadata, key, value, 0);
288			}
289		✗	}
290
291		✗	static void map_uninit(Map2D *map)
292			{
293		✗	av_freep(&map->value);
294		✗	}
295
296		✗	static int map_init(Map2D *map, int w, int h)
297			{
298		✗	map->value = av_calloc(h * w, sizeof(*map->value));
299		✗	if (!map->value)
300		✗	return AVERROR(ENOMEM);
301
302		✗	map->h = h;
303		✗	map->w = w;
304
305		✗	return 0;
306			}
307
308		✗	static void map_list_free(HeatmapList **pl)
309			{
310		✗	HeatmapList *l = *pl;
311
312		✗	while (l) {
313		✗	HeatmapList *next = l->next;
314		✗	map_uninit(&l->map);
315		✗	av_freep(&l);
316		✗	l = next;
317			}
318
319		✗	*pl = NULL;
320		✗	}
321
322		✗	static int map_alloc(HeatmapList **pl, int w, int h)
323			{
324			HeatmapList *l;
325			int ret;
326
327		✗	l = av_mallocz(sizeof(*l));
328		✗	if (!l)
329		✗	return AVERROR(ENOMEM);
330
331		✗	ret = map_init(&l->map, w, h);
332		✗	if (ret < 0) {
333		✗	av_freep(&l);
334		✗	return ret;
335			}
336
337		✗	*pl = l;
338		✗	return 0;
339			}
340
341			static void
342		✗	ssim360_4x4xn_16bit(const uint8_t *main8, ptrdiff_t main_stride,
343			const uint8_t *ref8, ptrdiff_t ref_stride,
344			int64_t (*sums)[4], int width)
345			{
346		✗	const uint16_t *main16 = (const uint16_t *)main8;
347		✗	const uint16_t *ref16 = (const uint16_t *)ref8;
348
349		✗	main_stride >>= 1;
350		✗	ref_stride >>= 1;
351
352		✗	for (int z = 0; z < width; z++) {
353		✗	uint64_t s1 = 0, s2 = 0, ss = 0, s12 = 0;
354
355		✗	for (int y = 0; y < 4; y++) {
356		✗	for (int x = 0; x < 4; x++) {
357		✗	unsigned a = main16[x + y * main_stride];
358		✗	unsigned b = ref16[x + y * ref_stride];
359
360		✗	s1 += a;
361		✗	s2 += b;
362		✗	ss += a*a;
363		✗	ss += b*b;
364		✗	s12 += a*b;
365			}
366			}
367
368		✗	sums[z][0] = s1;
369		✗	sums[z][1] = s2;
370		✗	sums[z][2] = ss;
371		✗	sums[z][3] = s12;
372		✗	main16 += 4;
373		✗	ref16 += 4;
374			}
375		✗	}
376
377			static void
378		✗	ssim360_4x4xn_8bit(const uint8_t *main, ptrdiff_t main_stride,
379			const uint8_t *ref, ptrdiff_t ref_stride,
380			int (*sums)[4], int width)
381			{
382		✗	for (int z = 0; z < width; z++) {
383		✗	uint32_t s1 = 0, s2 = 0, ss = 0, s12 = 0;
384
385		✗	for (int y = 0; y < 4; y++) {
386		✗	for (int x = 0; x < 4; x++) {
387		✗	int a = main[x + y * main_stride];
388		✗	int b = ref[x + y * ref_stride];
389
390		✗	s1 += a;
391		✗	s2 += b;
392		✗	ss += a*a;
393		✗	ss += b*b;
394		✗	s12 += a*b;
395			}
396			}
397
398		✗	sums[z][0] = s1;
399		✗	sums[z][1] = s2;
400		✗	sums[z][2] = ss;
401		✗	sums[z][3] = s12;
402		✗	main += 4;
403		✗	ref += 4;
404			}
405		✗	}
406
407		✗	static float ssim360_end1x(int64_t s1, int64_t s2, int64_t ss, int64_t s12, int max)
408			{
409		✗	int64_t ssim_c1 = (int64_t)(.01 * .01 * max * max * 64 + .5);
410		✗	int64_t ssim_c2 = (int64_t)(.03 * .03 * max * max * 64 * 63 + .5);
411
412		✗	int64_t fs1 = s1;
413		✗	int64_t fs2 = s2;
414		✗	int64_t fss = ss;
415		✗	int64_t fs12 = s12;
416		✗	int64_t vars = fss * 64 - fs1 * fs1 - fs2 * fs2;
417		✗	int64_t covar = fs12 * 64 - fs1 * fs2;
418
419		✗	return (float)(2 * fs1 * fs2 + ssim_c1) * (float)(2 * covar + ssim_c2)
420		✗	/ ((float)(fs1 * fs1 + fs2 * fs2 + ssim_c1) * (float)(vars + ssim_c2));
421			}
422
423		✗	static float ssim360_end1(int s1, int s2, int ss, int s12)
424			{
425			static const int ssim_c1 = (int)(.01*.01*255*255*64 + .5);
426			static const int ssim_c2 = (int)(.03*.03*255*255*64*63 + .5);
427
428		✗	int fs1 = s1;
429		✗	int fs2 = s2;
430		✗	int fss = ss;
431		✗	int fs12 = s12;
432		✗	int vars = fss * 64 - fs1 * fs1 - fs2 * fs2;
433		✗	int covar = fs12 * 64 - fs1 * fs2;
434
435		✗	return (float)(2 * fs1 * fs2 + ssim_c1) * (float)(2 * covar + ssim_c2)
436		✗	/ ((float)(fs1 * fs1 + fs2 * fs2 + ssim_c1) * (float)(vars + ssim_c2));
437			}
438
439			static double
440		✗	ssim360_endn_16bit(const int64_t (*sum0)[4], const int64_t (*sum1)[4],
441			int width, int max,
442			double *density_map, int map_width, double *total_weight)
443			{
444		✗	double ssim360 = 0.0, weight;
445
446		✗	for (int i = 0; i < width; i++) {
447		✗	weight = density_map ? density_map[(int) ((0.5 + i) / width * map_width)] : 1.0;
448		✗	ssim360 += weight * ssim360_end1x(
449		✗	sum0[i][0] + sum0[i + 1][0] + sum1[i][0] + sum1[i + 1][0],
450		✗	sum0[i][1] + sum0[i + 1][1] + sum1[i][1] + sum1[i + 1][1],
451		✗	sum0[i][2] + sum0[i + 1][2] + sum1[i][2] + sum1[i + 1][2],
452		✗	sum0[i][3] + sum0[i + 1][3] + sum1[i][3] + sum1[i + 1][3],
453			max);
454		✗	*total_weight += weight;
455			}
456		✗	return ssim360;
457			}
458
459			static double
460		✗	ssim360_endn_8bit(const int (*sum0)[4], const int (*sum1)[4], int width,
461			double *density_map, int map_width, double *total_weight)
462			{
463		✗	double ssim360 = 0.0, weight;
464
465		✗	for (int i = 0; i < width; i++) {
466		✗	weight = density_map ? density_map[(int) ((0.5 + i) / width * map_width)] : 1.0;
467		✗	ssim360 += weight * ssim360_end1(
468		✗	sum0[i][0] + sum0[i + 1][0] + sum1[i][0] + sum1[i + 1][0],
469		✗	sum0[i][1] + sum0[i + 1][1] + sum1[i][1] + sum1[i + 1][1],
470		✗	sum0[i][2] + sum0[i + 1][2] + sum1[i][2] + sum1[i + 1][2],
471		✗	sum0[i][3] + sum0[i + 1][3] + sum1[i][3] + sum1[i + 1][3]);
472		✗	*total_weight += weight;
473			}
474		✗	return ssim360;
475			}
476
477			static double
478		✗	ssim360_plane_16bit(uint8_t *main, int main_stride,
479			uint8_t *ref, int ref_stride,
480			int width, int height, void *temp,
481			int max, Map2D density)
482			{
483		✗	int z = 0;
484		✗	double ssim360 = 0.0;
485		✗	int64_t (*sum0)[4] = temp;
486		✗	int64_t (*sum1)[4] = sum0 + (width >> 2) + 3;
487		✗	double total_weight = 0.0;
488
489		✗	width >>= 2;
490		✗	height >>= 2;
491
492		✗	for (int y = 1; y < height; y++) {
493		✗	for (; z <= y; z++) {
494		✗	FFSWAP(void*, sum0, sum1);
495		✗	ssim360_4x4xn_16bit(&main[4 * z * main_stride], main_stride,
496		✗	&ref[4 * z * ref_stride], ref_stride,
497			sum0, width);
498			}
499		✗	ssim360 += ssim360_endn_16bit(
500			(const int64_t (*)[4])sum0, (const int64_t (*)[4])sum1,
501			width - 1, max,
502		✗	density.value ? density.value + density.w * ((int) ((z - 1.0) / height * density.h)) : NULL,
503			density.w, &total_weight);
504			}
505
506		✗	return (double) (ssim360 / total_weight);
507			}
508
509			static double
510		✗	ssim360_plane_8bit(uint8_t *main, int main_stride,
511			uint8_t *ref, int ref_stride,
512			int width, int height, void *temp,
513			int max, Map2D density)
514			{
515		✗	int z = 0;
516		✗	double ssim360 = 0.0;
517		✗	int (*sum0)[4] = temp;
518		✗	int (*sum1)[4] = sum0 + (width >> 2) + 3;
519		✗	double total_weight = 0.0;
520
521		✗	width >>= 2;
522		✗	height >>= 2;
523
524		✗	for (int y = 1; y < height; y++) {
525		✗	for (; z <= y; z++) {
526		✗	FFSWAP(void*, sum0, sum1);
527		✗	ssim360_4x4xn_8bit(
528		✗	&main[4 * z * main_stride], main_stride,
529		✗	&ref[4 * z * ref_stride], ref_stride,
530			sum0, width);
531			}
532		✗	ssim360 += ssim360_endn_8bit(
533			(const int (*)[4])sum0, (const int (*)[4])sum1, width - 1,
534		✗	density.value ? density.value + density.w * ((int) ((z - 1.0) / height * density.h)) : NULL,
535			density.w, &total_weight);
536			}
537
538		✗	return (double) (ssim360 / total_weight);
539			}
540
541		✗	static double ssim360_db(double ssim360, double weight)
542			{
543		✗	return 10 * log10(weight / (weight - ssim360));
544			}
545
546		✗	static int get_bilinear_sample(const uint8_t *data, BilinearMap *m, int max_value)
547			{
548			static const int fixed_point_half = 1 << (FIXED_POINT_PRECISION - 1);
549			static const int inv_byte_mask = UINT_MAX << 8;
550
551			int tl, tr, bl, br, v;
552
553		✗	if (max_value & inv_byte_mask) {
554		✗	uint16_t *data16 = (uint16_t *)data;
555		✗	tl = data16[m->tli];
556		✗	tr = data16[m->tri];
557		✗	bl = data16[m->bli];
558		✗	br = data16[m->bri];
559			} else {
560		✗	tl = data[m->tli];
561		✗	tr = data[m->tri];
562		✗	bl = data[m->bli];
563		✗	br = data[m->bri];
564			}
565
566		✗	v = m->tlf * tl +
567		✗	m->trf * tr +
568		✗	m->blf * bl +
569		✗	m->brf * br;
570
571			// Round by half, and revert the fixed-point offset
572		✗	return ((v + fixed_point_half) >> FIXED_POINT_PRECISION) & max_value;
573			}
574
575			static void
576		✗	ssim360_4x4x2_tape(const uint8_t *main, BilinearMap *main_maps,
577			const uint8_t *ref, BilinearMap *ref_maps,
578			int offset_y, int max_value, int (*sums)[4])
579			{
580		✗	int offset_x = 0;
581
582			// Two blocks along the width
583		✗	for (int z = 0; z < 2; z++) {
584		✗	int s1 = 0, s2 = 0, ss = 0, s12 = 0;
585
586			// 4 pixel block from (offset_x, offset_y)
587		✗	for (int y = offset_y; y < offset_y + 4; y++) {
588		✗	int y_stride = y << 3;
589		✗	for (int x = offset_x; x < offset_x + 4; x++) {
590		✗	int map_index = x + y_stride;
591		✗	int a = get_bilinear_sample(main, main_maps + map_index, max_value);
592		✗	int b = get_bilinear_sample(ref, ref_maps + map_index, max_value);
593
594		✗	s1 += a;
595		✗	s2 += b;
596		✗	ss += a*a;
597		✗	ss += b*b;
598		✗	s12 += a*b;
599			}
600			}
601
602		✗	sums[z][0] = s1;
603		✗	sums[z][1] = s2;
604		✗	sums[z][2] = ss;
605		✗	sums[z][3] = s12;
606
607		✗	offset_x += 4;
608			}
609		✗	}
610
611		✗	static float get_radius_between_negative_and_positive_pi(float theta)
612			{
613			int floor_theta_by_2pi, floor_theta_by_pi;
614
615			// Convert theta to range [0, 2*pi]
616		✗	floor_theta_by_2pi = (int)(theta / (2.0f * M_PI_F)) - (theta < 0.0f);
617		✗	theta -= 2.0f * M_PI_F * floor_theta_by_2pi;
618
619			// Convert theta to range [-pi, pi]
620		✗	floor_theta_by_pi = theta / M_PI_F;
621		✗	theta -= 2.0f * M_PI_F * floor_theta_by_pi;
622		✗	return FFMIN(M_PI_F, FFMAX(-M_PI_F, theta));
623			}
624
625		✗	static float get_heat(HeatmapList *heatmaps, float angular_resoluation, float norm_tape_pos)
626			{
627			float pitch, yaw, norm_pitch, norm_yaw;
628			int w, h;
629
630		✗	if (!heatmaps)
631		✗	return 1.0f;
632
633		✗	pitch = asinf(norm_tape_pos*2);
634		✗	yaw = M_PI_2_F * pitch / angular_resoluation;
635		✗	yaw = get_radius_between_negative_and_positive_pi(yaw);
636
637			// normalize into [0,1]
638		✗	norm_pitch = 1.0f - (pitch / M_PI_F + 0.5f);
639		✗	norm_yaw = yaw / 2.0f / M_PI_F + 0.5f;
640
641			// get heat on map
642		✗	w = FFMIN(heatmaps->map.w - 1, FFMAX(0, heatmaps->map.w * norm_yaw));
643		✗	h = FFMIN(heatmaps->map.h - 1, FFMAX(0, heatmaps->map.h * norm_pitch));
644		✗	return heatmaps->map.value[h * heatmaps->map.w + w];
645			}
646
647			static double
648		✗	ssim360_tape(uint8_t *main, BilinearMap *main_maps,
649			uint8_t *ref, BilinearMap *ref_maps,
650			int tape_length, int max_value, void *temp,
651			double *ssim360_hist, double *ssim360_hist_net,
652			float angular_resolution, HeatmapList *heatmaps)
653			{
654		✗	int horizontal_block_count = 2;
655		✗	int vertical_block_count = tape_length >> 2;
656
657		✗	int z = 0, y;
658			// Since the tape will be very long and we need to average over all 8x8 blocks, use double
659		✗	double ssim360 = 0.0;
660		✗	double sum_weight = 0.0;
661
662		✗	int (*sum0)[4] = temp;
663		✗	int (*sum1)[4] = sum0 + horizontal_block_count + 3;
664
665		✗	for (y = 1; y < vertical_block_count; y++) {
666			int fs1, fs2, fss, fs12, hist_index;
667			float norm_tape_pos, weight;
668			double sample_ssim360;
669
670		✗	for (; z <= y; z++) {
671		✗	FFSWAP(void*, sum0, sum1);
672		✗	ssim360_4x4x2_tape(main, main_maps, ref, ref_maps, z*4, max_value, sum0);
673			}
674
675			// Given we have only one 8x8 block, following sums fit within 26 bits even for 10bit videos
676		✗	fs1 = sum0[0][0] + sum0[1][0] + sum1[0][0] + sum1[1][0];
677		✗	fs2 = sum0[0][1] + sum0[1][1] + sum1[0][1] + sum1[1][1];
678		✗	fss = sum0[0][2] + sum0[1][2] + sum1[0][2] + sum1[1][2];
679		✗	fs12 = sum0[0][3] + sum0[1][3] + sum1[0][3] + sum1[1][3];
680
681		✗	if (max_value > 255) {
682			// Since we need high precision to multiply fss / fs12 by 64, use double
683		✗	double ssim_c1_d = .01*.01*64*max_value*max_value;
684		✗	double ssim_c2_d = .03*.03*64*63*max_value*max_value;
685
686		✗	double vars = 64. * fss - 1. * fs1 * fs1 - 1. * fs2 * fs2;
687		✗	double covar = 64. * fs12 - 1.*fs1 * fs2;
688		✗	sample_ssim360 = (2. * fs1 * fs2 + ssim_c1_d) * (2. * covar + ssim_c2_d)
689		✗	/ ((1. * fs1 * fs1 + 1. * fs2 * fs2 + ssim_c1_d) * (1. * vars + ssim_c2_d));
690			} else {
691			static const int ssim_c1 = (int)(.01*.01*255*255*64 + .5);
692			static const int ssim_c2 = (int)(.03*.03*255*255*64*63 + .5);
693
694		✗	int vars = fss * 64 - fs1 * fs1 - fs2 * fs2;
695		✗	int covar = fs12 * 64 - fs1 * fs2;
696		✗	sample_ssim360 = (double)(2 * fs1 * fs2 + ssim_c1) * (double)(2 * covar + ssim_c2)
697		✗	/ ((double)(fs1 * fs1 + fs2 * fs2 + ssim_c1) * (double)(vars + ssim_c2));
698			}
699
700		✗	hist_index = (int)(sample_ssim360 * ((double)SSIM360_HIST_SIZE - .5));
701		✗	hist_index = av_clip(hist_index, 0, SSIM360_HIST_SIZE - 1);
702
703		✗	norm_tape_pos = (y - 0.5f) / (vertical_block_count - 1.0f) - 0.5f;
704			// weight from an input heatmap if available, otherwise weight = 1.0
705		✗	weight = get_heat(heatmaps, angular_resolution, norm_tape_pos);
706		✗	ssim360_hist[hist_index] += weight;
707		✗	*ssim360_hist_net += weight;
708
709		✗	ssim360 += (sample_ssim360 * weight);
710		✗	sum_weight += weight;
711			}
712
713		✗	return ssim360 / sum_weight;
714			}
715
716		✗	static void compute_bilinear_map(SampleParams *p, BilinearMap *m, float x, float y)
717			{
718		✗	float fixed_point_scale = (float)(1 << FIXED_POINT_PRECISION);
719
720			// All operations in here will fit in the 22 bit mantissa of floating point,
721			// since the fixed point precision is well under 22 bits
722		✗	float x_image = av_clipf(x * p->x_image_range, 0, p->x_image_range) + p->x_image_offset;
723		✗	float y_image = av_clipf(y * p->y_image_range, 0, p->y_image_range) + p->y_image_offset;
724
725		✗	int x_floor = x_image;
726		✗	int y_floor = y_image;
727		✗	float x_diff = x_image - x_floor;
728		✗	float y_diff = y_image - y_floor;
729
730		✗	int x_ceil = x_floor + (x_diff > 1e-6);
731		✗	int y_ceil = y_floor + (y_diff > 1e-6);
732		✗	float x_inv_diff = 1.0f - x_diff;
733		✗	float y_inv_diff = 1.0f - y_diff;
734
735			// Indices of the 4 samples from source frame
736		✗	m->tli = x_floor + y_floor * p->sample_stride;
737		✗	m->tri = x_ceil + y_floor * p->sample_stride;
738		✗	m->bli = x_floor + y_ceil * p->sample_stride;
739		✗	m->bri = x_ceil + y_ceil * p->sample_stride;
740
741			// Scale to be applied to each of the 4 samples from source frame
742		✗	m->tlf = x_inv_diff * y_inv_diff * fixed_point_scale;
743		✗	m->trf = x_diff * y_inv_diff * fixed_point_scale;
744		✗	m->blf = x_inv_diff * y_diff * fixed_point_scale;
745		✗	m->brf = x_diff * y_diff * fixed_point_scale;
746		✗	}
747
748		✗	static void get_equirect_map(float phi, float theta, float *x, float *y)
749			{
750		✗	*x = 0.5f + theta / (2.0f * M_PI_F);
751			// y increases downwards
752		✗	*y = 0.5f - phi / M_PI_F;
753		✗	}
754
755		✗	static void get_barrel_map(float phi, float theta, float *x, float *y)
756			{
757		✗	float abs_phi = FFABS(phi);
758
759		✗	if (abs_phi <= M_PI_4_F) {
760			// Equirect region
761		✗	*x = 0.8f * (0.5f + theta / BARREL_THETA_RANGE);
762			// y increases downwards
763		✗	*y = 0.5f - phi / BARREL_PHI_RANGE;
764			} else {
765			// Radial ratio on a unit circle = cot(abs_phi) / (expansion_cefficient).
766			// Using cos(abs_phi)/sin(abs_phi) explicitly to avoid division by zero
767		✗	float radial_ratio = cosf(abs_phi) / (sinf(abs_phi) * DEFAULT_EXPANSION_COEF);
768		✗	float circle_x = radial_ratio * sinf(theta);
769		✗	float circle_y = radial_ratio * cosf(theta);
770		✗	float offset_y = 0.25f;
771		✗	if (phi < 0) {
772			// Bottom circle: theta increases clockwise, and front is upward
773		✗	circle_y *= -1.0f;
774		✗	offset_y += 0.5f;
775			}
776
777		✗	*x = 0.8f + 0.1f * (1.0f + circle_x);
778		✗	*y = offset_y + 0.25f * circle_y;
779			}
780		✗	}
781
782		✗	static void get_barrel_split_map(float phi, float theta, float expand_coef, float *x, float *y)
783			{
784		✗	float abs_phi = FFABS(phi);
785
786			// Front Face [-PI/2, PI/2] -> [0,1].
787			// Back Face [PI/2, PI] and [-PI, -PI/2] -> [1, 2]
788		✗	float radian_pi_theta = theta / M_PI_F + 0.5f;
789			int vFace;
790
791		✗	if (radian_pi_theta < 0.0f)
792		✗	radian_pi_theta += 2.0f;
793
794			// Front face at top (= 0), back face at bottom (= 1).
795		✗	vFace = radian_pi_theta >= 1.0f;
796
797		✗	if (abs_phi <= M_PI_4_F) {
798			// Equirect region
799		✗	*x = 2.0f / 3.0f * (0.5f + (radian_pi_theta - vFace - 0.5f) / expand_coef);
800			// y increases downwards
801		✗	*y = 0.25f + 0.5f * vFace - phi / (M_PI_F * expand_coef);
802			} else {
803			// Radial ratio on a unit circle = cot(abs_phi) / (expansion_cefficient).
804			// Using cos(abs_phi)/sin(abs_phi) explicitly to avoid division by zero
805		✗	float radial_ratio = cosf(abs_phi) / (sinf(abs_phi) * expand_coef);
806		✗	float circle_x = radial_ratio * sinf(theta);
807		✗	float circle_y = radial_ratio * cosf(theta);
808		✗	float offset_y = 0.25f;
809
810		✗	if (vFace == 1) {
811			// Back Face: Flip
812		✗	circle_x *= -1.0f;
813		✗	circle_y = (circle_y >= 0.0f) ? (1 - circle_y) : (-1 - circle_y);
814		✗	offset_y += 0.5f;
815
816			// Bottom circle: theta increases clockwise
817		✗	if (phi < 0)
818		✗	circle_y *= -1.0f;
819			} else {
820			// Front Face
821			// Bottom circle: theta increases clockwise
822		✗	if (phi < 0)
823		✗	circle_y *= -1.0f;
824			}
825
826		✗	*x = 2.0f / 3.0f + 0.5f / 3.0f * (1.0f + circle_x);
827		✗	*y = offset_y + 0.25f * circle_y / expand_coef; // y direction of expand_coeff (margin)
828			}
829		✗	}
830
831			// Returns cube face, and provided face_x & face_y will range from [0, 1]
832		✗	static int get_cubemap_face_map(float axis_vec_x, float axis_vec_y, float axis_vec_z, float *face_x, float *face_y)
833			{
834			// To check if phi, theta hits the top / bottom faces, we check the hit point of
835			// the axis vector on planes y = 1 and y = -1, and see if x & z are within [-1, 1]
836
837			// 0.577 < 1 / sqrt(3), which is less than the smallest sin(phi) falling on top/bottom faces
838			// This angle check will save computation from unnecessarily checking the top/bottom faces
839		✗	if (FFABS(axis_vec_y) > 0.577f) {
840		✗	float x_hit = axis_vec_x / FFABS(axis_vec_y);
841		✗	float z_hit = axis_vec_z / axis_vec_y;
842
843		✗	if (FFABS(x_hit) <= 1.f && FFABS(z_hit) <= 1.f) {
844		✗	*face_x = x_hit;
845			// y increases downwards
846		✗	*face_y = z_hit;
847		✗	return axis_vec_y > 0 ? TOP : BOTTOM;
848			}
849			}
850
851			// Check for left / right faces
852		✗	if (FFABS(axis_vec_x) > 0.577f) {
853		✗	float z_hit = -axis_vec_z / axis_vec_x;
854		✗	float y_hit = axis_vec_y / FFABS(axis_vec_x);
855
856		✗	if (FFABS(z_hit) <= 1.f && FFABS(y_hit) <= 1.f) {
857		✗	*face_x = z_hit;
858			// y increases downwards
859		✗	*face_y = -y_hit;
860		✗	return axis_vec_x > 0 ? RIGHT : LEFT;
861			}
862			}
863
864			// Front / back faces
865		✗	*face_x = axis_vec_x / axis_vec_z;
866			// y increases downwards
867		✗	*face_y = -axis_vec_y / FFABS(axis_vec_z);
868
869		✗	return axis_vec_z > 0 ? FRONT : BACK;
870			}
871
872		✗	static void get_cubemap32_map(float phi, float theta, float *x, float *y)
873			{
874			// face_projection_map maps each cube face to an index representing the face on the projection
875			// The indices 0->5 for cubemap 32 goes as:
876			// [0, 1, 2] as row 1, left to right
877			// [3, 4, 5] as row 2, left to right
878			static const int face_projection_map[] = {
879			[RIGHT] = 0, [LEFT] = 1, [TOP] = 2,
880			[BOTTOM] = 3, [FRONT] = 4, [BACK] = 5,
881			};
882
883		✗	float axis_vec_x = cosf(phi) * sinf(theta);
884		✗	float axis_vec_y = sinf(phi);
885		✗	float axis_vec_z = cosf(phi) * cosf(theta);
886		✗	float face_x = 0, face_y = 0;
887		✗	int face_index = get_cubemap_face_map(axis_vec_x, axis_vec_y, axis_vec_z, &face_x, &face_y);
888
889		✗	float x_offset = 1.f / 3.f * (face_projection_map[face_index] % 3);
890		✗	float y_offset = .5f * (face_projection_map[face_index] / 3);
891
892		✗	*x = x_offset + (face_x / DEFAULT_EXPANSION_COEF + 1.f) / 6.f;
893		✗	*y = y_offset + (face_y / DEFAULT_EXPANSION_COEF + 1.f) / 4.f;
894		✗	}
895
896		✗	static void get_rotated_cubemap_map(float phi, float theta, float expand_coef, float *x, float *y)
897			{
898			// face_projection_map maps each cube face to an index representing the face on the projection
899			// The indices 0->5 for rotated cubemap goes as:
900			// [0, 1] as row 1, left to right
901			// [2, 3] as row 2, left to right
902			// [4, 5] as row 3, left to right
903			static const int face_projection_map[] = {
904			[LEFT] = 0, [TOP] = 1,
905			[FRONT] = 2, [BACK] = 3,
906			[RIGHT] = 4, [BOTTOM] = 5,
907			};
908
909			float axis_yaw_vec_x, axis_yaw_vec_y, axis_yaw_vec_z;
910			float axis_pitch_vec_z, axis_pitch_vec_y;
911			float x_offset, y_offset;
912		✗	float face_x = 0, face_y = 0;
913			int face_index;
914
915			// Unrotate the cube and fix the face map:
916			// First undo the 45 degree yaw
917		✗	theta += M_PI_4_F;
918
919			// Now we are looking at the middle of an edge. So convert to axis vector & undo the pitch
920		✗	axis_yaw_vec_x = cosf(phi) * sinf(theta);
921		✗	axis_yaw_vec_y = sinf(phi);
922		✗	axis_yaw_vec_z = cosf(phi) * cosf(theta);
923
924			// The pitch axis is along +x, and has value of -45 degree. So, only y and z components change
925		✗	axis_pitch_vec_z = (axis_yaw_vec_z - axis_yaw_vec_y) / M_SQRT2_F;
926		✗	axis_pitch_vec_y = (axis_yaw_vec_y + axis_yaw_vec_z) / M_SQRT2_F;
927
928		✗	face_index = get_cubemap_face_map(axis_yaw_vec_x, axis_pitch_vec_y, axis_pitch_vec_z, &face_x, &face_y);
929
930			// Correct for the orientation of the axes on the faces
931		✗	if (face_index == LEFT || face_index == FRONT || face_index == RIGHT) {
932			// x increases downwards & y increases towards left
933		✗	float upright_y = face_y;
934		✗	face_y = face_x;
935		✗	face_x = -upright_y;
936		✗	} else if (face_index == TOP || face_index == BOTTOM) {
937			// turn the face upside-down for top and bottom
938		✗	face_x *= -1.f;
939		✗	face_y *= -1.f;
940			}
941
942		✗	x_offset = .5f * (face_projection_map[face_index] & 1);
943		✗	y_offset = 1.f / 3.f * (face_projection_map[face_index] >> 1);
944
945		✗	*x = x_offset + (face_x / expand_coef + 1.f) / 4.f;
946		✗	*y = y_offset + (face_y / expand_coef + 1.f) / 6.f;
947		✗	}
948
949		✗	static void get_projected_map(float phi, float theta, SampleParams *p, BilinearMap *m)
950			{
951		✗	float x = 0, y = 0;
952		✗	switch(p->projection) {
953			// TODO: Calculate for CDS
954		✗	case PROJECTION_CUBEMAP23:
955		✗	get_rotated_cubemap_map(phi, theta, p->expand_coef, &x, &y);
956		✗	break;
957		✗	case PROJECTION_CUBEMAP32:
958		✗	get_cubemap32_map(phi, theta, &x, &y);
959		✗	break;
960		✗	case PROJECTION_BARREL:
961		✗	get_barrel_map(phi, theta, &x, &y);
962		✗	break;
963		✗	case PROJECTION_BARREL_SPLIT:
964		✗	get_barrel_split_map(phi, theta, p->expand_coef, &x, &y);
965		✗	break;
966			// Assume PROJECTION_EQUIRECT as the default
967		✗	case PROJECTION_EQUIRECT:
968			default:
969		✗	get_equirect_map(phi, theta, &x, &y);
970		✗	break;
971			}
972		✗	compute_bilinear_map(p, m, x, y);
973		✗	}
974
975		✗	static int tape_supports_projection(int projection)
976			{
977		✗	switch(projection) {
978		✗	case PROJECTION_CUBEMAP23:
979			case PROJECTION_CUBEMAP32:
980			case PROJECTION_BARREL:
981			case PROJECTION_BARREL_SPLIT:
982			case PROJECTION_EQUIRECT:
983		✗	return 1;
984		✗	default:
985		✗	return 0;
986			}
987			}
988
989		✗	static float get_tape_angular_resolution(int projection, float expand_coef, int image_width, int image_height)
990			{
991			// NOTE: The angular resolution of a projected sphere is defined as
992			// the maximum possible horizontal angle of a pixel on the equator.
993			// We apply an intentional bias to the horizon as opposed to the meridian,
994			// since the view direction of most content is rarely closer to the poles
995
996		✗	switch(projection) {
997			// TODO: Calculate for CDS
998		✗	case PROJECTION_CUBEMAP23:
999			// Approximating atanf(pixel_width / (half_edge_width * sqrt2)) = pixel_width / (half_face_width * sqrt2)
1000		✗	return expand_coef / (M_SQRT2_F * image_width / 4.f);
1001		✗	case PROJECTION_CUBEMAP32:
1002			// Approximating atanf(pixel_width / half_face_width) = pixel_width / half_face_width
1003		✗	return DEFAULT_EXPANSION_COEF / (image_width / 6.f);
1004		✗	case PROJECTION_BARREL:
1005		✗	return FFMAX(BARREL_THETA_RANGE / (0.8f * image_width), BARREL_PHI_RANGE / image_height);
1006		✗	case PROJECTION_BARREL_SPLIT:
1007		✗	return FFMAX((expand_coef * M_PI_F) / (2.0f / 3.0f * image_width),
1008			expand_coef * M_PI_2_F / (image_height / 2.0f));
1009			// Assume PROJECTION_EQUIRECT as the default
1010		✗	case PROJECTION_EQUIRECT:
1011			default:
1012		✗	return FFMAX(2.0f * M_PI_F / image_width, M_PI_F / image_height);
1013			}
1014			}
1015
1016			static int
1017		✗	generate_eye_tape_map(SSIM360Context *s,
1018			int plane, int eye,
1019			SampleParams *ref_sample_params,
1020			SampleParams *main_sample_params)
1021			{
1022		✗	int ref_image_width = ref_sample_params->x_image_range + 1;
1023		✗	int ref_image_height = ref_sample_params->y_image_range + 1;
1024
1025			float angular_resolution =
1026		✗	get_tape_angular_resolution(s->ref_projection, 1.f + s->ref_pad,
1027			ref_image_width, ref_image_height);
1028
1029		✗	float conversion_factor = M_PI_2_F / (angular_resolution * angular_resolution);
1030		✗	float start_phi = -M_PI_2_F + 4.0f * angular_resolution;
1031		✗	float start_x = conversion_factor * sinf(start_phi);
1032		✗	float end_phi = M_PI_2_F - 3.0f * angular_resolution;
1033		✗	float end_x = conversion_factor * sinf(end_phi);
1034		✗	float x_range = end_x - start_x;
1035
1036			// Ensure tape length is a multiple of 4, for full SSIM block coverage
1037		✗	float tape_length_f = ROUNDED_DIV(x_range, 4);
1038			int tape_length;
1039
1040		✗	if (!(tape_length_f > 0.f) || tape_length_f > INT_MAX / 4.0f)
1041		✗	return AVERROR(EINVAL);
1042
1043		✗	tape_length = s->tape_length[plane] = (int)tape_length_f << 2;
1044
1045		✗	s->ref_tape_map[plane][eye] = av_malloc_array(tape_length, 8 * sizeof(BilinearMap));
1046		✗	s->main_tape_map[plane][eye] = av_malloc_array(tape_length, 8 * sizeof(BilinearMap));
1047		✗	if (!s->ref_tape_map[plane][eye] || !s->main_tape_map[plane][eye])
1048		✗	return AVERROR(ENOMEM);
1049
1050		✗	s->angular_resolution[plane][eye] = angular_resolution;
1051
1052			// For easy memory access, we navigate the tape lengthwise on y
1053		✗	for (int y_index = 0; y_index < tape_length; y_index ++) {
1054		✗	int y_stride = y_index << 3;
1055
1056		✗	float x = start_x + x_range * (y_index / (tape_length - 1.0f));
1057			// phi will be in range [-pi/2, pi/2]
1058		✗	float mid_phi = asinf(x / conversion_factor);
1059
1060		✗	float theta = mid_phi * M_PI_2_F / angular_resolution;
1061		✗	theta = get_radius_between_negative_and_positive_pi(theta);
1062
1063		✗	for (int x_index = 0; x_index < 8; x_index ++) {
1064		✗	float phi = mid_phi + angular_resolution * (3.0f - x_index);
1065		✗	int tape_index = y_stride + x_index;
1066		✗	get_projected_map(phi, theta, ref_sample_params, &s->ref_tape_map [plane][eye][tape_index]);
1067		✗	get_projected_map(phi, theta, main_sample_params, &s->main_tape_map[plane][eye][tape_index]);
1068			}
1069			}
1070
1071		✗	return 0;
1072			}
1073
1074		✗	static int generate_tape_maps(SSIM360Context *s, AVFrame *main, const AVFrame *ref)
1075			{
1076			// A tape is a long segment with 8 pixels thickness, with the angular center at the middle (below 4th pixel).
1077			// When it takes a full loop around a sphere, it will overlap the starting point at half the width from above.
1078		✗	int ref_stereo_format = s->ref_stereo_format;
1079		✗	int main_stereo_format = s->main_stereo_format;
1080		✗	int are_both_stereo = (main_stereo_format != STEREO_FORMAT_MONO) && (ref_stereo_format != STEREO_FORMAT_MONO);
1081		✗	int min_eye_count = 1 + are_both_stereo;
1082		✗	int bytes_per_sample = 1 + (s->max > 255);
1083			int ret;
1084
1085		✗	for (int i = 0; i < s->nb_components; i ++) {
1086		✗	int ref_width = s->ref_planewidth[i];
1087		✗	int ref_height = s->ref_planeheight[i];
1088		✗	int main_width = s->main_planewidth[i];
1089		✗	int main_height = s->main_planeheight[i];
1090
1091		✗	int is_ref_LR = (ref_stereo_format == STEREO_FORMAT_LR);
1092		✗	int is_ref_TB = (ref_stereo_format == STEREO_FORMAT_TB);
1093		✗	int is_main_LR = (main_stereo_format == STEREO_FORMAT_LR);
1094		✗	int is_main_TB = (main_stereo_format == STEREO_FORMAT_TB);
1095
1096		✗	int ref_image_width = is_ref_LR ? ref_width >> 1 : ref_width;
1097		✗	int ref_image_height = is_ref_TB ? ref_height >> 1 : ref_height;
1098		✗	int main_image_width = is_main_LR ? main_width >> 1 : main_width;
1099		✗	int main_image_height = is_main_TB ? main_height >> 1 : main_height;
1100
1101		✗	for (int eye = 0; eye < min_eye_count; eye ++) {
1102		✗	SampleParams ref_sample_params = {
1103		✗	.sample_stride = ref->linesize[i] / bytes_per_sample,
1104			.planewidth = ref_width,
1105			.planeheight = ref_height,
1106		✗	.x_image_range = ref_image_width - 1,
1107		✗	.y_image_range = ref_image_height - 1,
1108		✗	.x_image_offset = is_ref_LR * eye * ref_image_width,
1109		✗	.y_image_offset = is_ref_TB * eye * ref_image_height,
1110		✗	.projection = s->ref_projection,
1111		✗	.expand_coef = 1.f + s->ref_pad,
1112			};
1113
1114		✗	SampleParams main_sample_params = {
1115		✗	.sample_stride = main->linesize[i] / bytes_per_sample,
1116			.planewidth = main_width,
1117			.planeheight = main_height,
1118		✗	.x_image_range = main_image_width - 1,
1119		✗	.y_image_range = main_image_height - 1,
1120		✗	.x_image_offset = is_main_LR * eye * main_image_width,
1121		✗	.y_image_offset = is_main_TB * eye * main_image_height,
1122		✗	.projection = s->main_projection,
1123		✗	.expand_coef = 1.f + s->main_pad,
1124			};
1125
1126		✗	ret = generate_eye_tape_map(s, i, eye, &ref_sample_params, &main_sample_params);
1127		✗	if (ret < 0)
1128		✗	return ret;
1129			}
1130			}
1131
1132		✗	return 0;
1133			}
1134
1135		✗	static int do_ssim360(FFFrameSync *fs)
1136			{
1137		✗	AVFilterContext *ctx = fs->parent;
1138		✗	SSIM360Context *s = ctx->priv;
1139			AVFrame *master, *ref;
1140			AVDictionary **metadata;
1141		✗	double c[4], ssim360v = 0.0, ssim360p50 = 0.0;
1142			int ret;
1143		✗	int need_frame_skip = s->nb_net_frames % (s->frame_skip_ratio + 1);
1144		✗	HeatmapList* h_ptr = NULL;
1145
1146		✗	ret = ff_framesync_dualinput_get(fs, &master, &ref);
1147		✗	if (ret < 0)
1148		✗	return ret;
1149
1150		✗	s->nb_net_frames++;
1151
1152		✗	if (need_frame_skip)
1153		✗	return ff_filter_frame(ctx->outputs[0], master);
1154
1155		✗	metadata = &master->metadata;
1156
1157		✗	if (s->use_tape && !s->tape_length[0]) {
1158		✗	ret = generate_tape_maps(s, master, ref);
1159		✗	if (ret < 0)
1160		✗	return ret;
1161			}
1162
1163		✗	for (int i = 0; i < s->nb_components; i++) {
1164		✗	if (s->use_tape) {
1165		✗	c[i] = ssim360_tape(master->data[i], s->main_tape_map[i][0],
1166		✗	ref->data[i], s->ref_tape_map [i][0],
1167		✗	s->tape_length[i], s->max, s->temp,
1168			s->ssim360_hist[i], &s->ssim360_hist_net[i],
1169			s->angular_resolution[i][0], s->heatmaps);
1170
1171		✗	if (s->ref_tape_map[i][1]) {
1172		✗	c[i] += ssim360_tape(master->data[i], s->main_tape_map[i][1],
1173		✗	ref->data[i], s->ref_tape_map[i][1],
1174		✗	s->tape_length[i], s->max, s->temp,
1175			s->ssim360_hist[i], &s->ssim360_hist_net[i],
1176			s->angular_resolution[i][1], s->heatmaps);
1177		✗	c[i] /= 2.f;
1178			}
1179			} else {
1180		✗	c[i] = s->ssim360_plane(master->data[i], master->linesize[i],
1181		✗	ref->data[i], ref->linesize[i],
1182			s->ref_planewidth[i], s->ref_planeheight[i],
1183		✗	s->temp, s->max, s->density);
1184			}
1185
1186		✗	s->ssim360[i] += c[i];
1187		✗	ssim360v += s->coefs[i] * c[i];
1188			}
1189
1190		✗	s->nb_ssim_frames++;
1191		✗	if (s->heatmaps) {
1192		✗	map_uninit(&s->heatmaps->map);
1193		✗	h_ptr = s->heatmaps;
1194		✗	s->heatmaps = s->heatmaps->next;
1195		✗	av_freep(&h_ptr);
1196			}
1197		✗	s->ssim360_total += ssim360v;
1198
1199			// Record percentiles from histogram and attach metadata when using tape
1200		✗	if (s->use_tape) {
1201			int hist_indices[4];
1202			double hist_weight[4];
1203
1204		✗	for (int i = 0; i < s->nb_components; i++) {
1205		✗	hist_indices[i] = SSIM360_HIST_SIZE - 1;
1206		✗	hist_weight[i] = 0;
1207			}
1208
1209		✗	for (int p = 0; PERCENTILE_LIST[p] >= 0.0; p ++) {
1210		✗	for (int i = 0; i < s->nb_components; i++) {
1211			double target_weight, ssim360p;
1212
1213			// Target weight = total number of samples above the specified percentile
1214		✗	target_weight = (1. - PERCENTILE_LIST[p]) * s->ssim360_hist_net[i];
1215		✗	target_weight = FFMAX(target_weight, 1);
1216		✗	while(hist_indices[i] >= 0 && hist_weight[i] < target_weight) {
1217		✗	hist_weight[i] += s->ssim360_hist[i][hist_indices[i]];
1218		✗	hist_indices[i] --;
1219			}
1220
1221		✗	ssim360p = (double)(hist_indices[i] + 1) / (double)(SSIM360_HIST_SIZE - 1);
1222		✗	if (PERCENTILE_LIST[p] == 0.5)
1223		✗	ssim360p50 += s->coefs[i] * ssim360p;
1224		✗	s->ssim360_percentile_sum[i][p] += ssim360p;
1225			}
1226			}
1227
1228		✗	for (int i = 0; i < s->nb_components; i++) {
1229		✗	memset(s->ssim360_hist[i], 0, SSIM360_HIST_SIZE * sizeof(double));
1230		✗	s->ssim360_hist_net[i] = 0;
1231			}
1232
1233		✗	for (int i = 0; i < s->nb_components; i++) {
1234		✗	int cidx = s->is_rgb ? s->rgba_map[i] : i;
1235		✗	set_meta(metadata, "lavfi.ssim360.", s->comps[i], c[cidx]);
1236			}
1237
1238			// Use p50 as the aggregated value
1239		✗	set_meta(metadata, "lavfi.ssim360.All", 0, ssim360p50);
1240		✗	set_meta(metadata, "lavfi.ssim360.dB", 0, ssim360_db(ssim360p50, 1.0));
1241
1242		✗	if (s->stats_file) {
1243		✗	fprintf(s->stats_file, "n:%"PRId64" ", s->nb_ssim_frames);
1244
1245		✗	for (int i = 0; i < s->nb_components; i++) {
1246		✗	int cidx = s->is_rgb ? s->rgba_map[i] : i;
1247		✗	fprintf(s->stats_file, "%c:%f ", s->comps[i], c[cidx]);
1248			}
1249
1250		✗	fprintf(s->stats_file, "All:%f (%f)\n", ssim360p50, ssim360_db(ssim360p50, 1.0));
1251			}
1252			}
1253
1254		✗	return ff_filter_frame(ctx->outputs[0], master);
1255			}
1256
1257		✗	static int parse_heatmaps(void *logctx, HeatmapList **proot,
1258			const char *data, int w, int h)
1259			{
1260		✗	HeatmapList *root = NULL;
1261		✗	HeatmapList **next = &root;
1262
1263			int ret;
1264
1265			// skip video id line
1266		✗	data = strchr(data, '\n');
1267		✗	if (!data) {
1268		✗	av_log(logctx, AV_LOG_ERROR, "Invalid heatmap syntax\n");
1269		✗	return AVERROR(EINVAL);
1270			}
1271		✗	data++;
1272
1273		✗	while (*data) {
1274			HeatmapList *cur;
1275		✗	char *line = av_get_token(&data, "\n");
1276			char *saveptr, *val;
1277			int i;
1278
1279		✗	if (!line) {
1280		✗	ret = AVERROR(ENOMEM);
1281		✗	goto fail;
1282			}
1283
1284			// first value is frame id
1285		✗	av_strtok(line, ",", &saveptr);
1286
1287		✗	ret = map_alloc(next, w, h);
1288		✗	if (ret < 0)
1289		✗	goto line_fail;
1290
1291		✗	cur = *next;
1292		✗	next = &cur->next;
1293
1294		✗	i = 0;
1295		✗	while ((val = av_strtok(NULL, ",", &saveptr))) {
1296		✗	if (i >= w * h) {
1297		✗	av_log(logctx, AV_LOG_ERROR, "Too many entries in a heat map\n");
1298		✗	ret = AVERROR(EINVAL);
1299		✗	goto line_fail;
1300			}
1301
1302		✗	cur->map.value[i++] = atof(val);
1303			}
1304
1305		✗	line_fail:
1306		✗	av_freep(&line);
1307		✗	if (ret < 0)
1308		✗	goto fail;
1309			}
1310
1311		✗	*proot = root;
1312
1313		✗	return 0;
1314		✗	fail:
1315		✗	map_list_free(&root);
1316		✗	return ret;
1317			}
1318
1319		✗	static av_cold int init(AVFilterContext *ctx)
1320			{
1321		✗	SSIM360Context *s = ctx->priv;
1322			int err;
1323
1324		✗	if (s->stats_file_str) {
1325		✗	if (!strcmp(s->stats_file_str, "-")) {
1326		✗	s->stats_file = stdout;
1327			} else {
1328		✗	s->stats_file = avpriv_fopen_utf8(s->stats_file_str, "w");
1329		✗	if (!s->stats_file) {
1330		✗	err = AVERROR(errno);
1331		✗	av_log(ctx, AV_LOG_ERROR, "Could not open stats file %s: %s\n",
1332		✗	s->stats_file_str, av_err2str(err));
1333		✗	return err;
1334			}
1335			}
1336			}
1337
1338		✗	if (s->use_tape && s->heatmap_str) {
1339		✗	err = parse_heatmaps(ctx, &s->heatmaps, s->heatmap_str,
1340			s->default_heatmap_w, s->default_heatmap_h);
1341		✗	if (err < 0)
1342		✗	return err;
1343			}
1344
1345		✗	s->fs.on_event = do_ssim360;
1346		✗	return 0;
1347			}
1348
1349		✗	static int config_input_main(AVFilterLink *inlink)
1350			{
1351		✗	const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(inlink->format);
1352		✗	AVFilterContext *ctx = inlink->dst;
1353		✗	SSIM360Context *s = ctx->priv;
1354
1355		✗	s->main_planeheight[0] = inlink->h;
1356		✗	s->main_planeheight[3] = inlink->h;
1357		✗	s->main_planeheight[1] = AV_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
1358		✗	s->main_planeheight[2] = AV_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
1359
1360		✗	s->main_planewidth[0] = inlink->w;
1361		✗	s->main_planewidth[3] = inlink->w;
1362		✗	s->main_planewidth[1] = AV_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
1363		✗	s->main_planewidth[2] = AV_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
1364
1365			// If main projection is unindentified, assume it is same as reference
1366		✗	if (s->main_projection == PROJECTION_N)
1367		✗	s->main_projection = s->ref_projection;
1368
1369			// If main stereo format is unindentified, assume it is same as reference
1370		✗	if (s->main_stereo_format == STEREO_FORMAT_N)
1371		✗	s->main_stereo_format = s->ref_stereo_format;
1372
1373		✗	return 0;
1374			}
1375
1376		✗	static int generate_density_map(SSIM360Context *s, int w, int h)
1377			{
1378			double d, r_square, cos_square;
1379			int ow, oh, ret;
1380
1381		✗	ret = map_init(&s->density, w, h);
1382		✗	if (ret < 0)
1383		✗	return ret;
1384
1385		✗	switch (s->ref_stereo_format) {
1386		✗	case STEREO_FORMAT_TB:
1387		✗	h >>= 1;
1388		✗	break;
1389		✗	case STEREO_FORMAT_LR:
1390		✗	w >>= 1;
1391		✗	break;
1392			}
1393
1394		✗	switch (s->ref_projection) {
1395		✗	case PROJECTION_EQUIRECT:
1396		✗	for (int i = 0; i < h; i++) {
1397		✗	d = cos(((0.5 + i) / h - 0.5) * M_PI);
1398		✗	for (int j = 0; j < w; j++)
1399		✗	s->density.value[i * w + j] = d;
1400			}
1401		✗	break;
1402		✗	case PROJECTION_CUBEMAP32:
1403			// for one quater of a face
1404		✗	for (int i = 0; i < h / 4; i++) {
1405		✗	for (int j = 0; j < w / 6; j++) {
1406			// r = normalized distance to the face center
1407		✗	r_square =
1408		✗	(0.5 + i) / (h / 2) * (0.5 + i) / (h / 2) +
1409		✗	(0.5 + j) / (w / 3) * (0.5 + j) / (w / 3);
1410		✗	r_square /= DEFAULT_EXPANSION_COEF * DEFAULT_EXPANSION_COEF;
1411		✗	cos_square = 0.25 / (r_square + 0.25);
1412		✗	d = pow(cos_square, 1.5);
1413
1414		✗	for (int face = 0; face < 6; face++) {
1415			// center of a face
1416		✗	switch (face) {
1417		✗	case 0:
1418		✗	oh = h / 4;
1419		✗	ow = w / 6;
1420		✗	break;
1421		✗	case 1:
1422		✗	oh = h / 4;
1423		✗	ow = w / 6 + w / 3;
1424		✗	break;
1425		✗	case 2:
1426		✗	oh = h / 4;
1427		✗	ow = w / 6 + 2 * w / 3;
1428		✗	break;
1429		✗	case 3:
1430		✗	oh = h / 4 + h / 2;
1431		✗	ow = w / 6;
1432		✗	break;
1433		✗	case 4:
1434		✗	oh = h / 4 + h / 2;
1435		✗	ow = w / 6 + w / 3;
1436		✗	break;
1437		✗	case 5:
1438		✗	oh = h / 4 + h / 2;
1439		✗	ow = w / 6 + 2 * w / 3;
1440		✗	break;
1441			}
1442		✗	s->density.value[(oh - 1 - i) * w + ow - 1 - j] = d;
1443		✗	s->density.value[(oh - 1 - i) * w + ow + j] = d;
1444		✗	s->density.value[(oh + i) * w + ow - 1 - j] = d;
1445		✗	s->density.value[(oh + i) * w + ow + j] = d;
1446			}
1447			}
1448			}
1449		✗	break;
1450		✗	case PROJECTION_CUBEMAP23:
1451			// for one quater of a face
1452		✗	for (int i = 0; i < h / 6; i++) {
1453		✗	for (int j = 0; j < w / 4; j++) {
1454			// r = normalized distance to the face center
1455		✗	r_square =
1456		✗	(0.5 + i) / (h / 3) * (0.5 + i) / (h / 3) +
1457		✗	(0.5 + j) / (w / 2) * (0.5 + j) / (w / 2);
1458		✗	r_square /= (1.f + s->ref_pad) * (1.f + s->ref_pad);
1459		✗	cos_square = 0.25 / (r_square + 0.25);
1460		✗	d = pow(cos_square, 1.5);
1461
1462		✗	for (int face = 0; face < 6; face++) {
1463			// center of a face
1464		✗	switch (face) {
1465		✗	case 0:
1466		✗	ow = w / 4;
1467		✗	oh = h / 6;
1468		✗	break;
1469		✗	case 1:
1470		✗	ow = w / 4;
1471		✗	oh = h / 6 + h / 3;
1472		✗	break;
1473		✗	case 2:
1474		✗	ow = w / 4;
1475		✗	oh = h / 6 + 2 * h / 3;
1476		✗	break;
1477		✗	case 3:
1478		✗	ow = w / 4 + w / 2;
1479		✗	oh = h / 6;
1480		✗	break;
1481		✗	case 4:
1482		✗	ow = w / 4 + w / 2;
1483		✗	oh = h / 6 + h / 3;
1484		✗	break;
1485		✗	case 5:
1486		✗	ow = w / 4 + w / 2;
1487		✗	oh = h / 6 + 2 * h / 3;
1488		✗	break;
1489			}
1490		✗	s->density.value[(oh - 1 - i) * w + ow - 1 - j] = d;
1491		✗	s->density.value[(oh - 1 - i) * w + ow + j] = d;
1492		✗	s->density.value[(oh + i) * w + ow - 1 - j] = d;
1493		✗	s->density.value[(oh + i) * w + ow + j] = d;
1494			}
1495			}
1496			}
1497		✗	break;
1498		✗	case PROJECTION_BARREL:
1499			// side face
1500		✗	for (int i = 0; i < h; i++) {
1501		✗	d = cos(((0.5 + i) / h - 0.5) * DEFAULT_EXPANSION_COEF * M_PI_2);
1502		✗	for (int j = 0; j < w * 4 / 5; j++) {
1503		✗	s->density.value[i * w + j] = d * d * d;
1504			}
1505			}
1506			// top and bottom
1507		✗	for (int i = 0; i < h; i++) {
1508		✗	for (int j = w * 4 / 5; j < w; j++) {
1509		✗	double dx = DEFAULT_EXPANSION_COEF * (0.5 + j - w * 0.90) / (w * 0.10);
1510		✗	double dx_squared = dx * dx;
1511
1512		✗	double top_dy = DEFAULT_EXPANSION_COEF * (0.5 + i - h * 0.25) / (h * 0.25);
1513		✗	double top_dy_squared = top_dy * top_dy;
1514
1515		✗	double bottom_dy = DEFAULT_EXPANSION_COEF * (0.5 + i - h * 0.75) / (h * 0.25);
1516		✗	double bottom_dy_squared = bottom_dy * bottom_dy;
1517
1518			// normalized distance to the circle center
1519		✗	r_square = (i < h / 2 ? top_dy_squared : bottom_dy_squared) + dx_squared;
1520		✗	if (r_square > 1.0)
1521		✗	continue;
1522
1523		✗	cos_square = 1.0 / (r_square + 1.0);
1524		✗	d = pow(cos_square, 1.5);
1525		✗	s->density.value[i * w + j] = d;
1526			}
1527			}
1528		✗	break;
1529		✗	default:
1530			// TODO: SSIM360_v1
1531		✗	for (int i = 0; i < h; i++) {
1532		✗	for (int j = 0; j < w; j++)
1533		✗	s->density.value[i * w + j] = 0;
1534			}
1535			}
1536
1537		✗	switch (s->ref_stereo_format) {
1538		✗	case STEREO_FORMAT_TB:
1539		✗	for (int i = 0; i < h; i++) {
1540		✗	for (int j = 0; j < w; j++)
1541		✗	s->density.value[(i + h) * w + j] = s->density.value[i * w + j];
1542			}
1543		✗	break;
1544		✗	case STEREO_FORMAT_LR:
1545		✗	for (int i = 0; i < h; i++) {
1546		✗	for (int j = 0; j < w; j++)
1547		✗	s->density.value[i * w + j + w] = s->density.value[i * w + j];
1548			}
1549			}
1550
1551		✗	return 0;
1552			}
1553
1554		✗	static int config_input_ref(AVFilterLink *inlink)
1555			{
1556		✗	const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(inlink->format);
1557		✗	AVFilterContext *ctx = inlink->dst;
1558		✗	SSIM360Context *s = ctx->priv;
1559		✗	int sum = 0;
1560
1561		✗	s->nb_components = desc->nb_components;
1562
1563		✗	s->ref_planeheight[0] = inlink->h;
1564		✗	s->ref_planeheight[3] = inlink->h;
1565		✗	s->ref_planeheight[1] = AV_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
1566		✗	s->ref_planeheight[2] = AV_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
1567
1568		✗	s->ref_planewidth[0] = inlink->w;
1569		✗	s->ref_planewidth[3] = inlink->w;
1570		✗	s->ref_planewidth[1] = AV_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
1571		✗	s->ref_planewidth[2] = AV_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
1572
1573		✗	s->is_rgb = ff_fill_rgba_map(s->rgba_map, inlink->format) >= 0;
1574		✗	s->comps[0] = s->is_rgb ? 'R' : 'Y';
1575		✗	s->comps[1] = s->is_rgb ? 'G' : 'U';
1576		✗	s->comps[2] = s->is_rgb ? 'B' : 'V';
1577		✗	s->comps[3] = 'A';
1578
1579			// If chroma computation is disabled, and the format is YUV, skip U & V channels
1580		✗	if (!s->is_rgb && !s->compute_chroma)
1581		✗	s->nb_components = 1;
1582
1583		✗	s->max = (1 << desc->comp[0].depth) - 1;
1584
1585		✗	s->ssim360_plane = desc->comp[0].depth > 8 ? ssim360_plane_16bit : ssim360_plane_8bit;
1586
1587		✗	for (int i = 0; i < s->nb_components; i++)
1588		✗	sum += s->ref_planeheight[i] * s->ref_planewidth[i];
1589		✗	for (int i = 0; i < s->nb_components; i++)
1590		✗	s->coefs[i] = (double) s->ref_planeheight[i] * s->ref_planewidth[i] / sum;
1591
1592		✗	return 0;
1593			}
1594
1595		✗	static int config_output(AVFilterLink *outlink)
1596			{
1597		✗	AVFilterContext *ctx = outlink->src;
1598		✗	SSIM360Context *s = ctx->priv;
1599		✗	AVFilterLink *mainlink = ctx->inputs[0];
1600		✗	AVFilterLink *reflink = ctx->inputs[0];
1601		✗	FilterLink *il = ff_filter_link(mainlink);
1602		✗	FilterLink *ol = ff_filter_link(outlink);
1603		✗	const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(outlink->format);
1604			int ret;
1605
1606			// Use tape algorithm if any of frame sizes, projections or stereo format are not equal
1607		✗	if (ctx->inputs[0]->w != ctx->inputs[1]->w || ctx->inputs[0]->h != ctx->inputs[1]->h ||
1608		✗	s->ref_projection != s->main_projection || s->ref_stereo_format != s->main_stereo_format)
1609		✗	s->use_tape = 1;
1610
1611			// Finally, if we have decided to / forced to use tape, check if tape supports both input and output projection
1612		✗	if (s->use_tape &&
1613		✗	!(tape_supports_projection(s->main_projection) &&
1614		✗	tape_supports_projection(s->ref_projection))) {
1615		✗	av_log(ctx, AV_LOG_ERROR, "Projection is unsupported for the tape based algorithm\n");
1616		✗	return AVERROR(EINVAL);
1617			}
1618
1619		✗	if (s->use_tape) {
1620			// s->temp will be allocated for the tape width = 8. The tape is long downwards
1621		✗	s->temp = av_malloc_array((2 * 8 + 12), sizeof(*s->temp));
1622		✗	if (!s->temp)
1623		✗	return AVERROR(ENOMEM);
1624
1625		✗	memset(s->ssim360_percentile_sum, 0, sizeof(s->ssim360_percentile_sum));
1626
1627		✗	for (int i = 0; i < s->nb_components; i++) {
1628		✗	FF_ALLOCZ_TYPED_ARRAY(s->ssim360_hist[i], SSIM360_HIST_SIZE);
1629		✗	if (!s->ssim360_hist[i])
1630		✗	return AVERROR(ENOMEM);
1631			}
1632			} else {
1633		✗	s->temp = av_malloc_array((2 * reflink->w + 12), sizeof(*s->temp) * (1 + (desc->comp[0].depth > 8)));
1634		✗	if (!s->temp)
1635		✗	return AVERROR(ENOMEM);
1636
1637		✗	if (!s->density.value) {
1638		✗	ret = generate_density_map(s, reflink->w, reflink->h);
1639		✗	if (ret < 0)
1640		✗	return ret;
1641			}
1642			}
1643
1644		✗	ret = ff_framesync_init_dualinput(&s->fs, ctx);
1645		✗	if (ret < 0)
1646		✗	return ret;
1647
1648		✗	outlink->w = mainlink->w;
1649		✗	outlink->h = mainlink->h;
1650		✗	outlink->time_base = mainlink->time_base;
1651		✗	outlink->sample_aspect_ratio = mainlink->sample_aspect_ratio;
1652		✗	ol->frame_rate = il->frame_rate;
1653
1654		✗	s->fs.opt_shortest = 1;
1655		✗	s->fs.opt_repeatlast = 1;
1656
1657		✗	ret = ff_framesync_configure(&s->fs);
1658		✗	if (ret < 0)
1659		✗	return ret;
1660
1661		✗	return 0;
1662			}
1663
1664		✗	static int activate(AVFilterContext *ctx)
1665			{
1666		✗	SSIM360Context *s = ctx->priv;
1667		✗	return ff_framesync_activate(&s->fs);
1668			}
1669
1670		✗	static av_cold void uninit(AVFilterContext *ctx)
1671			{
1672		✗	SSIM360Context *s = ctx->priv;
1673
1674		✗	if (s->nb_ssim_frames > 0) {
1675			char buf[256];
1676		✗	buf[0] = 0;
1677			// Log average SSIM360 values
1678		✗	for (int i = 0; i < s->nb_components; i++) {
1679		✗	int c = s->is_rgb ? s->rgba_map[i] : i;
1680		✗	av_strlcatf(buf, sizeof(buf), " %c:%f (%f)", s->comps[i], s->ssim360[c] / s->nb_ssim_frames,
1681		✗	ssim360_db(s->ssim360[c], s->nb_ssim_frames));
1682			}
1683		✗	av_log(ctx, AV_LOG_INFO, "SSIM360%s All:%f (%f)\n", buf,
1684		✗	s->ssim360_total / s->nb_ssim_frames, ssim360_db(s->ssim360_total, s->nb_ssim_frames));
1685
1686			// Log percentiles from histogram when using tape
1687		✗	if (s->use_tape) {
1688		✗	for (int p = 0; PERCENTILE_LIST[p] >= 0.0; p++) {
1689		✗	buf[0] = 0;
1690		✗	for (int i = 0; i < s->nb_components; i++) {
1691		✗	int c = s->is_rgb ? s->rgba_map[i] : i;
1692		✗	double ssim360p = s->ssim360_percentile_sum[i][p] / (double)(s->nb_ssim_frames);
1693		✗	av_strlcatf(buf, sizeof(buf), " %c:%f (%f)", s->comps[c], ssim360p, ssim360_db(ssim360p, 1));
1694			}
1695		✗	av_log(ctx, AV_LOG_INFO, "SSIM360_p%d%s\n", (int)(PERCENTILE_LIST[p] * 100.), buf);
1696			}
1697			}
1698			}
1699
1700			// free density map
1701		✗	map_uninit(&s->density);
1702
1703		✗	map_list_free(&s->heatmaps);
1704
1705		✗	for (int i = 0; i < s->nb_components; i++) {
1706		✗	for (int eye = 0; eye < 2; eye++) {
1707		✗	av_freep(&s->ref_tape_map[i][eye]);
1708		✗	av_freep(&s->main_tape_map[i][eye]);
1709			}
1710		✗	av_freep(&s->ssim360_hist[i]);
1711			}
1712
1713		✗	ff_framesync_uninit(&s->fs);
1714
1715		✗	if (s->stats_file && s->stats_file != stdout)
1716		✗	fclose(s->stats_file);
1717
1718		✗	av_freep(&s->temp);
1719		✗	}
1720
1721			#define PF(suf) AV_PIX_FMT_YUV420##suf, AV_PIX_FMT_YUV422##suf, AV_PIX_FMT_YUV444##suf, AV_PIX_FMT_GBR##suf
1722			static const enum AVPixelFormat ssim360_pixfmts[] = {
1723			AV_PIX_FMT_GRAY8,
1724			AV_PIX_FMT_YUV420P, AV_PIX_FMT_YUV422P, AV_PIX_FMT_YUV444P,
1725			AV_PIX_FMT_YUV440P, AV_PIX_FMT_YUV411P, AV_PIX_FMT_YUV410P,
1726			AV_PIX_FMT_YUVJ411P, AV_PIX_FMT_YUVJ420P, AV_PIX_FMT_YUVJ422P,
1727			AV_PIX_FMT_YUVJ440P, AV_PIX_FMT_YUVJ444P,
1728			AV_PIX_FMT_GBRP,
1729			PF(P9), PF(P10), PF(P12), PF(P14), PF(P16),
1730			AV_PIX_FMT_NONE
1731			};
1732			#undef PF
1733
1734			static const AVFilterPad ssim360_inputs[] = {
1735			{
1736			.name = "main",
1737			.type = AVMEDIA_TYPE_VIDEO,
1738			.config_props = config_input_main,
1739			},
1740			{
1741			.name = "reference",
1742			.type = AVMEDIA_TYPE_VIDEO,
1743			.config_props = config_input_ref,
1744			},
1745			};
1746
1747			static const AVFilterPad ssim360_outputs[] = {
1748			{
1749			.name = "default",
1750			.type = AVMEDIA_TYPE_VIDEO,
1751			.config_props = config_output,
1752			},
1753			};
1754
1755			const FFFilter ff_vf_ssim360 = {
1756			.p.name = "ssim360",
1757			.p.description = NULL_IF_CONFIG_SMALL("Calculate the SSIM between two 360 video streams."),
1758			.p.priv_class = &ssim360_class,
1759			.preinit = ssim360_framesync_preinit,
1760			.init = init,
1761			.uninit = uninit,
1762			.activate = activate,
1763			.priv_size = sizeof(SSIM360Context),
1764			FILTER_INPUTS(ssim360_inputs),
1765			FILTER_OUTPUTS(ssim360_outputs),
1766			FILTER_PIXFMTS_ARRAY(ssim360_pixfmts),
1767			};
1768