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