← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to contrib/wavelet/src/wa_pers.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/wavelet/src/wa_pers.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence 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
**
 
30
**    UNIT
 
31
**
 
32
**    Version: 19.1
 
33
**
 
34
**    Author: Jean-Luc Starck 
 
35
**
 
36
**    Date:  03/02/25
 
37
**    
 
38
**    File: wa_pers.c
 
39
**
 
40
.VERSION
 
41
 090810         last modif
 
42
 
 
43
*******************************************************************************
 
44
**
 
45
**    DESCRIPTION creation of an image from a wavelet transform
 
46
**    ----------- 
 
47
**                 
 
48
**    PARAMETRES    
 
49
**    ----------   
 
50
** 
 
51
**      File_Name_Transform (keyword: IN_A):
 
52
**             File name of the input wavelet transform
 
53
**
 
54
**      File_Name_Imag (keyword: OUT_A):
 
55
**             File name of the output image
 
56
**
 
57
**      Increment (keyword: INPUTI(1)): (default is 1)
 
58
**             we keep one line on Increment
 
59
**
 
60
**      Visualisation Mode (keyword: IN_B): CO or BW (default is CO)
 
61
**               CO =  Color 
 
62
**               BW = black and white
 
63
**
 
64
**      Threshold value (keyword: INPUTR(1)): (default is 5)
 
65
**               all the wavelet coefficients > Treshold*Sigma_Scale
 
66
**               are set to Treshold*Sigma_Scale
 
67
**
 
68
**    RESULTS : create an image from a wavelet file   
 
69
**    -------   the image is a representation of the 
 
70
**              wavelet coefficicients with a 3 dimension 
 
71
**              visualisation 
 
72
******************************************************************************/
 
73
 
 
74
 
 
75
#include <stdlib.h>
 
76
#include <stdio.h>
 
77
#include <math.h>
 
78
#include <string.h>
 
79
 
 
80
#include <midas_def.h>
 
81
 
 
82
#include "Def_Math.h"
 
83
#include "Def_Mem.h"
 
84
#include "Def_Wavelet.h"
 
85
 
 
86
extern float lib_mat_ecart_type();
 
87
 
 
88
extern void wave_io_read(), wavelet_extract_plan(), io_write_pict_f_to_file();
 
89
 
 
90
 
 
91
#define MAXI 255
 
92
 
 
93
#define COULEUR_1 1
 
94
#define COULEUR_2 1
 
95
#define BLACK_WHITE 0
 
96
 
 
97
int MAXX = 512;
 
98
int MAXY = 512;
 
99
 
 
100
/*****************************************************************************/
 
101
 
 
102
void ligne (x1,y1,x2,y2,Couleur,Nl,Nc,Tab_Caval)
 
103
/* 
 
104
   Draw a line in the image Tab_Caval between the two points:
 
105
    I1,J1 and I2,J2. The intensity of the line is defined by Couleur
 
106
*/
 
107
   
 
108
int x1,y1,x2,y2,Nl,Nc;
 
109
float Couleur;
 
110
float *Tab_Caval;
 
111
{
 
112
    int j,pasx,x,y,dx,dy,cumul;
 
113
    float pasy;
 
114
    
 
115
    if (x1 < x2) pasx = 1;
 
116
    else pasx = -1;
 
117
    if (y1 < y2) pasy = 1;
 
118
    else pasy = -1;
 
119
    dx = abs (x2 - x1);
 
120
    dy = abs (y2 - y1);
 
121
    x= x1;
 
122
    y = y1;
 
123
    Tab_Caval [x1*Nc+y1] = Couleur;
 
124
    
 
125
    if (dx > dy)
 
126
    {
 
127
        cumul = dx / 2;
 
128
        for (j = 0; j < dx; j++)
 
129
        {
 
130
            x += pasx;
 
131
            cumul += dy;
 
132
            if (cumul > dx)
 
133
            {
 
134
                cumul -= dx;
 
135
                y += pasy;
 
136
            }
 
137
            Tab_Caval [x*Nc+y] = Couleur;
 
138
        }
 
139
    }
 
140
    else
 
141
    {
 
142
        
 
143
        cumul = dy / 2;
 
144
        for (j = 0; j < dy; j++)
 
145
        {
 
146
            y += pasy;
 
147
            cumul += dx;
 
148
            if (cumul > dy)
 
149
            {
 
150
                cumul -= dy;
 
151
                x += pasx;
 
152
            }
 
153
            Tab_Caval [x*Nc+y] = Couleur;
 
154
        } 
 
155
    }
 
156
}  
 
