1 /*
2 * Copyright (c) 2017 Paul B Mahol
3 *
4 * This file is part of FFmpeg.
5 *
6 * FFmpeg is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * FFmpeg is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with FFmpeg; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21 #include "config_components.h"
22
24
28
32
33 #define MAX_THREADS 16
34
38
41
44
48
51
54
63
70
73
75 int w,
int h,
int n,
int plane,
float scale);
77
80
81 #define OFFSET(x) offsetof(ConvolveContext, x)
82 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
83
87 {
"first",
"process only first impulse, ignore rest", 0,
AV_OPT_TYPE_CONST, {.i64=0}, 0, 0,
FLAGS, .unit =
"impulse" },
91 };
92
112 };
113
115 {
120
122 s->planewidth[0] =
s->planewidth[3] =
w;
124 s->planeheight[0] =
s->planeheight[3] =
h;
125
126 s->nb_planes =
desc->nb_components;
127 s->depth =
desc->comp[0].depth;
128
129 for (
int i = 0;
i <
s->nb_planes;
i++) {
130 int w =
s->planewidth[
i];
131 int h =
s->planeheight[
i];
133
134 s->fft_len[
i] = 1 << (
av_log2(2 * n - 1));
135
138
141
144
147
150
153
156
159 }
160
161 return 0;
162 }
163
165 {
167
168 if (
ctx->inputs[0]->w !=
ctx->inputs[1]->w ||
169 ctx->inputs[0]->h !=
ctx->inputs[1]->h) {
172 }
173
174 return 0;
175 }
176
182
184 {
189 const int plane =
td->plane;
191 int start = (n * jobnr) / nb_jobs;
192 int end = (n * (jobnr+1)) / nb_jobs;
193 int y;
194
195 for (y = start; y < end; y++) {
196 s->tx_fn[plane](
s->fft[plane][jobnr], hdata_out + y * n, hdata_in + y * n,
sizeof(
AVComplexFloat));
197 }
198
199 return 0;
200 }
201
202 #define SQR(x) ((x) * (x))
203
207 int n,
int plane,
float scale)
208 {
209 float sum = 0.f;
211 int y, x;
212
214 for (y = 0; y <
h; y++) {
216
217 for (x = 0; x <
w; x++)
219 }
220
222 sum = 0.f;
223 for (y = 0; y <
h; y++) {
225
226 for (x = 0; x <
w; x++)
228 }
229
232 for (y = 0; y <
h; y++) {
234
235 for (x = 0; x <
w; x++) {
237 fft_hdata[y * n + x].
im = 0;
238 }
239
240 for (x =
w; x < n; x++) {
241 fft_hdata[y * n + x].
re = 0;
242 fft_hdata[y * n + x].
im = 0;
243 }
244 }
245
246 for (y =
h; y < n; y++) {
247 for (x = 0; x < n; x++) {
248 fft_hdata[y * n + x].
re = 0;
249 fft_hdata[y * n + x].
im = 0;
250 }
251 }
252 } else {
253 for (y = 0; y <
h; y++) {
254 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
255
256 for (x = 0; x <
w; x++)
258 }
259
261 sum = 0.f;
262 for (y = 0; y <
h; y++) {
263 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
264
265 for (x = 0; x <
w; x++)
267 }
268
271 for (y = 0; y <
h; y++) {
272 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
273
274 for (x = 0; x <
w; x++) {
276 fft_hdata[y * n + x].
im = 0;
277 }
278
279 for (x =
w; x < n; x++) {
280 fft_hdata[y * n + x].
re = 0;
281 fft_hdata[y * n + x].
im = 0;
282 }
283 }
284
285 for (y =
h; y < n; y++) {
286 for (x = 0; x < n; x++) {
287 fft_hdata[y * n + x].
re = 0;
288 fft_hdata[y * n + x].
im = 0;
289 }
290 }
291 }
292 }
293
296 {
297 const int iw = (n -
w) / 2, ih = (n -
h) / 2;
298 int y, x;
299
301 for (y = 0; y <
h; y++) {
303
304 for (x = 0; x <
w; x++) {
305 fft_hdata[(y + ih) * n + iw + x].re =
src[x] *
scale;
306 fft_hdata[(y + ih) * n + iw + x].im = 0;
307 }
308
309 for (x = 0; x < iw; x++) {
310 fft_hdata[(y + ih) * n + x].re = fft_hdata[(y + ih) * n + iw].
re;
311 fft_hdata[(y + ih) * n + x].im = 0;
312 }
313
314 for (x = n - iw; x < n; x++) {
315 fft_hdata[(y + ih) * n + x].re = fft_hdata[(y + ih) * n + n - iw - 1].
re;
316 fft_hdata[(y + ih) * n + x].im = 0;
317 }
318 }
319
320 for (y = 0; y < ih; y++) {
321 for (x = 0; x < n; x++) {
322 fft_hdata[y * n + x].
re = fft_hdata[ih * n + x].
re;
323 fft_hdata[y * n + x].
im = 0;
324 }
325 }
326
327 for (y = n - ih; y < n; y++) {
328 for (x = 0; x < n; x++) {
329 fft_hdata[y * n + x].
re = fft_hdata[(n - ih - 1) * n + x].re;
330 fft_hdata[y * n + x].
im = 0;
331 }
332 }
333 } else {
334 for (y = 0; y <
h; y++) {
335 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
336
337 for (x = 0; x <
w; x++) {
338 fft_hdata[(y + ih) * n + iw + x].re =
src[x] *
scale;
339 fft_hdata[(y + ih) * n + iw + x].im = 0;
340 }
341
342 for (x = 0; x < iw; x++) {
343 fft_hdata[(y + ih) * n + x].re = fft_hdata[(y + ih) * n + iw].
re;
344 fft_hdata[(y + ih) * n + x].im = 0;
345 }
346
347 for (x = n - iw; x < n; x++) {
348 fft_hdata[(y + ih) * n + x].re = fft_hdata[(y + ih) * n + n - iw - 1].
re;
349 fft_hdata[(y + ih) * n + x].im = 0;
350 }
351 }
352
353 for (y = 0; y < ih; y++) {
354 for (x = 0; x < n; x++) {
355 fft_hdata[y * n + x].
re = fft_hdata[ih * n + x].
re;
356 fft_hdata[y * n + x].
im = 0;
357 }
358 }
359
360 for (y = n - ih; y < n; y++) {
361 for (x = 0; x < n; x++) {
362 fft_hdata[y * n + x].
re = fft_hdata[(n - ih - 1) * n + x].re;
363 fft_hdata[y * n + x].
im = 0;
364 }
365 }
366 }
367 }
368
370 {
376 const int plane =
td->plane;
378 int start = (n * jobnr) / nb_jobs;
379 int end = (n * (jobnr+1)) / nb_jobs;
380 int y, x;
381
382 for (y = start; y < end; y++) {
383 for (x = 0; x < n; x++) {
384 vdata_in[y * n + x].
re = hdata[x * n + y].
re;
385 vdata_in[y * n + x].
im = hdata[x * n + y].
im;
386 }
387
388 s->tx_fn[plane](
s->fft[plane][jobnr], vdata_out + y * n, vdata_in + y * n,
sizeof(
AVComplexFloat));
389 }
390
391 return 0;
392 }
393
395 {
401 const int plane =
td->plane;
403 int start = (n * jobnr) / nb_jobs;
404 int end = (n * (jobnr+1)) / nb_jobs;
405 int y, x;
406
407 for (y = start; y < end; y++) {
408 s->itx_fn[plane](
s->ifft[plane][jobnr], vdata_out + y * n, vdata_in + y * n,
sizeof(
AVComplexFloat));
409
410 for (x = 0; x < n; x++) {
411 hdata[x * n + y].
re = vdata_out[y * n + x].
re;
412 hdata[x * n + y].
im = vdata_out[y * n + x].
im;
413 }
414 }
415
416 return 0;
417 }
418
420 {
425 const int plane =
td->plane;
427 int start = (n * jobnr) / nb_jobs;
428 int end = (n * (jobnr+1)) / nb_jobs;
429 int y;
430
431 for (y = start; y < end; y++) {
432 s->itx_fn[plane](
s->ifft[plane][jobnr], hdata_out + y * n, hdata_in + y * n,
sizeof(
AVComplexFloat));
433 }
434
435 return 0;
436 }
437
439 int w,
int h,
int n,
int plane,
float scale)
440 {
441 const int imax = (1 <<
s->depth) - 1;
442
445 for (
int y = 0; y <
h; y++) {
446 uint8_t *dst =
out->data[plane] + y *
out->linesize[plane];
447 for (
int x = 0; x <
w; x++)
449 }
450 } else {
451 for (
int y = 0; y <
h; y++) {
452 uint16_t *dst = (uint16_t *)(
out->data[plane] + y *
out->linesize[plane]);
453 for (
int x = 0; x <
w; x++)
455 }
456 }
457 }
458
460 int w,
int h,
int n,
int plane,
float scale)
461 {
462 const int max = (1 <<
s->depth) - 1;
463 const int hh =
h / 2;
464 const int hw =
w / 2;
465 int y, x;
466
468 for (y = 0; y < hh; y++) {
469 uint8_t *dst =
out->data[plane] + (y + hh) *
out->linesize[plane] + hw;
470 for (x = 0; x < hw; x++)
472 }
473 for (y = 0; y < hh; y++) {
474 uint8_t *dst =
out->data[plane] + (y + hh) *
out->linesize[plane];
475 for (x = 0; x < hw; x++)
477 }
478 for (y = 0; y < hh; y++) {
479 uint8_t *dst =
out->data[plane] + y *
out->linesize[plane] + hw;
480 for (x = 0; x < hw; x++)
482 }
483 for (y = 0; y < hh; y++) {
484 uint8_t *dst =
out->data[plane] + y *
out->linesize[plane];
485 for (x = 0; x < hw; x++)
487 }
488 } else {
489 for (y = 0; y < hh; y++) {
490 uint16_t *dst = (uint16_t *)(
out->data[plane] + (y + hh) *
out->linesize[plane] + hw * 2);
491 for (x = 0; x < hw; x++)
493 }
494 for (y = 0; y < hh; y++) {
495 uint16_t *dst = (uint16_t *)(
out->data[plane] + (y + hh) *
out->linesize[plane]);
496 for (x = 0; x < hw; x++)
498 }
499 for (y = 0; y < hh; y++) {
500 uint16_t *dst = (uint16_t *)(
out->data[plane] + y *
out->linesize[plane] + hw * 2);
501 for (x = 0; x < hw; x++)
503 }
504 for (y = 0; y < hh; y++) {
505 uint16_t *dst = (uint16_t *)(
out->data[plane] + y *
out->linesize[plane]);
506 for (x = 0; x < hw; x++)
508 }
509 }
510 }
511
513 {
518 const float noise =
s->noise;
520 int start = (n * jobnr) / nb_jobs;
521 int end = (n * (jobnr+1)) / nb_jobs;
522 int y, x;
523
524 for (y = start; y < end; y++) {
525 int yn = y * n;
526
527 for (x = 0; x < n; x++) {
528 float re, im, ire, iim;
529
530 re =
input[yn + x].re;
531 im =
input[yn + x].im;
534
535 input[yn + x].re = ire * re - iim * im;
536 input[yn + x].im = iim * re + ire * im;
537 }
538 }
539
540 return 0;
541 }
542
544 {
549 const float scale = 1.f / (n * n);
550 int start = (n * jobnr) / nb_jobs;
551 int end = (n * (jobnr+1)) / nb_jobs;
552
553 for (int y = start; y < end; y++) {
554 int yn = y * n;
555
556 for (int x = 0; x < n; x++) {
557 float re, im, ire, iim;
558
559 re =
input[yn + x].re;
560 im =
input[yn + x].im;
563
564 input[yn + x].re = ire * re - iim * im;
565 input[yn + x].im = iim * re + ire * im;
566 }
567 }
568
569 return 0;
570 }
571
573 {
578 const float noise =
s->noise;
580 int start = (n * jobnr) / nb_jobs;
581 int end = (n * (jobnr+1)) / nb_jobs;
582 int y, x;
583
584 for (y = start; y < end; y++) {
585 int yn = y * n;
586
587 for (x = 0; x < n; x++) {
588 float re, im, ire, iim, div;
589
590 re =
input[yn + x].re;
591 im =
input[yn + x].im;
594 div = ire * ire + iim * iim +
noise;
595
596 input[yn + x].re = (ire * re + iim * im) / div;
597 input[yn + x].im = (ire * im - iim * re) / div;
598 }
599 }
600
601 return 0;
602 }
603
605 {
607 const int n =
s->fft_len[plane];
608 const int w =
s->secondarywidth[plane];
609 const int h =
s->secondaryheight[plane];
611 float total = 0;
612
614 for (
int y = 0; y <
h; y++) {
615 const uint8_t *
src = (
const uint8_t *)(impulsepic->
data[plane] + y * impulsepic->
linesize[plane]) ;
616 for (
int x = 0; x <
w; x++) {
618 }
619 }
620 } else {
621 for (
int y = 0; y <
h; y++) {
622 const uint16_t *
src = (
const uint16_t *)(impulsepic->
data[plane] + y * impulsepic->
linesize[plane]) ;
623 for (
int x = 0; x <
w; x++) {
625 }
626 }
627 }
628 total =
FFMAX(1, total);
629
630 s->get_input(
s,
s->fft_hdata_impulse_in[plane], impulsepic,
w,
h, n, plane, 1.f / total);
631
634 td.hdata_in =
s->fft_hdata_impulse_in[plane];
635 td.vdata_in =
s->fft_vdata_impulse_in[plane];
636 td.hdata_out =
s->fft_hdata_impulse_out[plane];
637 td.vdata_out =
s->fft_vdata_impulse_out[plane];
638
643
644 s->got_impulse[plane] = 1;
645 }
646
648 {
650 const int n =
s->fft_len[plane];
652
653 s->get_input(
s,
s->fft_hdata_impulse_in[plane], secondary,
654 s->secondarywidth[plane],
655 s->secondaryheight[plane],
656 n, plane, 1.f);
657
660 td.hdata_in =
s->fft_hdata_impulse_in[plane];
661 td.vdata_in =
s->fft_vdata_impulse_in[plane];
662 td.hdata_out =
s->fft_hdata_impulse_out[plane];
663 td.vdata_out =
s->fft_vdata_impulse_out[plane];
664
669
670 s->got_impulse[plane] = 1;
671 }
672
674 {
680
684 if (!impulsepic)
686
687 for (plane = 0; plane <
s->nb_planes; plane++) {
690 const int n =
s->fft_len[plane];
691 const int w =
s->primarywidth[plane];
692 const int h =
s->primaryheight[plane];
693 const int ow =
s->planewidth[plane];
694 const int oh =
s->planeheight[plane];
696
697 if (!(
s->planes & (1 << plane))) {
698 continue;
699 }
700
701 td.plane = plane,
td.n = n;
702 s->get_input(
s,
s->fft_hdata_in[plane], mainpic,
w,
h, n, plane, 1.f);
703
704 td.hdata_in =
s->fft_hdata_in[plane];
705 td.vdata_in =
s->fft_vdata_in[plane];
706 td.hdata_out =
s->fft_hdata_out[plane];
707 td.vdata_out =
s->fft_vdata_out[plane];
708
713
714 if ((!
s->impulse && !
s->got_impulse[plane]) ||
s->impulse) {
715 s->prepare_impulse(
ctx, impulsepic, plane);
716 }
717
720
723
724 td.hdata_in =
s->fft_hdata_out[plane];
725 td.vdata_in =
s->fft_vdata_out[plane];
726 td.hdata_out =
s->fft_hdata_in[plane];
727 td.vdata_out =
s->fft_vdata_in[plane];
728
731
732 td.hdata_out =
s->fft_hdata_out[plane];
733 td.hdata_in =
s->fft_hdata_in[plane];
734
737
738 s->get_output(
s,
s->fft_hdata_out[plane], mainpic, ow, oh, n, plane, 1.f / (n * n));
739 }
740
742 }
743
745 {
752
754 s->primarywidth[0] =
s->primarywidth[3] = mainlink->
w;
756 s->primaryheight[0] =
s->primaryheight[3] = mainlink->
h;
757
759 s->secondarywidth[0] =
s->secondarywidth[3] = secondlink->
w;
761 s->secondaryheight[0] =
s->secondaryheight[3] = secondlink->
h;
762
767 outlink->
w = mainlink->
w;
768 outlink->
h = mainlink->
h;
772
775
776 for (
i = 0;
i <
s->nb_planes;
i++) {
779
786 }
787 }
788
789 return 0;
790 }
791
793 {
796 }
797
799 {
801
802 if (!strcmp(
ctx->filter->name,
"convolve")) {
807 }
else if (!strcmp(
ctx->filter->name,
"xcorrelate")) {
812 }
else if (!strcmp(
ctx->filter->name,
"deconvolve")) {
817 } else {
819 }
820
821 return 0;
822 }
823
825 {
828
829 for (
i = 0;
i < 4;
i++) {
838
842 }
843 }
844
846 }
847
849 {
853 },{
854 .name = "impulse",
857 },
858 };
859
861 {
865 },
866 };
867
869
870 #if CONFIG_CONVOLVE_FILTER
871
873
877 .preinit = convolve_framesync_preinit,
882 .priv_class = &convolve_class,
887 };
888
889 #endif /* CONFIG_CONVOLVE_FILTER */
890
891 #if CONFIG_DECONVOLVE_FILTER
892
893 static const AVOption deconvolve_options[] = {
896 {
"first",
"process only first impulse, ignore rest", 0,
AV_OPT_TYPE_CONST, {.i64=0}, 0, 0,
FLAGS, .unit =
"impulse" },
900 };
901
903
905 .
name =
"deconvolve",
907 .preinit = convolve_framesync_preinit,
912 .priv_class = &deconvolve_class,
917 };
918
919 #endif /* CONFIG_DECONVOLVE_FILTER */
920
921 #if CONFIG_XCORRELATE_FILTER
922
923 static const AVOption xcorrelate_options[] = {
925 {
"secondary",
"when to process secondary frame",
OFFSET(impulse),
AV_OPT_TYPE_INT, {.i64=1}, 0, 1,
FLAGS, .unit =
"impulse" },
926 {
"first",
"process only first secondary frame, ignore rest", 0,
AV_OPT_TYPE_CONST, {.i64=0}, 0, 0,
FLAGS, .unit =
"impulse" },
927 {
"all",
"process all secondary frames", 0,
AV_OPT_TYPE_CONST, {.i64=1}, 0, 0,
FLAGS, .unit =
"impulse" },
929 };
930
932
934 {
936
937 if (
ctx->inputs[0]->w <=
ctx->inputs[1]->w ||
938 ctx->inputs[0]->h <=
ctx->inputs[1]->h) {
939 av_log(
ctx,
AV_LOG_ERROR,
"Width and height of second input videos must be less than first input.\n");
941 }
942
943 return 0;
944 }
945
947 {
951 },{
952 .name = "secondary",
954 .config_props = config_input_secondary,
955 },
956 };
957
958 #define xcorrelate_outputs convolve_outputs
959
961 .
name =
"xcorrelate",
962 .description =
NULL_IF_CONFIG_SMALL(
"Cross-correlate first video stream with second video stream."),
963 .preinit = convolve_framesync_preinit,
968 .priv_class = &xcorrelate_class,
973 };
974
975 #endif /* CONFIG_XCORRELATE_FILTER */