smooth.c 6.38 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
/*-----------------------------------------------------------------------------------------
Florian Wittkamp's avatar
Florian Wittkamp committed
2
 * Copyright (C) 2016  For the list of authors, see file AUTHORS.
Tilman Steinweg's avatar
Tilman Steinweg committed
3
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
4
 * This file is part of IFOS.
5
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
6
 * IFOS is free software: you can redistribute it and/or modify
Tilman Steinweg's avatar
Tilman Steinweg committed
7 8
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 2.0 of the License only.
9
 *
Florian Wittkamp's avatar
Florian Wittkamp committed
10
 * IFOS is distributed in the hope that it will be useful,
Tilman Steinweg's avatar
Tilman Steinweg committed
11 12 13
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
14
 *
Tilman Steinweg's avatar
Tilman Steinweg committed
15
 * You should have received a copy of the GNU General Public License
Florian Wittkamp's avatar
Florian Wittkamp committed
16
 * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
17
 -----------------------------------------------------------------------------------------*/
Tilman Steinweg's avatar
Tilman Steinweg committed
18 19 20 21 22 23 24

/*
 * Smoothing gradient / model with a median filter
 */

#include "fd.h"

25
void smooth(float ** mat, int sws, int filter, float Vs_avg, float F_LOW_PASS)
Tilman Steinweg's avatar
Tilman Steinweg committed
26
{
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
    
    /* extern variables */
    
    extern float DH, A;
    extern int FREE_SURF, NX, NY, NXG, NYG;
    extern int NPROCX, NPROCY, MYID, POS[3],WAVETYPE,VERBOSE;
    extern char JACOBIAN[STRING_SIZE], INV_MODELFILE[STRING_SIZE];
    extern FILE *FP;
    extern int FILT_SIZE_GRAD, FILT_SIZE, MODEL_FILTER, GRAD_FILTER, TIME_FILT,GRAD_FILT_WAVELENGTH;
    extern int VERBOSE;
    /* local variables */
    int i, j, ii, jj;
    int i1, j1, filtsize, hfs;
    
    float **model_tmp, **model_med, **filterpart, grad, normgauss, smooth_meter;
    
    char jac_tmp[STRING_SIZE];
    
    FILE *model;
    
    char modfile[STRING_SIZE];
    
49 50 51
    float ** global_matrix;
    
    global_matrix=get_global_from_local_matrix(mat);
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
    switch (filter){
        case 1:
            if((GRAD_FILT_WAVELENGTH==1)&&(TIME_FILT==1)){
                if(VERBOSE) printf("\n -------------------------------------------------------------------------- \n");
                if(VERBOSE) printf("\n Calculating a wavelength dependent filter size for smoothing the gradient: \n");
                FILT_SIZE_GRAD = (int)(Vs_avg/F_LOW_PASS*A/DH);
                if(VERBOSE) printf("\n FILT_SIZE_GRAD = Vs_avg = %4.2f m/s / F_LOW_PASS = %4.2f Hz * weighting factor A = %4.2f / grid spacing DH = %4.2f m  \n",Vs_avg,F_LOW_PASS,A,DH);
                if(VERBOSE) printf("\n New FILT_SIZE_GRAD = %d (grid points) is used (-> %4.2f m).                \n",FILT_SIZE_GRAD,FILT_SIZE_GRAD*DH);
            }
            if (FILT_SIZE_GRAD==0)	return;
            if (!(FILT_SIZE_GRAD % 2)) {
                if (FILT_SIZE_GRAD > 0)	FILT_SIZE_GRAD += 1;
                else			FILT_SIZE_GRAD -= 1;
            }
            hfs = abs(FILT_SIZE_GRAD)/2;
            if(VERBOSE) printf("\n ----------------------------------------------------------------\n");
            if(VERBOSE) printf("\n Filter size is %d gridpoints, half filter size is %d gridpoints.\n",FILT_SIZE_GRAD,hfs);
            filterpart=matrix(1,abs(FILT_SIZE_GRAD),1,abs(FILT_SIZE_GRAD));
            model_tmp = matrix(-hfs+1,NYG+hfs,-hfs+1,NXG+hfs);
            break;
            
        case 2:
            if (FILT_SIZE==0)	return;
            if (!(FILT_SIZE % 2)) {
                if (FILT_SIZE > 0)	FILT_SIZE += 1;
                else			FILT_SIZE -= 1;
            }
            hfs = abs(FILT_SIZE)/2;
            if(VERBOSE) printf("\n ----------------------------------------------------------------\n");
            if(VERBOSE) printf("\n Filter size is %d gridpoints, half filter size is %d gridpoints.\n",FILT_SIZE,hfs);
            filterpart=matrix(1,abs(FILT_SIZE),1,abs(FILT_SIZE));
            model_tmp = matrix(-hfs+1,NYG+hfs,-hfs+1,NXG+hfs);
            break;
86
    }
87
    model_med = matrix(1,NYG,1,NXG);
88
    
89 90 91 92 93
    /* load merged model */
    for (i=1;i<=NXG;i++){
        for (j=1;j<=NYG;j++){
            model_tmp[j][i]=global_matrix[j][i];
        }
94 95 96
    }
    
    
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
    /* apply 2D-Gaussian filter on vp and vs model */
    /* extrapolate array */
    /* left/right boundary */
    for (j=1;j<=NYG;j++){
        for (i=-hfs+1;i<=0;i++){
            model_tmp[j][i] = model_tmp[j][1];}
        for (i=NXG+1;i<=NXG+hfs;i++){
            model_tmp[j][i] = model_tmp[j][NXG];}
    }
    /* top/bottom boundary incl. corners */
    for (j=-hfs+1;j<=0;j++){
        for (i=-hfs+1;i<=NXG+hfs;i++){
            model_tmp[j][i] = model_tmp[1][i];}
    }
    for (j=NYG+1;j<=NYG+hfs;j++){
        for (i=-hfs+1;i<=NXG+hfs;i++){
            model_tmp[j][i] = model_tmp[NYG][i];}
    }
    
    /* filter */
    for (j=1;j<=NYG;j++){
118
        for (i=1;i<=NXG;i++){
119 120 121 122 123
            /* create a filtersize x filtersize matrix */
            for (ii=-hfs;ii<=hfs;ii++){
                for (jj=-hfs;jj<=hfs;jj++){
                    
                    filterpart[ii+hfs+1][jj+hfs+1] = model_tmp[j+jj][i+ii];
124 125
                }
            }
126 127 128 129 130 131 132 133
            /* filter */
            switch (filter){
                case 1:
                    model_med[j][i] = median2d(filterpart,abs(FILT_SIZE_GRAD),abs(FILT_SIZE_GRAD));
                    break;
                case 2:
                    model_med[j][i] = median2d(filterpart,abs(FILT_SIZE),abs(FILT_SIZE));
                    break;
134
            }
135
            
136
        }
137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
    }
    
    get_local_from_global_matrix(model_med, mat);
    
    switch (filter){
        case 1:
            free_matrix(model_tmp,-hfs+1,NYG+hfs,-hfs+1,NXG+hfs);
            free_matrix(filterpart,1,abs(FILT_SIZE_GRAD),1,abs(FILT_SIZE_GRAD));
            break;
        case 2:
            free_matrix(model_tmp,-hfs+1,NYG+hfs,-hfs+1,NXG+hfs);
            free_matrix(filterpart,1,abs(FILT_SIZE),1,abs(FILT_SIZE));
            break;
    }
    free_matrix(model_med,1,NXG,1,NYG);
    
    
154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
    
    MPI_Barrier(MPI_COMM_WORLD);
    
    switch (filter){
        case 1:
            smooth_meter=FILT_SIZE_GRAD*DH;
            if(VERBOSE) fprintf(FP,"\n Gradient %s is smoothed with filter length of %4.2f meter.\n",jac_tmp,smooth_meter);
            break;
        case 2:
            smooth_meter=FILT_SIZE*DH;
            if(VERBOSE) fprintf(FP,"\n Model %s is smoothed with filter length of %4.2f meter.\n",jac_tmp,smooth_meter);
            break ;
    }
    
    
    if(VERBOSE) fprintf(FP,"\n Smoothed gradient / model is distributed on computational nodes.\n");
    if(VERBOSE) fprintf(FP,"\n ----------------------------------------------------------------\n");
Tilman Steinweg's avatar
Tilman Steinweg committed
171
}/* end of smoothing */