157
 
 
158
/*****************************************************************************/
 
159
 
 
160
void lib_caval_pict (Nl, Nc, Maxy, Maxx, Tab_Img, Tab_Caval, Inc, Mod_Visu)
 
161
int Nl,Nc,Maxy,Maxx,Inc;
 
162
float *Tab_Img,*Tab_Caval; 
 
163
int Mod_Visu;
 
164
{
 
165
    int j,Size,Indi,Indj;
 
166
    float i,Seuil,X_Pas,Increment,Nb_Line;
 
167
    float *Tab_Max;
 
168
    int Borne;
 
169
    int X1,Y1,X2,Y2;
 
170
    int Couleur;
 
171
    float Val,H_Facteur,L_Facteur,u,Interp;
 
172
    float Level;
 
173
 
 
174
    Size = Maxx * Maxy;    
 
175
    Tab_Max = f_vector_alloc (Maxx);
 
176
    for (j = 0; j < Maxx; j++) Tab_Max [j] = 0.;
 
177
    for (j = 0; j < Size; j++) Tab_Caval [j] = 0.;
 
178
 
 
179
    Increment  = (float) Inc;
 
180
    {
 
181
        float Min, Max;
 
182
 
 
183
        Min = Max = Tab_Img[0];
 
184
        for (j = 0; j < Nl*Nc; j++)
 
185
        {
 
186
            if (Min > Tab_Img[j]) Min = Tab_Img[j];
 
187
            if (Max < Tab_Img[j]) Max = Tab_Img[j];
 
188
        }
 
189
        for (j = 0; j < Nl*Nc; j++) Tab_Img[j] = (Tab_Img[j]-Min) / (Max - Min);
 
190
    }
 
191
 
 
192
    /* Color threshold  */
 
193
    Seuil = 1. / 5.;
 
194
 
 
195
    /* scale parameter realted to the visualisation angle*/
 
196
    H_Facteur = (float) Maxy / 2.;
 
197
    L_Facteur = (float) Maxx / (float) Nc;
 
198
    X_Pas = 1. / L_Facteur;
 
199
    Nb_Line = (float) Nl / ((float) Maxy / 2.);
 
200
    Borne = (Maxy /2. > Nl)? Nl : Maxy / 2.;
 
201
    
 
202
    for (i = 0; i < Borne; i += Increment)
 
203
    {
 
204
        /* First point position */
 
205
        u = i * Nb_Line;
 
206
        Indi = (int) u * Nc;
 
207
 
 
208
        if (Indi > Nl*Nc)
 
209
        {
 
210
             printf ("(%f,%d)",Nb_Line,Indi);
 
211
        }
 
212
        Val = H_Facteur * Tab_Img [Indi] +  i;
 
213
        if (Val > Tab_Max [0]) Tab_Max [0] = Val;
 
214
        X1 = 0;
 
215
        Y1 = Val;
 
216
        u = 0.;
 
217
        
 
218
        for (j = 1; j < Maxx; j ++)
 
219
        {
 
220
            /* ordinate position */
 
221
            u += X_Pas;
 
222
            u = ((float) Nc / (float) Maxx) * (float) j;
 
223
            Indj = u;
 
224
            
 
225
            if ((L_Facteur > 1.) && ((Indj + 1) < Nc))
 
226
            {
 
227
                Interp = u - Indj;
 
228
                Val  = Tab_Img [Indi+Indj];
 
229
                Val = H_Facteur * (Val + (Tab_Img [Indi+Indj+1] 
 
230
                                             - Val) * Interp) +  i;
 
231
                if ((Val >= Maxy) || (Val < 0.))
 
232
                {
 
233
                    printf ("zoom val (%d,%d)=%f ",Indi,Indj,Val);
 
234
                    exit (0);
 
235
                }
 
236
            }
 
237
            else
 
238
            {
 
239
                Val = H_Facteur * Tab_Img [Indi+Indj] + i;
 
240
                if ((Val >= Maxy) || (Val < 0.))
 
241
                {
 
242
                    printf ("calcul val %5.1f:(%d,%d)=%f ",i,Indi,Indj,Val);
 
243
                    exit (0);
 
244
                }
 
245
            }
 
246
                
 
247
            /* Tab_Max = Z-array: contains the max for each columns   */
 
248
            if (Val > Tab_Max [j])
 
249
            {
 
250
                Tab_Max [j] = Val;
 
251
                X2 = j;
 
252
               /* Y2 = Maxy - 1 - (Tab_Max [j] + 0.5);*/
 
253
                Y2 = Val;
 
254
                if (Y2 < 0)
 
255
                {
 
256
                    printf ("Y2 < 0 : %d",Y2);
 
257
                    exit (0);
 
258
                    Y2 = 0;
 
259
                }
 
260
                if (Tab_Img [Indi+Indj] < Seuil) Couleur = COULEUR_1;
 
261
                else Couleur = COULEUR_2;
 
262
                
 
263
                if ((Couleur == COULEUR_1) 
 
264
                            && ((Tab_Caval [Y1*Maxx+X1] == COULEUR_2)
 
265
                                || (Tab_Caval [Y2*Maxx+X2] == COULEUR_2)))
 
266
                {
 
267
                    X1 = X1+1;
 
268
                    Y1 = Y2;
 
269
                }
 
270
                else
 
271
                {
 
272
                    Level = Tab_Img [Indi+Indj];
 
273
                    if (Y1 < 0)
 
274
                    {
 
275
                        printf ("Y1(%d,%5.1f)",j,Tab_Max [j]);
 
276
                        exit (-1);
 
277
                        Y1 = 0;
 
278
                    }
 
279
                    if (Mod_Visu == BLACK_WHITE)
 
280
                    {
 
281
                        Level = Couleur;
 
282
                    }
 
283
                    ligne (Y1,X1,Y2,X2,Level,Maxy,Maxx,Tab_Caval);
 
284
                    X1 = X2;
 
285
                    Y1 = Y2;
 
286
                }
 
287
            }
 
288
            else
 
289
            {
 
290
                X1 = j;
 
291
               /*  Y1 = Maxy - Tab_Max [j];*/
 
292
                Y1 = Tab_Max [j];
 
293
                if (Y1 < 0)
 
294
                {
 
295
                    printf ("Y1(%d,%f)",j,Tab_Max [j]);
 
296
                    exit (-1);
 
297
                    Y1 = 0;
 
298
                }
 
299
            }
 
300
        }
 
301
    }
 
302
    
 
303
    free ((char *)Tab_Max);
 
304
}
 
