1 /**
2  *  @file    smooth.c
3  *  @author  Nathan Baker
4  *  @brief   Convolve grid data with various filters
5  *  @version $Id$
6  */
7 
8 /*
9  *  Last update: 08/29/2016 by Leighton Wilson
10  *  Description: Added ability to read in and output binary DX files
11  */
12 
13 #include "apbs.h"
14 
15 #define IJK(i,j,k)  (((k)*(nx)*(ny))+((j)*(nx))+(i))
16 #define ERRRC 2
17 
18 typedef enum Smooth_Filter {
19     SM_GAUSSIAN  /**< Gaussian filter */
20 } Smooth_Filter;
21 
22 
23 VEMBED(rcsid="$Id$")
24 
25 int gaussian(Vgrid *grid, double stddev, double bandwidth);
26 
usage(int rc)27 int usage(int rc) {
28 
29     char *usage = "\n\n\
30     ----------------------------------------------------------------------\n\
31     This driver program convolves grid data with various filters.  It is\n\
32     invoked as:\n\
33       smooth <args>\n\n\
34     where <args> is a list of the following arguments:\n\
35       REQUIRED GENERAL ARGUMENTS:\n\
36       --format=<format>  where <format> specifies the data format and is one\n\
37                           of the following:\n\
38                                 dx (standard OpenDX)\n\
39                                 dxbin (binary OpenDX)\n\
40       --input=<file>     where <file> is the input mesh data in the\n\
41                          specified format\n\
42       --output=<file>    where <file> is the output mesh data in the\n\
43                          specified format\n\
44       --filter=<filter>  where <filter> is the filter with which the data\n\
45                          will be convolved and is one of the following:\n\
46                          gaussian (Gaussian filter)\n\
47       REQUIRED FILTER-SPECIFIC ARGUMENTS:\n\
48         Gaussian filter:\n\
49         --stddev=<n>     the standard deviation of the filter (in A)\n\
50         --bandwidth=<n>  the bandwith of the filter (in units of stddev)\n\
51     ----------------------------------------------------------------------\n\n";
52 
53     Vnm_print(2, usage);
54 
55     exit(rc);
56 
57     return 0;
58 }
59 
main(int argc,char ** argv)60 int main(int argc, char **argv) {
61 
62     /* *************** VARIABLES ******************* */
63     int i;
64     Vgrid *grid = VNULL;
65     /* Input parameters */
66     Vdata_Format format; int gotFormat = 0;
67     char inPath[VMAX_BUFSIZE]; int gotInPath = 0;
68     char outPath[VMAX_BUFSIZE]; int gotOutPath = 0;
69     Smooth_Filter filter; int gotFilter = 0;
70     double stddev; int gotStddev = 0;
71     double bandwidth; int gotBandwidth = 0;
72 
73     char *header = "\n\n\
74     ----------------------------------------------------------------------\n\
75     ----------------------------------------------------------------------\n\
76     \n\n";
77 
78     /* *************** CHECK INVOCATION ******************* */
79     Vio_start();
80     Vnm_redirect(1);
81     Vnm_print(1, "%s", header);
82     for (i=1; i<argc; i++) {
83         Vnm_print(1, "Parsing: %s...\n", argv[i]);
84         if (strstr(argv[i], "--format") != NULL) {
85             if (strstr(argv[i], "dxbin") != NULL) {
86                 gotFormat = 1;
87                 format = VDF_DXBIN;
88             } else if (strstr(argv[i], "dx") != NULL) {
89                 gotFormat = 1;
90                 format = VDF_DX;
91             } else {
92                 Vnm_print(2, "Error:  %s\n", argv[i]);
93                 usage(2);
94             }
95         } else if (strstr(argv[i], "--input") != NULL) {
96             if (sscanf(argv[i], "--input=%s", inPath) == 1) gotInPath = 1;
97             else {
98                 Vnm_print(2, "Error:  %s\n", argv[i]);
99                 usage(2);
100             }
101         } else if (strstr(argv[i], "--output") != NULL) {
102             if (sscanf(argv[i], "--output=%s", outPath) == 1) gotOutPath = 1;
103             else {
104                 Vnm_print(2, "Error:  %s\n", argv[i]);
105                 usage(2);
106             }
107         } else if (strstr(argv[i], "--filter") != NULL) {
108             if (strstr(argv[i], "gaussian") != NULL) {
109                 gotFilter = 1;
110                 filter = SM_GAUSSIAN;
111             } else {
112                 Vnm_print(2, "Error:  %s\n", argv[i]);
113                 usage(2);
114             }
115         } else if (strstr(argv[i], "--stddev") != NULL) {
116             if (sscanf(argv[i], "--stddev=%lf", &stddev) == 1) gotStddev = 1;
117             else {
118                 Vnm_print(2, "Error:  %s\n", argv[i]);
119                 usage(2);
120             }
121         } else if (strstr(argv[i], "--bandwidth") != NULL) {
122             if (sscanf(argv[i], "--bandwidth=%lf", &bandwidth) == 1)
123                 gotBandwidth = 1;
124             else {
125                 Vnm_print(2, "Error:  %s\n", argv[i]);
126                 usage(2);
127             }
128         } else {
129             Vnm_print(2, "Error:  %s\n", argv[i]);
130             usage(2);
131         }
132     }
133     if (!gotFormat) {
134         Vnm_print(2, "Error:  --format not specified!\n");
135         usage(2);
136     }
137     if (!gotInPath) {
138         Vnm_print(2, "Error:  --input not specified!\n");
139         usage(2);
140     }
141     if (!gotOutPath) {
142         Vnm_print(2, "Error:  --output not specified!\n");
143         usage(2);
144     }
145     if (!gotFilter) {
146         Vnm_print(2, "Error:  --filter not specified!\n");
147         usage(2);
148     }
149     if (filter == SM_GAUSSIAN) {
150         if (!gotStddev) {
151             Vnm_print(2, "Error:  --stddev not specified!\n");
152             usage(2);
153         }
154         if (!gotBandwidth) {
155             Vnm_print(2, "Error:  --bandwidth not specified!\n");
156             usage(2);
157         }
158     }
159 
160     /* *************** READ DATA ******************* */
161     Vnm_print(1, "main:  Reading data from %s...\n", inPath);
162     grid = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
163     if (format == VDF_DX) {
164         if (!Vgrid_readDX(grid, "FILE", "ASC", VNULL, inPath)) {
165             Vnm_print(2, "main:  Problem reading standard OpenDX-format grid from %s\n",
166               inPath);
167             return ERRRC;
168         }
169     } else if (format == VDF_DXBIN) {
170         if (!Vgrid_readDXBIN(grid, "FILE", "ASC", VNULL, inPath)) {
171             Vnm_print(2, "main:  Problem reading binary OpenDX-format grid from %s\n",
172               inPath);
173             return ERRRC;
174         }
175     }
176 
177     /* *************** SMOOTH ******************* */
178     switch(filter) {
179         case SM_GAUSSIAN:
180            Vnm_print(1, "Smoothing data with Gaussian filter...\n");
181            gaussian(grid, stddev, bandwidth);
182            break;
183         default:
184            Vnm_print(2, "Invalid format (%d)!\n", format);
185            usage(2);
186     }
187 
188     /* *************** READ DATA ******************* */
189     Vnm_print(1, "main:  Writing data to %s...\n", outPath);
190     if (format == VDF_DX)
191       Vgrid_writeDX(grid, "FILE", "ASC", VNULL, outPath, "Smoothed data",
192         VNULL);
193     else if (format == VDF_DXBIN)
194       Vgrid_writeDXBIN(grid, "FILE", "ASC", VNULL, outPath, "Smoothed data",
195         VNULL);
196 
197     return 0;
198 }
199 
gaussian(Vgrid * grid,double stddev,double bandwidth)200 int gaussian(Vgrid *grid, double stddev, double bandwidth) {
201 
202     int nx, ny, nz, iband, jband, kband, i, j, k, ii, jj, kk;
203     int kkmin, jjmin, iimin, kkmax, jjmax, iimax;
204     double hx, hy, hzed, xmin, ymin, zmin, *newData, *oldData, scal, norm;
205     double u, ga, dist2;
206 
207     Vnm_print(1, "Gaussian filter:  std. dev. = %g A, bandwidth = %g A.\n",
208       stddev, bandwidth*stddev);
209 
210     nx = grid->nx; ny = grid->ny; nz = grid->nz;
211     hx = grid->hx; hy = grid->hy; hzed = grid->hzed;
212     xmin = grid->xmin; ymin = grid->ymin; zmin = grid->zmin;
213     Vnm_print(1, "Grid:  %d x %d x %d points\n", nx, ny, nz);
214     Vnm_print(1, "Grid:  %g, %g, %g A spacing\n", hx, hy, hzed);
215     Vnm_print(1, "Grid:  (%g, %g, %g) A origin\n", xmin, ymin, zmin);
216 
217     /* Convert HALF bandwidth to grid units */
218     iband = (int)(stddev*bandwidth/hx);
219     jband = (int)(stddev*bandwidth/hy);
220     kband = (int)(stddev*bandwidth/hzed);
221     /* Special handling for iband, jband and kband, they are non-zero positive integers */
222     if (iband == 0) iband = 1;
223     if (jband == 0) jband = 1;
224     if (kband == 0) kband = 1;
225     Vnm_print(1, "Bandwidth converted to %d x %d x %d grid units.\n",
226       iband, jband, kband);
227     Vnm_print(1, "This means any non-zero data within (%g, %g, %g) of the\n",
228       (iband+1)*hx, (jband+1)*hy, (kband+1)*hzed);
229     Vnm_print(1, "domain boundary will be convolved differently.\n");
230 
231     /* Get exponent scaling factor */
232     scal = 2.0 * stddev * stddev;
233     VASSERT(scal > 0);
234     scal = 1.0/scal;
235 
236     /* Get data */
237     oldData = grid->data;
238     newData = Vmem_malloc(VNULL, (nx*ny*nz), sizeof(double));
239 
240     /* Copy over old data.  All but the boundary values will be replaced in the next step so this is more copying than is strictly
241        necessary... */
242     for (i=0; i<(nx*ny*nz); i++) newData[i] = oldData[i];
243 
244     /* Apply filter */
245     for (k=1; k<(nz-1); k++) {
246         kkmin = VMAX2(1, (k - kband));
247         kkmax = VMIN2((nz-1), (k + kband));
248         for (j=1; j<(ny-1); j++) {
249             jjmin = VMAX2(1, (j - jband));
250             jjmax = VMIN2((ny-1), (j + jband));
251             for (i=1; i<(nx-1); i++) {
252                 iimin = VMAX2(1, (i - iband));
253                 iimax = VMIN2((nx-1), (i + iband));
254                 u = 0;
255                 /* We normalize within the loop to conserve densities */
256                 norm = 0;
257                 for (kk=kkmin; kk<kkmax; kk++) {
258                     for (jj=jjmin; jj<jjmax; jj++) {
259                         for (ii=iimin; ii<iimax; ii++) {
260                             dist2 = VSQR(hx*(i-ii)) + VSQR(hy*(j-jj)) +
261                               VSQR(hzed*(k-kk));
262                             ga = VEXP(-dist2*scal);
263                             u += (ga*oldData[IJK(ii,jj,kk)]);
264                             norm += ga;
265                         } /* ii loop */
266                     } /* jj loop */
267                 } /* kk loop */
268                 if (norm > 0) newData[IJK(i,j,k)] = u/norm;
269                 else newData[IJK(i,j,k)] = 0;
270             } /* i loop */
271         } /* j loop */
272     } /* k loop */
273 
274     /* Replace data */
275     for (i=0; i<(nx*ny*nz); i++) grid->data[i] = newData[i];
276     Vmem_free(VNULL, (nx*ny*nz), sizeof(double), (void **)&newData);
277 
278     return 0;
279 }
280