annotate Meerwald/dwt.c @ 1:647ab0f05aae v0.4

drop, not needed
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Tue, 14 Aug 2007 19:58:25 +0200
parents be303a3f5ea8
children acb6967ee76d
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
1 #include "wm.h"
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
2 #include "dwt.h"
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
3
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
4 char filter_file[MAXPATHLEN] = "";
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
5 AllFilters dwt_allfilters;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
6 FilterGH *dwt_filters = NULL;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
7 int dwt_method;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
8 int dwt_cols;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
9 int dwt_rows;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
10 int dwt_levels;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
11 int dwt_filter;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
12
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
13 void init_dwt(int cols, int rows, const char *filter_name, int filter, int level, int method) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
14 int i;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
15
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
16 if (strcmp(filter_file, filter_name)) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
17 if (filter_name)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
18 strcpy(filter_file, filter_name);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
19 else
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
20 strcpy(filter_file, "filter.dat");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
21 /* memory leak here - there is no function unload_filters() */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
22 dwt_allfilters = load_filters(filter_file);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
23 if (!dwt_allfilters) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
24 fprintf(stderr, "init_dwt(): unable to open filter definition file %s\n", filter_file);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
25 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
26 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
27 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
28
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
29 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
30 if (level <= 0 || level > rint(log(MIN(cols, rows))/log(2.0)) - 2) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
31 fprintf(stderr, "init_dwt(): level parameter does not match image width/height\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
32 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
33 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
34 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
35
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
36 if (dwt_filters && level != dwt_levels) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
37 free(dwt_filters);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
38 dwt_filters = NULL;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
39 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
40
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
41 dwt_levels = level;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
42
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
43 if (!dwt_filters)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
44 dwt_filters = calloc(level + 1, sizeof(FilterGH));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
45
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
46 for (i = 0; i < level + 1; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
47 dwt_filters[i] = (dwt_allfilters->filter)[filter];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
48
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
49 dwt_filter = filter;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
50 dwt_method = method;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
51 dwt_cols = cols;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
52 dwt_rows = rows;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
53 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
54
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
55 Image_tree fdwt(gray **pixels) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
56 Image image;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
57 Image_tree tree;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
58 int i, j;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
59
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
60 image = new_image(dwt_cols, dwt_rows);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
61
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
62 for (i = 0; i < dwt_rows; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
63 for (j = 0; j < dwt_cols; j++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
64 set_pixel(image, j, i, pixels[i][j]);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
65
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
66 tree = wavelettransform(image, dwt_levels, dwt_filters, dwt_method);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
67 free_image(image);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
68
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
69 return tree;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
70 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
71
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
72 Image_tree fdwt_wp(gray **pixels) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
73 Image image;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
74 Image_tree tree;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
75 int i, j;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
76
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
77 image = new_image(dwt_cols, dwt_rows);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
78
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
79 for (i = 0; i < dwt_rows; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
80 for (j = 0; j < dwt_cols; j++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
81 set_pixel(image, j, i, pixels[i][j]);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
82
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
83 tree = wavelettransform_wp(image, dwt_levels, dwt_filters, dwt_method);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
84 free_image(image);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
85
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
86 return tree;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
87 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
88
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
89 void idwt(Image_tree dwts, gray **pixels) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
90 Image image;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
91 int i, j;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
92
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
93 image = inv_transform(dwts, dwt_filters, dwt_method + 1);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
94
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
95 for (i = 0; i < dwt_rows; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
96 for (j = 0; j < dwt_cols; j++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
97 pixels[i][j] = PIXELRANGE((int) (get_pixel(image, j, i) + 0.5));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
98
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
99 free_image(image);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
100 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
101
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
102 void idwt_wp(Image_tree dwts, gray **pixels) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
103 Image image;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
104 int i, j;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
105
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
106 image = inv_transform(dwts, dwt_filters, dwt_method + 1);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
107
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
108 for (i = 0; i < dwt_rows; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
109 for (j = 0; j < dwt_cols; j++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
110 pixels[i][j] = PIXELRANGE((int) (get_pixel(image, j, i) + 0.5));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
111
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
112 free_image(image);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
113 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
114
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
115 int gen_pollen_filter(double *filter, double alpha, double beta, int which) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
116 int i, j, k, filterlength;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
117 double tf[6];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
118
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
119 /* parameter alpha, beta have to be in range -Pi .. Pi */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
120 if (alpha < -M_PI || alpha >= M_PI) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
121 fprintf(stderr, "alpha %f out of range\n", alpha);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
122 return -1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
123 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
124
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
125 if (beta < -M_PI || beta >= M_PI) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
126 fprintf(stderr, "beta %f out of range\n", beta);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
127 return -1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
128 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
129
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
130 /* generate Pollen filter coefficients, see http://www.dfw.net/~cody for details */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
131 tf[0] = ((1.0 + cos(alpha) + sin(alpha)) * (1.0 - cos(beta) - sin(beta)) + 2.0 * sin(beta) * cos(alpha)) / 4.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
132 tf[1] = ((1.0 - cos(alpha) + sin(alpha)) * (1.0 + cos(beta) - sin(beta)) - 2.0 * sin(beta) * cos(alpha)) / 4.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
133 tf[2] = (1.0 + cos(alpha - beta) + sin(alpha - beta)) / 2.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
134 tf[3] = (1.0 + cos(alpha - beta) - sin(alpha - beta)) / 2.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
135 tf[4] = 1.0 - tf[0] - tf[2];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
136 tf[5] = 1.0 - tf[1] - tf[3];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
137
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
138 /* set close-to-zero filter coefficients to zero */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
139 for (i = 0; i < 6; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
140 if (fabs(tf[i]) < 1.0e-15) tf[i] = 0.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
141
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
142 /* find the first non-zero wavelet coefficient */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
143 i = 0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
144 while (tf[i] == 0.0) i++;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
145
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
146 /* find the last non-zero wavelet coefficient */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
147 j = 5;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
148 while (tf[j] == 0.0) j--;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
149
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
150 filterlength = j - i + 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
151 for (k = 0; k < filterlength; k++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
152 switch (which) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
153 case FILTERH:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
154 filter[k] = tf[j--] / 2.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
155 break;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
156 case FILTERG:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
157 filter[k] = (double) (((i & 0x01) * 2) - 1) * tf[i++] / 2.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
158 break;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
159 case FILTERHi:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
160 filter[k] = tf[j--];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
161 break;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
162 case FILTERGi:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
163 filter[k] = (double) (((i & 0x01) * 2) - 1) * tf[i++];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
164 break;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
165 default:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
166 return -1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
167 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
168
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
169 while (k < 6)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
170 filter[k++] = 0.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
171
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
172 return filterlength;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
173 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
174
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
175 void dwt_pollen_filter(double alpha, double beta) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
176 FilterGH filter;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
177 int i;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
178
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
179 filter = malloc(sizeof(struct FilterGHStruct));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
180 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
181 if (!filter) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
182 fprintf(stderr, "dwt_pollen_filter(): malloc failed()\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
183 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
184 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
185 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
186
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
187 filter->type = FTOther;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
188 filter->name = "pollen";
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
189
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
190 filter->g = new_filter(6);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
191 filter->g->type = FTSymm;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
192 filter->g->hipass = 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
193 filter->g->len = gen_pollen_filter(filter->g->data, alpha, beta, FILTERG);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
194 filter->g->start = -filter->g->len / 2;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
195 filter->g->end = filter->g->len / 2 - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
196
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
197 filter->h = new_filter(6);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
198 filter->h->type = FTSymm;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
199 filter->h->hipass = 0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
200 filter->h->len = gen_pollen_filter(filter->h->data, alpha, beta, FILTERH);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
201 filter->h->start = -filter->h->len / 2;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
202 filter->h->end = filter->h->len / 2 - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
203
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
204 filter->gi = new_filter(6);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
205 filter->gi->type = FTSymm;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
206 filter->gi->hipass = 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
207 filter->gi->len = gen_pollen_filter(filter->gi->data, alpha, beta, FILTERGi);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
208 filter->gi->start = -filter->gi->len / 2;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
209 filter->gi->end = filter->gi->len / 2 - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
210
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
211 filter->hi = new_filter(6);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
212 filter->hi->type = FTSymm;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
213 filter->hi->hipass = 0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
214 filter->hi->len = gen_pollen_filter(filter->hi->data, alpha, beta, FILTERHi);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
215 filter->hi->start = -filter->hi->len / 2;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
216 filter->hi->end = filter->hi->len / 2 - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
217
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
218 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
219 if (dwt_levels <= 0) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
220 fprintf(stderr, "dwt_pollen_filter(): level invalid - set to zero\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
221 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
222 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
223 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
224
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
225 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
226 if (!dwt_filters) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
227 fprintf(stderr, "dwt_pollen_filter(): wm_dwt not initialized, call init_dwt() first\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
228 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
229 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
230 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
231
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
232 for (i = 0; i < dwt_levels + 1; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
233 dwt_filters[i] = filter;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
234 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
235
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
236 int gen_param_filter(double *filter, int n, double alpha[], int which) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
237 int i, j, k, filterlength;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
238 double *tf, *t;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
239
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
240 tf = malloc(2 * (n + 1) * sizeof(double));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
241 t = malloc(2 * (n + 1) * sizeof(double));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
242 if (!tf) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
243 fprintf(stderr, "gen_param_filter(): malloc() failed\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
244 return -1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
245 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
246
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
247 tf[0] = 1.0 / sqrt(2.0);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
248 tf[1] = 1.0 / sqrt(2.0);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
249
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
250 for (k = 0; k < n; k++) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
251 for (i = 0; i < 2 * (k + 2); i++) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
252
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
253 #define H(X) (((X) < 0 || (X) >= 2 * (k + 1)) ? 0.0 : tf[X])
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
254
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
255 t[i] = 0.5 * (H(i - 2) + H(i) +
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
256 cos(alpha[k]) * (H(i - 2) - H(i)) +
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
257 (i & 1 ? -1.0 : 1.0) * sin(alpha[k]) * (H(2 * (k + 2) - i - 1) - H(2 * (k + 2) - i - 3)));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
258 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
259 for (i = 0; i < 2 * (k + 2); i++) tf[i] = t[i];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
260 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
261
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
262 /* set close-to-zero filter coefficients to zero */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
263 for (i = 0; i < 2 * (n + 1) ; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
264 if (fabs(tf[i]) < 1.0e-15) tf[i] = 0.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
265
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
266 /* find the first non-zero wavelet coefficient */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
267 i = 0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
268 while (tf[i] == 0.0) i++;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
269
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
270 /* find the last non-zero wavelet coefficient */
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
271 j = 2 * (n + 1) - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
272 while (tf[j] == 0.0) j--;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
273
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
274 filterlength = j - i + 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
275 for (k = 0; k < filterlength; k++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
276 switch (which) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
277 case FILTERG:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
278 case FILTERGi:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
279 filter[k] = (double) (((i+1 & 0x01) * 2) - 1) * tf[i++];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
280 break;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
281 case FILTERH:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
282 case FILTERHi:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
283 filter[k] = tf[j--];
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
284 break;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
285 default:
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
286 return -1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
287 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
288
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
289 while (k < 2 * (n + 1))
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
290 filter[k++] = 0.0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
291
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
292 return filterlength;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
293 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
294
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
295 void dwt_param_filter(double alpha[], int n) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
296 FilterGH filter;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
297 int i;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
298
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
299 filter = malloc(sizeof(struct FilterGHStruct));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
300 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
301 if (!filter) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
302 fprintf(stderr, "dwt_param_filter(): malloc failed()\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
303 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
304 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
305 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
306
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
307 filter->type = FTOrtho;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
308 filter->name = "param";
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
309
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
310 filter->g = new_filter(2 * (n + 1));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
311 filter->g->type = FTSymm;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
312 filter->g->hipass = 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
313 filter->g->len = gen_param_filter(filter->g->data, n, alpha, FILTERG);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
314 filter->g->start = -filter->g->len / 2;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
315 filter->g->end = filter->g->len / 2 - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
316
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
317 filter->h = new_filter(2 * (n + 1));
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
318 filter->h->type = FTSymm;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
319 filter->h->hipass = 0;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
320 filter->h->len = gen_param_filter(filter->h->data, n, alpha, FILTERH);
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
321 filter->h->start = -filter->h->len / 2;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
322 filter->h->end = filter->h->len / 2 - 1;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
323
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
324 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
325 if (dwt_levels <= 0) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
326 fprintf(stderr, "dwt_pollen_filter(): level invalid - set to zero\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
327 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
328 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
329 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
330
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
331 #ifdef DEBUG
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
332 if (!dwt_filters) {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
333 fprintf(stderr, "dwt_pollen_filter(): wm_dwt not initialized, call init_dwt() first\n");
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
334 return;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
335 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
336 #endif
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
337
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
338 for (i = 0; i < dwt_levels + 1; i++)
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
339 dwt_filters[i] = filter;
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
340 }
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
341
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
342 void done_dwt() {
Peter Meerwald <pmeerw@cosy.sbg.ac.at>
parents:
diff changeset
343 }

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.