1 /* @(#)transform.c 19.1 (ES0-DMD) 02/25/03 13:34:40 */
2 /*===========================================================================
3 Copyright (C) 1995 European Southern Observatory (ESO)
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 2 of
8 the License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public
16 License along with this program;
17 If not, see <http://www.gnu.org/licenses/>.
18
19 Corresponding concerning ESO-MIDAS should be addressed as follows:
20 Internet e-mail: midas@eso.org
21 Postal address: European Southern Observatory
22 Data Management Division
23 Karl-Schwarzschild-Strasse 2
24 D 85748 Garching bei Muenchen
25 GERMANY
26 ===========================================================================*/
27
28 /******************************************************************************
29 ** Copyright (C) 1993 by European Southern Observatory
30 *******************************************************************************
31 **
32 ** UNIT
33 **
34 ** Version: 19.1
35 **
36 ** Author: Jean-Luc Starck
37 **
38 ** Date: 03/02/25
39 **
40 ** File: transform.c
41 **
42 *******************************************************************************
43 **
44 ** DESCRIPTION This module contains wavelet transform routines
45 ** -----------
46 ******************************************************************************
47 **
48 ** wavelet_transform_file (File_Name_Imag, File_Name_Transform,
49 ** Type_Transform, Fc, Nbr_Plan)
50 ** char *File_Name_Imag, *File_Name_Transform;
51 ** int Type_Transform;
52 ** float Fc;
53 ** int Nbr_Plan;
54 **
55 ** Computes the wavelet transform of the image with a file name File_Name_Imag
56 ** and write the result in the file File_Name_Transform
57 **
58 ** File_Name_Imag = File name of the imput image
59 ** File_Name_Transform = File name of the output wavelet transform
60 ** Type_Transform = wavelet transform algorithm number
61 ** Fc = cut-off frequency if the algorithm use the FFT
62 ** Nbr_Plan = number of scales
63 **
64 ******************************************************************************
65 **
66 ** wavelet_transform_data (Imag, Nl, Nc, Wavelet, Type_Transform, Fc, Nbr_Plan)
67 ** float *Imag;
68 ** int Nl, Nc;
69 ** wave_transf_des *Wavelet;
70 ** int Type_Transform;
71 ** float Fc;
72 ** int Nbr_Plan;
73 **
74 ** Computes the wavelet transform of the image Imag
75 ** and write the result in the structure Wavelet
76 **
77 ** Imag = INPUT:image
78 ** Wavelet = OUTPUT:wavelet
79 ** Nl,Nc = INPUT: number of lines and columns
80 ** Type_Transform = INPUT:wavelet transform algorithm number
81 ** Which_Algorithm = 1...7
82 ** 1 ==> a trous algorithm with a linear scaling function
83 ** the wavelet function is the difference betwwen two resolutions
84 ** 2 ==> a trous algorithm with a B3-spline scaling function
85 ** the wavelet function is the difference betwwen two resolutions
86 ** 3 ==> algorithm using the Fourier transform:
87 ** without any reduction of the samples between two scales
88 ** the Fourier transform of the scaling function is a b3-spline
89 ** the wavelet function is the difference between two resolutions
90 ** 4 ==> pyramidal algorithm in the direct space, with a linear
91 ** scaling function
92 ** 5 ==> pyramidal algorithm in the direct space, with a b3-spline
93 ** scaling function
94 ** 6 ==> algorithm using the Fourier transform:
95 ** with a reduction of the samples between two scales
96 ** the Fourier transform of the scaling function is a b3-spline
97 ** the wavelet function is the difference between two resolutions
98 ** 7 ==> algorithm using the Fourier transform:
99 ** with a reduction of the samples between two scales
100 ** the Fourier transform of the scaling function is a b3-spline
101 ** the wavelet function is the difference between the square of
102 ** two resolutions
103 ** 8 ==> Mallat's Algorithm with biorthogonal filters.
104 ** Fc = INPUT:cut-off frequency if the algorithm use the FFT
105 ** Nbr_Plan = INPUT:number of scales
106 **
107 ******************************************************************************/
108
109 #include <stdio.h>
110 #include <math.h>
111 #include <string.h>
112
113 #include "core/siril.h"
114 #include "core/siril_log.h"
115 #include "gui/utils.h"
116 #include "gui/progress_and_log.h"
117 #include "algos/Def_Math.h"
118 #include "algos/Def_Mem.h"
119 #include "algos/Def_Wavelet.h"
120
prepare_rawdata(float * Imag,int Nl,int Nc,WORD * buf)121 int prepare_rawdata(float *Imag, int Nl, int Nc, WORD *buf) {
122 float *im = Imag;
123 int i;
124
125 #ifdef _OPENMP
126 #pragma omp parallel for num_threads(com.max_thread) private(i) schedule(static)
127 #endif
128 for (i = 0; i < (Nl) * (Nc); ++i) {
129 im[i] = (float) buf[i];
130 }
131 return 0;
132 }
133
134 /****************************************************************************/
135
f_vector_alloc(Nbr_Elem)136 float *f_vector_alloc(Nbr_Elem)
137 /* allocates a vector of float */
138 int Nbr_Elem; {
139 float *Vector;
140
141 Vector = (float*) calloc(Nbr_Elem, sizeof(float));
142 if (Vector == NULL) {
143 PRINT_ALLOC_ERR;
144 }
145 return (Vector);
146 }
147
148 /*****************************************************************************/
149
wavelet_transform_file(float * Imag,int Nl,int Nc,char * File_Name_Transform,int Type_Transform,int Nbr_Plan,WORD * data)150 int wavelet_transform_file(float *Imag, int Nl, int Nc,
151 char *File_Name_Transform, int Type_Transform, int Nbr_Plan, WORD *data) {
152 wave_transf_des Wavelet;
153 memset(&Wavelet, 0, sizeof(wave_transf_des));
154
155 /* read the input image */
156 prepare_rawdata(Imag, Nl, Nc, data);
157 snprintf(Wavelet.Name_Imag, MAX_SIZE_NAME_IMAG - 1, "%s",
158 File_Name_Transform);
159 Wavelet.Name_Imag[MAX_SIZE_NAME_IMAG - 1] = '\0';
160
161 if (wavelet_transform_data(Imag, Nl, Nc, &Wavelet, Type_Transform,
162 Nbr_Plan)) {
163 return 1;
164 }
165 if (wave_io_write(File_Name_Transform, &Wavelet)) {
166 return 1;
167 }
168
169 wave_io_free(&Wavelet);
170 return 0;
171 }
172
wavelet_transform_file_float(float * Imag,int Nl,int Nc,char * File_Name_Transform,int Type_Transform,int Nbr_Plan)173 int wavelet_transform_file_float(float *Imag, int Nl, int Nc,
174 char *File_Name_Transform, int Type_Transform, int Nbr_Plan) {
175 wave_transf_des Wavelet;
176 memset(&Wavelet, 0, sizeof(wave_transf_des));
177
178 /* read the input image */
179 snprintf(Wavelet.Name_Imag, MAX_SIZE_NAME_IMAG - 1, "%s",
180 File_Name_Transform);
181 Wavelet.Name_Imag[MAX_SIZE_NAME_IMAG - 1] = '\0';
182
183 if (wavelet_transform_data(Imag, Nl, Nc, &Wavelet, Type_Transform,
184 Nbr_Plan)) {
185 return 1;
186 }
187 if (wave_io_write(File_Name_Transform, &Wavelet)) {
188 return 1;
189 }
190
191 wave_io_free(&Wavelet);
192 return 0;
193 }
194
195 /*****************************************************************************/
196
197 // alternative to wavelet_transform_file that does not use a file
198 // call wave_io_free() on the returned Wavelet if retval is 0
wavelet_transform(float * Imag,int Nl,int Nc,wave_transf_des * Wavelet,int Type_Transform,int Nbr_Plan,WORD * data)199 int wavelet_transform(float *Imag, int Nl, int Nc,
200 wave_transf_des *Wavelet, int Type_Transform, int Nbr_Plan, WORD *data) {
201
202 memset(Wavelet, 0, sizeof(wave_transf_des));
203
204 /* read the input image */
205 prepare_rawdata(Imag, Nl, Nc, data);
206
207 if (wavelet_transform_data(Imag, Nl, Nc, Wavelet, Type_Transform,
208 Nbr_Plan)) {
209 wave_io_free(Wavelet);
210 return 1;
211 }
212
213 return 0;
214 }
215
wavelet_transform_float(float * Imag,int Nl,int Nc,wave_transf_des * Wavelet,int Type_Transform,int Nbr_Plan)216 int wavelet_transform_float(float *Imag, int Nl, int Nc,
217 wave_transf_des *Wavelet, int Type_Transform, int Nbr_Plan) {
218
219 memset(Wavelet, 0, sizeof(wave_transf_des));
220
221 if (wavelet_transform_data(Imag, Nl, Nc, Wavelet, Type_Transform,
222 Nbr_Plan)) {
223 wave_io_free(Wavelet);
224 return 1;
225 }
226 return 0;
227 }
228 /*****************************************************************************/
229
230
wavelet_transform_data(float * Imag,int Nl,int Nc,wave_transf_des * Wavelet,int Type_Transform,int Nbr_Plan)231 int wavelet_transform_data(float *Imag, int Nl, int Nc,
232 wave_transf_des *Wavelet, int Type_Transform, int Nbr_Plan) {
233 float *Pave;
234 double Exp;
235 int Size, Min, temp;
236
237 Wavelet->Nbr_Ligne = Nl;
238 Wavelet->Nbr_Col = Nc;
239 Wavelet->Nbr_Plan = Nbr_Plan;
240 Wavelet->Type_Wave_Transform = Type_Transform;
241
242 /* test if the number of planes is not too high */
243 Min = (Nl < Nc) ? Nl : Nc;
244 Exp = (double) Nbr_Plan + 2.;
245 temp = pow(2., Exp) + 0.5;
246 if (Min < temp) {
247 siril_log_message(_("wavelet_transform_data: bad plane number\n"));
248 return 1;
249 }
250
251 switch (Type_Transform) {
252 case TO_PAVE_LINEAR:
253 case TO_PAVE_BSPLINE:
254 Size = Nl * Nc * Nbr_Plan;
255 Wavelet->Pave.Data = f_vector_alloc(Size);
256 if (Wavelet->Pave.Data == NULL) {
257 PRINT_ALLOC_ERR;
258 return 1;
259 }
260 Pave = Wavelet->Pave.Data;
261 pave_2d_tfo(Imag, Pave, Nl, Nc, Nbr_Plan, Type_Transform);
262 break;
263 default:
264 printf("wavelet_transform_data: wrong transform type\n");
265 return 1;
266 }
267 return 0;
268 }
269
270 /*****************************************************************************/
271
272