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
22
27
33
34 #define MAX_THREADS 16
35
39
42
45
49
52
55
64
71
74
76 int w,
int h,
int n,
int plane,
float scale);
78
81
82 #define OFFSET(x) offsetof(ConvolveContext, x)
83 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
84
88 {
"first",
"process only first impulse, ignore rest", 0,
AV_OPT_TYPE_CONST, {.i64=0}, 0, 0,
FLAGS,
"impulse" },
92 };
93
113 };
114
116 {
121
123 s->planewidth[0] =
s->planewidth[3] =
w;
125 s->planeheight[0] =
s->planeheight[3] =
h;
126
127 s->nb_planes =
desc->nb_components;
128 s->depth =
desc->comp[0].depth;
129
130 for (
int i = 0;
i <
s->nb_planes;
i++) {
131 int w =
s->planewidth[
i];
132 int h =
s->planeheight[
i];
134
135 s->fft_len[
i] = 1 << (
av_log2(2 * n - 1));
136
139
142
145
148
151
154
157
160 }
161
162 return 0;
163 }
164
166 {
168
169 if (
ctx->inputs[0]->w !=
ctx->inputs[1]->w ||
170 ctx->inputs[0]->h !=
ctx->inputs[1]->h) {
173 }
174
175 return 0;
176 }
177
183
185 {
190 const int plane =
td->plane;
192 int start = (n * jobnr) / nb_jobs;
193 int end = (n * (jobnr+1)) / nb_jobs;
194 int y;
195
196 for (y = start; y < end; y++) {
197 s->tx_fn[plane](
s->fft[plane][jobnr], hdata_out + y * n, hdata_in + y * n,
sizeof(float));
198 }
199
200 return 0;
201 }
202
203 #define SQR(x) ((x) * (x))
204
208 int n,
int plane,
float scale)
209 {
210 float sum = 0.f;
212 int y, x;
213
215 for (y = 0; y <
h; y++) {
217
218 for (x = 0; x <
w; x++)
220 }
221
223 sum = 0.f;
224 for (y = 0; y <
h; y++) {
226
227 for (x = 0; x <
w; x++)
229 }
230
231 dev = sqrtf(sum / (
w *
h));
233 for (y = 0; y <
h; y++) {
235
236 for (x = 0; x <
w; x++) {
238 fft_hdata[y * n + x].
im = 0;
239 }
240
241 for (x =
w; x < n; x++) {
242 fft_hdata[y * n + x].
re = 0;
243 fft_hdata[y * n + x].
im = 0;
244 }
245 }
246
247 for (y =
h; y < n; y++) {
248 for (x = 0; x < n; x++) {
249 fft_hdata[y * n + x].
re = 0;
250 fft_hdata[y * n + x].
im = 0;
251 }
252 }
253 } else {
254 for (y = 0; y <
h; y++) {
255 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
256
257 for (x = 0; x <
w; x++)
259 }
260
262 sum = 0.f;
263 for (y = 0; y <
h; y++) {
264 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
265
266 for (x = 0; x <
w; x++)
268 }
269
270 dev = sqrtf(sum / (
w *
h));
272 for (y = 0; y <
h; y++) {
273 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
274
275 for (x = 0; x <
w; x++) {
277 fft_hdata[y * n + x].
im = 0;
278 }
279
280 for (x =
w; x < n; x++) {
281 fft_hdata[y * n + x].
re = 0;
282 fft_hdata[y * n + x].
im = 0;
283 }
284 }
285
286 for (y =
h; y < n; y++) {
287 for (x = 0; x < n; x++) {
288 fft_hdata[y * n + x].
re = 0;
289 fft_hdata[y * n + x].
im = 0;
290 }
291 }
292 }
293 }
294
297 {
298 const int iw = (n -
w) / 2, ih = (n -
h) / 2;
299 int y, x;
300
302 for (y = 0; y <
h; y++) {
304
305 for (x = 0; x <
w; x++) {
306 fft_hdata[(y + ih) * n + iw + x].
re =
src[x] *
scale;
307 fft_hdata[(y + ih) * n + iw + x].
im = 0;
308 }
309
310 for (x = 0; x < iw; x++) {
311 fft_hdata[(y + ih) * n + x].
re = fft_hdata[(y + ih) * n + iw].
re;
312 fft_hdata[(y + ih) * n + x].
im = 0;
313 }
314
315 for (x = n - iw; x < n; x++) {
316 fft_hdata[(y + ih) * n + x].
re = fft_hdata[(y + ih) * n + n - iw - 1].
re;
317 fft_hdata[(y + ih) * n + x].
im = 0;
318 }
319 }
320
321 for (y = 0; y < ih; y++) {
322 for (x = 0; x < n; x++) {
323 fft_hdata[y * n + x].
re = fft_hdata[ih * n + x].
re;
324 fft_hdata[y * n + x].
im = 0;
325 }
326 }
327
328 for (y = n - ih; y < n; y++) {
329 for (x = 0; x < n; x++) {
330 fft_hdata[y * n + x].
re = fft_hdata[(n - ih - 1) * n + x].
re;
331 fft_hdata[y * n + x].
im = 0;
332 }
333 }
334 } else {
335 for (y = 0; y <
h; y++) {
336 const uint16_t *
src = (
const uint16_t *)(in->
data[plane] + in->
linesize[plane] * y);
337
338 for (x = 0; x <
w; x++) {
339 fft_hdata[(y + ih) * n + iw + x].
re =
src[x] *
scale;
340 fft_hdata[(y + ih) * n + iw + x].
im = 0;
341 }
342
343 for (x = 0; x < iw; x++) {
344 fft_hdata[(y + ih) * n + x].
re = fft_hdata[(y + ih) * n + iw].
re;
345 fft_hdata[(y + ih) * n + x].
im = 0;
346 }
347
348 for (x = n - iw; x < n; x++) {
349 fft_hdata[(y + ih) * n + x].
re = fft_hdata[(y + ih) * n + n - iw - 1].
re;
350 fft_hdata[(y + ih) * n + x].
im = 0;
351 }
352 }
353
354 for (y = 0; y < ih; y++) {
355 for (x = 0; x < n; x++) {
356 fft_hdata[y * n + x].
re = fft_hdata[ih * n + x].
re;
357 fft_hdata[y * n + x].
im = 0;
358 }
359 }
360
361 for (y = n - ih; y < n; y++) {
362 for (x = 0; x < n; x++) {
363 fft_hdata[y * n + x].
re = fft_hdata[(n - ih - 1) * n + x].
re;
364 fft_hdata[y * n + x].
im = 0;
365 }
366 }
367 }
368 }
369
371 {
377 const int plane =
td->plane;
379 int start = (n * jobnr) / nb_jobs;
380 int end = (n * (jobnr+1)) / nb_jobs;
381 int y, x;
382
383 for (y = start; y < end; y++) {
384 for (x = 0; x < n; x++) {
385 vdata_in[y * n + x].
re = hdata[x * n + y].
re;
386 vdata_in[y * n + x].
im = hdata[x * n + y].
im;
387 }
388
389 s->tx_fn[plane](
s->fft[plane][jobnr], vdata_out + y * n, vdata_in + y * n,
sizeof(float));
390 }
391
392 return 0;
393 }
394
396 {
402 const int plane =
td->plane;
404 int start = (n * jobnr) / nb_jobs;
405 int end = (n * (jobnr+1)) / nb_jobs;
406 int y, x;
407
408 for (y = start; y < end; y++) {
409 s->itx_fn[plane](
s->ifft[plane][jobnr], vdata_out + y * n, vdata_in + y * n,
sizeof(float));
410
411 for (x = 0; x < n; x++) {
412 hdata[x * n + y].
re = vdata_out[y * n + x].
re;
413 hdata[x * n + y].
im = vdata_out[y * n + x].
im;
414 }
415 }
416
417 return 0;
418 }
419
421 {
426 const int plane =
td->plane;
428 int start = (n * jobnr) / nb_jobs;
429 int end = (n * (jobnr+1)) / nb_jobs;
430 int y;
431
432 for (y = start; y < end; y++) {
433 s->itx_fn[plane](
s->ifft[plane][jobnr], hdata_out + y * n, hdata_in + y * n,
sizeof(float));
434 }
435
436 return 0;
437 }
438
440 int w,
int h,
int n,
int plane,
float scale)
441 {
442 const int imax = (1 <<
s->depth) - 1;
443
446 for (
int y = 0; y <
h; y++) {
447 uint8_t *dst =
out->data[plane] + y *
out->linesize[plane];
448 for (
int x = 0; x <
w; x++)
450 }
451 } else {
452 for (
int y = 0; y <
h; y++) {
453 uint16_t *dst = (uint16_t *)(
out->data[plane] + y *
out->linesize[plane]);
454 for (
int x = 0; x <
w; x++)
456 }
457 }
458 }
459
461 int w,
int h,
int n,
int plane,
float scale)
462 {
463 const int max = (1 <<
s->depth) - 1;
464 const int hh =
h / 2;
465 const int hw =
w / 2;
466 int y, x;
467
469 for (y = 0; y < hh; y++) {
470 uint8_t *dst =
out->data[plane] + (y + hh) *
out->linesize[plane] + hw;
471 for (x = 0; x < hw; x++)
473 }
474 for (y = 0; y < hh; y++) {
475 uint8_t *dst =
out->data[plane] + (y + hh) *
out->linesize[plane];
476 for (x = 0; x < hw; x++)
478 }
479 for (y = 0; y < hh; y++) {
480 uint8_t *dst =
out->data[plane] + y *
out->linesize[plane] + hw;
481 for (x = 0; x < hw; x++)
483 }
484 for (y = 0; y < hh; y++) {
485 uint8_t *dst =
out->data[plane] + y *
out->linesize[plane];
486 for (x = 0; x < hw; x++)
488 }
489 } else {
490 for (y = 0; y < hh; y++) {
491 uint16_t *dst = (uint16_t *)(
out->data[plane] + (y + hh) *
out->linesize[plane] + hw * 2);
492 for (x = 0; x < hw; x++)
494 }
495 for (y = 0; y < hh; y++) {
496 uint16_t *dst = (uint16_t *)(
out->data[plane] + (y + hh) *
out->linesize[plane]);
497 for (x = 0; x < hw; x++)
499 }
500 for (y = 0; y < hh; y++) {
501 uint16_t *dst = (uint16_t *)(
out->data[plane] + y *
out->linesize[plane] + hw * 2);
502 for (x = 0; x < hw; x++)
504 }
505 for (y = 0; y < hh; y++) {
506 uint16_t *dst = (uint16_t *)(
out->data[plane] + y *
out->linesize[plane]);
507 for (x = 0; x < hw; x++)
509 }
510 }
511 }
512
514 {
519 const float noise =
s->noise;
521 int start = (n * jobnr) / nb_jobs;
522 int end = (n * (jobnr+1)) / nb_jobs;
523 int y, x;
524
525 for (y = start; y < end; y++) {
526 int yn = y * n;
527
528 for (x = 0; x < n; x++) {
529 float re,
im, ire, iim;
530
535
538 }
539 }
540
541 return 0;
542 }
543
545 {
550 const float scale = 1.f / (n * n);
551 int start = (n * jobnr) / nb_jobs;
552 int end = (n * (jobnr+1)) / nb_jobs;
553
554 for (int y = start; y < end; y++) {
555 int yn = y * n;
556
557 for (int x = 0; x < n; x++) {
558 float re,
im, ire, iim;
559
564
567 }
568 }
569
570 return 0;
571 }
572
574 {
579 const float noise =
s->noise;
581 int start = (n * jobnr) / nb_jobs;
582 int end = (n * (jobnr+1)) / nb_jobs;
583 int y, x;
584
585 for (y = start; y < end; y++) {
586 int yn = y * n;
587
588 for (x = 0; x < n; x++) {
589 float re,
im, ire, iim, div;
590
595 div = ire * ire + iim * iim +
noise;
596
597 input[yn + x].re = (ire *
re + iim *
im) / div;
598 input[yn + x].im = (ire *
im - iim *
re) / div;
599 }
600 }
601
602 return 0;
603 }
604
606 {
608 const int n =
s->fft_len[plane];
609 const int w =
s->secondarywidth[plane];
610 const int h =
s->secondaryheight[plane];
612 float total = 0;
613
615 for (
int y = 0; y <
h; y++) {
616 const uint8_t *
src = (
const uint8_t *)(impulsepic->
data[plane] + y * impulsepic->
linesize[plane]) ;
617 for (
int x = 0; x <
w; x++) {
619 }
620 }
621 } else {
622 for (
int y = 0; y <
h; y++) {
623 const uint16_t *
src = (
const uint16_t *)(impulsepic->
data[plane] + y * impulsepic->
linesize[plane]) ;
624 for (
int x = 0; x <
w; x++) {
626 }
627 }
628 }
629 total =
FFMAX(1, total);
630
631 s->get_input(
s,
s->fft_hdata_impulse_in[plane], impulsepic,
w,
h, n, plane, 1.f / total);
632
635 td.hdata_in =
s->fft_hdata_impulse_in[plane];
636 td.vdata_in =
s->fft_vdata_impulse_in[plane];
637 td.hdata_out =
s->fft_hdata_impulse_out[plane];
638 td.vdata_out =
s->fft_vdata_impulse_out[plane];
639
644
645 s->got_impulse[plane] = 1;
646 }
647
649 {
651 const int n =
s->fft_len[plane];
653
654 s->get_input(
s,
s->fft_hdata_impulse_in[plane], secondary,
655 s->secondarywidth[plane],
656 s->secondaryheight[plane],
657 n, plane, 1.f);
658
661 td.hdata_in =
s->fft_hdata_impulse_in[plane];
662 td.vdata_in =
s->fft_vdata_impulse_in[plane];
663 td.hdata_out =
s->fft_hdata_impulse_out[plane];
664 td.vdata_out =
s->fft_vdata_impulse_out[plane];
665
670
671 s->got_impulse[plane] = 1;
672 }
673
675 {
681
685 if (!impulsepic)
687
688 for (plane = 0; plane <
s->nb_planes; plane++) {
691 const int n =
s->fft_len[plane];
692 const int w =
s->primarywidth[plane];
693 const int h =
s->primaryheight[plane];
694 const int ow =
s->planewidth[plane];
695 const int oh =
s->planeheight[plane];
697
698 if (!(
s->planes & (1 << plane))) {
699 continue;
700 }
701
702 td.plane = plane,
td.n = n;
703 s->get_input(
s,
s->fft_hdata_in[plane], mainpic,
w,
h, n, plane, 1.f);
704
705 td.hdata_in =
s->fft_hdata_in[plane];
706 td.vdata_in =
s->fft_vdata_in[plane];
707 td.hdata_out =
s->fft_hdata_out[plane];
708 td.vdata_out =
s->fft_vdata_out[plane];
709
714
715 if ((!
s->impulse && !
s->got_impulse[plane]) ||
s->impulse) {
716 s->prepare_impulse(
ctx, impulsepic, plane);
717 }
718
721
724
725 td.hdata_in =
s->fft_hdata_out[plane];
726 td.vdata_in =
s->fft_vdata_out[plane];
727 td.hdata_out =
s->fft_hdata_in[plane];
728 td.vdata_out =
s->fft_vdata_in[plane];
729
732
733 td.hdata_out =
s->fft_hdata_out[plane];
734 td.hdata_in =
s->fft_hdata_in[plane];
735
738
739 s->get_output(
s,
s->fft_hdata_out[plane], mainpic, ow, oh, n, plane, 1.f / (n * n));
740 }
741
743 }
744
746 {
753
755 s->primarywidth[0] =
s->primarywidth[3] = mainlink->
w;
757 s->primaryheight[0] =
s->primaryheight[3] = mainlink->
h;
758
760 s->secondarywidth[0] =
s->secondarywidth[3] = secondlink->
w;
762 s->secondaryheight[0] =
s->secondaryheight[3] = secondlink->
h;
763
768 outlink->
w = mainlink->
w;
769 outlink->
h = mainlink->
h;
773
776
777 for (
i = 0;
i <
s->nb_planes;
i++) {
780
787 }
788 }
789
790 return 0;
791 }
792
794 {
797 }
798
800 {
802
803 if (!strcmp(
ctx->filter->name,
"convolve")) {
808 }
else if (!strcmp(
ctx->filter->name,
"xcorrelate")) {
813 }
else if (!strcmp(
ctx->filter->name,
"deconvolve")) {
818 } else {
820 }
821
822 return 0;
823 }
824
826 {
829
830 for (
i = 0;
i < 4;
i++) {
839
843 }
844 }
845
847 }
848
850 {
854 },{
855 .name = "impulse",
858 },
859 };
860
862 {
866 },
867 };
868
870
871 #if CONFIG_CONVOLVE_FILTER
872
874
878 .preinit = convolve_framesync_preinit,
883 .priv_class = &convolve_class,
888 };
889
890 #endif /* CONFIG_CONVOLVE_FILTER */
891
892 #if CONFIG_DECONVOLVE_FILTER
893
894 static const AVOption deconvolve_options[] = {
897 {
"first",
"process only first impulse, ignore rest", 0,
AV_OPT_TYPE_CONST, {.i64=0}, 0, 0,
FLAGS,
"impulse" },
901 };
902
904
906 .
name =
"deconvolve",
908 .preinit = convolve_framesync_preinit,
913 .priv_class = &deconvolve_class,
918 };
919
920 #endif /* CONFIG_DECONVOLVE_FILTER */
921
922 #if CONFIG_XCORRELATE_FILTER
923
924 static const AVOption xcorrelate_options[] = {
927 {
"first",
"process only first secondary frame, ignore rest", 0,
AV_OPT_TYPE_CONST, {.i64=0}, 0, 0,
FLAGS,
"impulse" },
930 };
931
933
935 {
937
938 if (
ctx->inputs[0]->w <=
ctx->inputs[1]->w ||
939 ctx->inputs[0]->h <=
ctx->inputs[1]->h) {
940 av_log(
ctx,
AV_LOG_ERROR,
"Width and height of second input videos must be less than first input.\n");
942 }
943
944 return 0;
945 }
946
948 {
952 },{
953 .name = "secondary",
955 .config_props = config_input_secondary,
956 },
957 };
958
960 {
964 },
965 };
966
968 .
name =
"xcorrelate",
969 .description =
NULL_IF_CONFIG_SMALL(
"Cross-correlate first video stream with second video stream."),
970 .preinit = convolve_framesync_preinit,
975 .priv_class = &xcorrelate_class,
980 };
981
982 #endif /* CONFIG_XCORRELATE_FILTER */