suseislet/allp2.c at Main · vandanue/suseislet · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
/*
 * allp2.c  --  2D all-pass plane-wave destruction filter
 */

#include <stdlib.h>     /* malloc, free                               */
#include <stdbool.h>    /* bool, true, false                          */

#include "su.h"         /* ealloc1float, err                          */
#include "suutils.h"    /* adjnull, SF_NINT                           */
#include "allp2.h"
#include "apfilt.h"     /* apfilt_init, apfilt_close, passfilter, aderfilter */

/* Full struct definition (private to this file) */
struct Allpass2 {
    int    nx, ny, nw, nj;
    bool   drift;
    float *flt;    /* filter coefficients buffer [2*nw+1] */
    float **pp;    /* dip field [ny][nx] */
};

static allpas2 ap2;    /* global for allpass21_lop / allpass22_init */


/* ================================================================== */
allpas2 allpass2_init(int nw, int nj, int nx, int ny,
                      bool drift, float **pp)
/*< initialize the 2D all-pass filter >*/
{
    allpas2 ap;

    ap = (allpas2) malloc(sizeof(*ap));   /* was: sf_alloc(1,sizeof(*ap)) */
    if (!ap) err("allpass2_init: malloc failed");

    ap->nw    = nw;
    ap->nj    = nj;
    ap->nx    = nx;
    ap->ny    = ny;
    ap->drift = drift;
    ap->pp    = pp;
    ap->flt   = ealloc1float(2*nw+1);    /* was: sf_floatalloc(2*nw+1)   */

    apfilt_init(nw);    /* from apfilt.c */
    return ap;
}


/* ================================================================== */
void allpass2_close(allpas2 ap)
/*< free allocated storage >*/
{
    apfilt_close();     /* from apfilt.c */
    free(ap->flt);
    free(ap);
}


/* ================================================================== */
void allpass22_init(allpas2 ap1)
/*< set global pointer for allpass21_lop >*/
{
    ap2 = ap1;
}


/* ================================================================== */
void allpass21_lop(bool adj, bool add, int n1, int n2,
                   float *xx, float *yy)
/*< PWD as a linear operator (uses global ap2 set by allpass22_init) >*/
{
    int   i, ix, iy, iw, is, nx, ny, id;

    adjnull(adj, add, n1, n2, xx, yy);   /* was: sf_adjnull */

    nx = ap2->nx;
    ny = ap2->ny;

    for (iy = 0; iy < ny-1; iy++) {
        for (ix = ap2->nw*ap2->nj; ix < nx-ap2->nw*ap2->nj; ix++) {
            i = ix + iy*nx;

            if (ap2->drift) {
                id = SF_NINT(ap2->pp[iy][ix]);
                if (ix - ap2->nw*ap2->nj - id < 0 ||
                    ix + ap2->nw*ap2->nj - id >= nx) continue;

                passfilter(ap2->pp[iy][ix]-id, ap2->flt);

                for (iw = 0; iw <= 2*ap2->nw; iw++) {
                    is = (iw - ap2->nw) * ap2->nj;
                    if (adj) {
                        xx[i+is+nx] += yy[i] * ap2->flt[iw];
                        xx[i-is-id] -= yy[i] * ap2->flt[iw];
                    } else {
                        yy[i] += (xx[i+is+nx] - xx[i-is-id]) * ap2->flt[iw];
                    }
                }
            } else {
                passfilter(ap2->pp[iy][ix], ap2->flt);

                for (iw = 0; iw <= 2*ap2->nw; iw++) {
                    is = (iw - ap2->nw) * ap2->nj;
                    if (adj) {
                        xx[i+is+nx] += yy[i] * ap2->flt[iw];
                        xx[i-is]    -= yy[i] * ap2->flt[iw];
                    } else {
                        yy[i] += (xx[i+is+nx] - xx[i-is]) * ap2->flt[iw];
                    }
                }
            }
        }
    }
}


/* ================================================================== */
void allpass21(bool der, const allpas2 ap, float **xx, float **yy)
/*< 2D plane-wave destruction (pure algorithm, unchanged from Madagascar) >*/
{
    int ix, iy, iw, is, id;

    for (iy = 0; iy < ap->ny; iy++)
        for (ix = 0; ix < ap->nx; ix++)
            yy[iy][ix] = 0.0f;

    for (iy = 0; iy < ap->ny-1; iy++) {
        for (ix = ap->nw*ap->nj; ix < ap->nx-ap->nw*ap->nj; ix++) {
            if (ap->drift) {
                id = SF_NINT(ap->pp[iy][ix]);
                if (ix - ap->nw*ap->nj - id < 0 ||
                    ix + ap->nw*ap->nj - id >= ap->nx) continue;

                if (der) aderfilter(ap->pp[iy][ix]-id, ap->flt);
                else     passfilter(ap->pp[iy][ix]-id, ap->flt);

                for (iw = 0; iw <= 2*ap->nw; iw++) {
                    is = (iw - ap->nw) * ap->nj;
                    yy[iy][ix] += (xx[iy+1][ix+is] -
                                   xx[iy  ][ix-is-id]) * ap->flt[iw];
                }
            } else {
                if (der) aderfilter(ap->pp[iy][ix], ap->flt);
                else     passfilter(ap->pp[iy][ix], ap->flt);

                for (iw = 0; iw <= 2*ap->nw; iw++) {
                    is = (iw - ap->nw) * ap->nj;
                    yy[iy][ix] += (xx[iy+1][ix+is] -
                                   xx[iy  ][ix-is]) * ap->flt[iw];
                }
            }
        }
    }
}