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);}