laplace_filter.c (@Copyright (C) 2004 University of Texas at Austin: 2-D Laplacian operator)

/* 2-D finite-difference Laplacian operation. */

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<string.h>

void laplac2_lop(int adj, int nz, int nx, float *in, float *out)
/*< linear operator >*/
{
    int iz,ix,j;

    for (ix=0; ix < nx; ix++) {
	for (iz=0; iz < nz; iz++) {
	    j = iz+ix*nz;

	    if (iz > 0) {
		if (adj) {
		    out[j-1] -= in[j];
		    out[j]   += in[j];
		} else {
		    out[j] += in[j] - in[j-1];
		}
	    }
	    if (iz < nz-1) {
		if (adj) {
		    out[j+1] -= in[j];
		    out[j]   += in[j];
		} else {
		    out[j] += in[j] - in[j+1];
		}
	    }

	    if (ix > 0) {
		if (adj) {
		    out[j-nz] -= in[j];
		    out[j]    += in[j];
		} else {
		    out[j] += in[j] - in[j-nz];
		}
	    }
	    if (ix < nx-1) {
		if (adj) {
		    out[j+nz] -= in[j];
		    out[j]    += in[j];
		} else {
		    out[j] += in[j] - in[j+nz];
		}
	    }
	}
    }
}

int main(int argc, char* argv[])
{
    int nz,nx,n3,i3;
    float *in,*out;
    int adj;
    char FN1[250]={"mig_ns.dat"};
    char FN2[250]={"mig_ns_laplace.dat"};

    FILE *fp1, *fp2;
    fp1  = fopen (FN1,"rb");
    fp2 = fopen (FN2,"wb");

    nz=301;
    nx=601;
    n3 = 1;

    adj=1; /* adjoint flag */

     in = (float*)malloc(nx*nz*sizeof(float));
    out = (float*)malloc(nx*nz*sizeof(float));

    fread(in,sizeof(float),nx*nz,fp1);
    for (i3=0; i3<n3; i3++) {
       laplac2_lop ( adj, nz, nx, in,  out );
    }
    fwrite(out,sizeof(float),nx*nz,fp2);
    exit(0);



}