305
 
 
306
/*****************************************************************************/
 
307
 
 
308
static void copy_plan_to_pict(Tab_Caval, Nl1, Nc1, Pict_Caval, Maxy, Maxx, Num_Etap)
 
309
float *Tab_Caval,*Pict_Caval;
 
310
int Nl1, Nc1, Maxy, Maxx;
 
311
int Num_Etap;
 
312
{
 
313
    int Offsx,Offsy;
 
314
    int i,j;
 
315
    int ind;
 
316
    
 
317
    Offsx = (Maxx - Nc1) / 2;
 
318
    Offsy = Nl1 * Num_Etap * Maxx;
 
319
    
 
320
    for (i = 0; i < Nl1; i++)
 
321
    {
 
322
        for (j = 0; j < Nc1; j++)
 
323
        {
 
324
            ind = Offsy + Offsx + i * Maxx  + j;
 
325
            Pict_Caval [ind] = Tab_Caval [i * Nc1 + j] * MAXI;
 
326
        }
 
327
    }
 
328
}
 
329
 
 
330
/**************************************************************/
 
331
 
 
332
int main()
 
333
{
 
334
    int i,j;
 
335
    int Nl,Nc;
 
336
    wave_transf_des Wavelet;
 
337
    float *ImaSynt, *Imag, Sigma;
 
338
    string File_Name_Imag, File_Name_Transform;
 
339
    int Actvals, Maxvals, Felem, Unit, Null;
 
340
    int Stat, Buffer_Int;
 
341
    int Nl1,Nc1,Size;
 
342
    int Inc = 1;
 
343
    int Mod_Visu = 1;
 
344
    float Zoom, *Pict_Caval, K_Sigma;
 
345
    string Buffer_Visu;
 
346
    int Size_Plan;
 
347
 
 
348
    /* Initialisation */
 
349
    SCSPRO("visu");
 
350
 
 
351
    /* read the wavelet transform file name */
 
352
    Felem = 1;
 
353
    Maxvals = 60;
 
354
    Stat = SCKGETC("IN_A", Felem, Maxvals, &Actvals, File_Name_Transform);
 
355
 
 
356
    /* read the image output file name */
 
357
    Stat = SCKGETC("OUT_A", Felem, Maxvals, &Actvals, File_Name_Imag);
 
358
 
 
359
    /* read the increment */
 
360
    Maxvals = 1;
 
361
    Stat = SCKRDI("INPUTI",Felem, Maxvals, &Actvals, &Buffer_Int, &Unit, &Null);
 
362
    Inc = Buffer_Int;
 
363
 
 
364
    /* read the visualisation mode (Black and White (BW) or Color (CO))*/
 
365
    Maxvals = 60;
 
366
    Stat = SCKGETC("IN_B", Felem, Maxvals, &Actvals, Buffer_Visu);
 
367
    if (   ((Buffer_Visu[0] == 'B') || (Buffer_Visu[0] == 'b'))
 
368
        && ((Buffer_Visu[1] == 'W') || (Buffer_Visu[1] == 'w'))) Mod_Visu = 0;
 
369
 
 
370
    /* read the threshold value */
 
371
    Maxvals = 1;
 
372
    Stat = SCKRDR("INPUTR", Felem, Maxvals, &Actvals, &K_Sigma, &Unit, &Null);
 
373
 
 
374
   /* read the wavelet */
 
375
    wave_io_read (File_Name_Transform, &Wavelet);
 
376
 
 
377
    Nl = Wavelet.Nbr_Ligne;
 
378
    Nc = Wavelet.Nbr_Col;
 
379
    Zoom = (float) MAXX / (float) Nc;
 
380
    Size_Plan = (float) MAXY / (float) Wavelet.Nbr_Plan;
 
381
    Size = MAXX * MAXY;    
 
382
    Pict_Caval = f_vector_alloc (Size);
 
383
    ImaSynt = f_vector_alloc (Size);
 
384
 
 
385
    for (i = 0; i < Size; i++) ImaSynt[i] = 0.;
 
386
 
 
387
    for (i = 0; i <Wavelet.Nbr_Plan ; i++)
 
388
    {
 
389
        wavelet_extract_plan (&Wavelet, &Imag, &Nl, &Nc, i+1);
 
390
        Sigma = K_Sigma * lib_mat_ecart_type(Imag,Nl,Nc);
 
391
 
 
392
        Nl1 = Size_Plan;
 
393
        Nc1 = (((float) Nc * Zoom) > MAXX) ? MAXX : ((float)Nc * Zoom);
 
394
        if (i != Wavelet.Nbr_Plan-1)
 
395
        for (j = 0; j < Nl*Nc; j++)
 
396
        {
 
397
            if (Imag [j] > Sigma) Imag [j] = Sigma;
 
398
            if (Imag [j] < - Sigma) Imag [j] = - Sigma;
 
399
        }
 
400
       lib_caval_pict  (Nl, Nc, Nl1, Nc1, Imag, Pict_Caval, Inc, Mod_Visu);
 
401
       copy_plan_to_pict(Pict_Caval,Nl1,Nc1,ImaSynt,MAXY,MAXX, i);
 
402
       free ((char *) Imag); 
 
403
    }
 
404
 
 
405
   io_write_pict_f_to_file (File_Name_Imag, ImaSynt, MAXY, MAXX);
 
406
 
 
407
   /* End */
 
408
 return  SCSEPI();
 
409
}
 
410
 
 
411
 
 
412
 
 
413
 
 
414