source: Sophya/trunk/SophyaLib/NTools/imageop.cc@ 1009

Last change on this file since 1009 was 832, checked in by ansari, 26 years ago

pb compil linux int k ds boucle rz+cmv 6/4/00
File size: 41.1 KB
Line 
1#include "machdefs.h"
2#include <stdio.h>
3#include <stdlib.h>
4
5#include "simplesort.h"
6#include "nbmath.h"
7#include "histos.h"
8#include "datatypes.h"
9#include "imageop.h"
10#include "nbtri.h"
11
12
13//++
14// Module Operations sur les images (C++)
15// Lib Images++
16// include imageop.h
17//
18// Operations sur les images.
19//--
20
21
22// Outils prives.
23template <class T>
24float RzCmvSigmaCiel(T const * pix,int nx,int ny,int ndata,int bin,
25 int taille,float frac,int deg,
26 float low,float high,int deb);
27
28#ifdef IMGRGCHECK
29#define CHECKIMG2(_img_,_x_,_y_) \
30 if ((_x_ >= (_img_).siz_x) || (_y_ >= (_img_).siz_y) || \
31 (_x_ < 0) || (_y_ < 0)) THROW(rangeCheckErr); \
32 if (!(_img_).ImagePtr()) THROW(nullPtrErr)
33#else
34#if defined(IMGVOIDCHECK)
35#define CHECKIMG2(_img_,_x_,_y_) \
36 if (!(_img_).ImagePtr()) THROW(NullPtrErr)
37#else
38#define CHECKIMG2(_img_,_x_,_y_)
39#endif
40#endif /* IMGRGCHECK */
41
42// Les fonctions
43
44/* Nouvelle-Fonction */
45
46template <class T>
47void CleanImagesOld(int nImg, Image<T>** imgTab, double nSig,
48 double lowbad, double highbad)
49
50/* Nettoie une serie d'image a partir de l'histo de la */
51/* serie de valeurs pour chaque pixel (Voir E. Aubourg) */
52
53{
54 if (lowbad < MinRange((T)0))
55 lowbad = MinRange((T)0);
56
57 if (highbad > MaxRange((T)0))
58 highbad = MaxRange((T)0);
59
60 if (lowbad >= highbad)
61 {lowbad = MinRange((T)0); highbad = MaxRange((T)0);}
62 Histo h(lowbad, highbad, 500);
63
64 for (int i=0; i<imgTab[0]->XSize(); i++)
65 for (int j=0; j<imgTab[0]->YSize(); j++) {
66 h.Zero();
67 int k;
68 for (k=0; k<nImg; k++)
69 h.Add((*imgTab[k])(i,j));
70 float x0, s0=0;
71 x0 = h.CleanedMean(s0);
72 float low = x0 - s0*nSig; if (low < lowbad) low = lowbad;
73 float hi = x0 + s0*nSig; // if (hi > highbad) hi = highbad;
74 for (k=0; k<nImg; k++)
75 if ((*imgTab[k])(i,j) > hi || (*imgTab[k])(i,j) < low)
76 (*imgTab[k])(i,j) = 0;
77 }
78}
79
80/* Nouvelle-Fonction */
81
82template <class T>
83//++
84void CleanImages(int nImg, Image<T>** imgTab, double nSig,
85 double lowbad, double highbad)
86//
87// Nettoie une serie de `nImg' images `imgTab' a partir
88// d'une moyenne tronquee a `nSig' sigma
89// de la serie de valeurs pour chaque pixel.
90// Les valeurs des pixels sont prises entre `lowbad' et `highbad'.
91//--
92{
93 if (lowbad < MinRange((T)0))
94 lowbad = MinRange((T)0);
95
96 if (highbad > MaxRange((T)0))
97 highbad = MaxRange((T)0);
98
99 if (lowbad >= highbad)
100 {lowbad = MinRange((T)0); highbad = MaxRange((T)0);}
101
102 double sx, sx2, sig, mean;
103 int n;
104
105 for (int i=0; i<imgTab[0]->XSize(); i++)
106 for (int j=0; j<imgTab[0]->YSize(); j++) {
107 sig = 1e20; mean = 0;
108 for (int l = 0; l<3; l++) {
109 double low = mean - nSig*sig;
110 if (low < lowbad) low = lowbad;
111 double high = mean + nSig*sig;
112 if (high > highbad) high = highbad;
113 sx = sx2 = n = 0;
114 for (int k=0; k<nImg; k++) {
115 double p = (*imgTab[k])(i,j);
116 if (p > low && p < high) {
117 sx += p;
118 sx2 += p*p;
119 n++;
120 }
121 }
122 if (n>0) {
123 mean = sx / n;
124 sig = sqrt( sx2/n - (sx/n)*(sx/n) );
125 } // sinon on laisse tel quel, ca sera coupe...
126 }
127 double low = mean - nSig*sig;
128 if (low < lowbad) low = lowbad;
129 double high = mean + nSig*sig;
130 if (high > highbad) high = highbad;
131 for (int k=0; k<nImg; k++)
132 if ((*imgTab[k])(i,j) > high || (*imgTab[k])(i,j) < low)
133 (*imgTab[k])(i,j) = 0;
134 }
135}
136
137/* Nouvelle-Fonction */
138
139template <class T>
140//++
141float ImgFondCielHisto(Image<T> const& img, float binWidth, int degree, float frac,
142 float low, float high, int debug)
143//
144// Calcul du fond de ciel par histogrammation. L'histogramme
145// va de `low' a `high' avec une largeur de bin `binWidth'.
146// La valeur la plus probalbe est fittee par un polynome de degre
147// `degree' en utilisant les bins entourant le maximum
148// a `frac' pourcent de sa hauteur.
149//--
150{
151 // CHECKIMG2(img,0,0);
152
153 if (degree < 2 || degree > 3)
154 degree = 3;
155
156 if (frac <= 0 || frac >= 1)
157 frac = 0.5;
158
159 if (low < MinRange((T)0)+2)
160 low = MinRange((T)0)+2;
161
162 if (high > MaxRange((T)0)-2)
163 high = MaxRange((T)0)-2;
164
165 if (low >= high)
166 {low = MinRange((T)0); high = MaxRange((T)0);}
167
168// verifier binwidth
169
170// if ( binWidth < 65536.0/img.XSize()*img.YSize() )
171// binWidth = 65536.0/img.XSize()*img.YSize();
172
173 if (binWidth<0)
174 binWidth = 1;
175
176 if (debug > 1)
177 cout << __FUNCTION__ << " binWidth =" << binWidth <<
178 ", degree = " << degree <<
179 ", frac = " << frac <<
180 ", low = " << low <<
181 ", high = " << high << endl;
182
183 int nBins = int(((double)high-(double)low)/(double)binWidth + 1);
184 Histo hist(low, high, nBins);
185
186 int nPix = img.XSize()*img.YSize();
187 for (int i=0; i<nPix; i++)
188 hist.Add(img[i]);
189
190 hist(0) = 0;
191
192 float fd = hist.FitMax(degree, frac, debug);
193
194 (float&) img.fond = fd;
195
196 return fd;
197}
198
199
200/* Nouvelle-Fonction */
201
202template <class T>
203float ImgSigmaCiel(Image<T> const& img, int ndata, int binWidth,
204 int taille, float frac, int degree,
205 float low, float high, int debug)
206{
207 (float&) img.sigmaFond = RzCmvSigmaCiel(img.ImagePtr(), img.XSize(), img.YSize(),
208 ndata, binWidth, taille, frac,
209 degree, low, high, debug);
210 return img.sigmaFond;
211}
212
213
214/* Nouvelle-Fonction */
215
216template <class T>
217//++
218float RzCmvSigmaCiel(T const* pix,int nx,int ny,int ndata,int bin,
219 int taille,float frac,int deg,
220 float low,float high,int deb)
221//
222// Calcul du fond de ciel par histogrammation
223//
224//| INPUT:
225//| - pix : tableau pix[nx][ny]
226//| - nx,ny : dimensions du tableau pix
227//| - ndata : nombre de donnees a utiliser
228//| - bin : largeur du bin en ADU
229//| - taille: 1/2 dimension de l'histogramme en ADU
230//| - frac : fraction du maximum au dessus de laquelle on fite
231//| - deg : degre fit polynome sous la gaussienne (0 ou 2)
232//| - low : valeur basse limite en ADU
233//| - high : valeur haute limite en ADU
234//| - deb : flag de debug
235//--
236//++
237//| OUTPUT:
238//| - retourne : le sigma du fond de ciel si OK
239//| -1 si probleme definition
240//| -2 si probleme allocation
241//| -3 si probleme remplissage histogramme
242//| -4 si probleme recherche maximum
243//| -5 si probleme calcul fraction du maximum
244//| -6 si probleme nombre de bin
245//| -7 si probleme nombre de fit
246//|
247//--
248// // EA seul changement par rapport a la version C :
249// // pix est un T*, en template. -> C++ necessaire.
250{
251FLOATDOUBLE X,Y,EY;
252double (*Fun) (double ,double *,double *);
253int_4 npar;
254double par[6], epar[6],minpar[6],maxpar[6],stepar[6],stochi2;
255int_4 i,j,i1,j1,j2,n,ip,nbin,inc,ibin0,max,imax,lmax,b1,b2;
256float *histo,*xbin,*ehisto,fwhm,rc;
257double vp,pxv;
258
259/* set reasonable parameters */
260Fun = Gauss1DPolF ;
261if ( deg<=1 ) deg=0; else deg=2;
262low = ( low < 0 ) ? 0 : low;
263high = ( high > 65530 ) ? 65530 : high;
264if ( low >= high ) { low=0; high = 65530;}
265bin = ( bin <= 0 ) ? 1 : bin;
266nbin = int(2.*taille /bin + 1);
267if ( nbin < 20 ) nbin = 20;
268if ( nbin >= 65530 ) nbin = 65530;
269ibin0 = int(-nbin*bin/2.);
270inc = int(nx*ny / ndata); inc = int(sqrt( (double) inc ));
271inc = (inc<1) ? 1 : inc;
272if(deb>1)
273 { printf("Sigma_Ciel: pix[%d,%d] ndata=%d bin=%d taille=%d"
274 ,nx,ny,ndata,bin,taille);
275 printf(" frac=%f low=%g high=%g nbin=%d ibin0=%d inc=%d\n"
276 ,frac,low,high,nbin,ibin0,inc); }
277if ( nx < 5 || ny < 5 ) return (-1.);
278
279/* allocate place for histogramme */
280if ( (histo = (float *) calloc(nbin, sizeof(float)) ) == NULL )
281 return(-2.);
282if ( (xbin = (float *) calloc(nbin, sizeof(float)) ) == NULL )
283 {free(histo); return(-2.);}
284if ( (ehisto = (float *) calloc(nbin, sizeof(float)) ) == NULL )
285 {free(histo); free(xbin); return(-2.);}
286
287/* fill histogramme */
288ndata=0;
289for (j=1;j<ny-1;j=j+inc) { for (i=1;i<nx-1;i=i+inc) {
290 ip = (int)(*(pix+j*nx+i));
291 if ( ip < low || ip > high ) continue;
292 vp = (double) ip;
293 pxv = 0.; n=0;
294 for(j1=j-1;j1<=j+1;j1++) { for(i1=i-1;i1<=i+1;i1++) {
295 if ( i1 == i && j1 == j ) continue;
296 ip = (int)(*(pix+j1*nx+i1));
297 if ( ip < low || ip > high ) continue;
298 pxv += (double) ip;
299 n++;
300 } }
301 if ( n != 8 ) continue;
302 vp = vp - pxv/n;
303 ip = int((vp-ibin0)/bin);
304 if( ip<0 || ip >= nbin ) continue;
305 histo[ip]++;
306 ndata++;
307} }
308for(i=0;i<nbin;i++)
309 { xbin[i] = (float) ibin0 + ( (float) i + 0.5 ) * (float) bin;
310 ehisto[i] = (histo[i]<1.) ? -1. : sqrt((double) histo[i]);}
311
312if(deb>1) printf(" nombre de data corrects = %d\n",ndata);
313if( ndata < 1 ) {rc= -3.; goto END;}
314
315/* find maximum */
316max = -1; imax = -1;
317for(i=nbin-1;i>=0;i--) if(histo[i]>max) { max=int(histo[i]); imax=i;}
318if(deb>1) printf(" maxi histo = %d (bin=%d)\n",max,imax);
319if ( imax == -1 || imax == 0 || imax == nbin-1 ) {rc= -4.; goto END;}
320if(deb>2) {
321 for(i=imax-10;i<=imax+10;i++) if(i>=0 && i<nbin) printf(" %5f",xbin[i]);
322 printf("\n");
323 for(i=imax-10;i<=imax+10;i++) if(i>=0 && i<nbin) printf(" %5d",(int) histo[i]);
324 printf("\n"); }
325
326/* find limites a half maximum */
327lmax = (int) ( (float) max * 0.5 );
328for(b1=imax;b1>0;b1--) if(histo[b1]<=lmax) break;
329for(b2=imax;b2<nbin;b2++) if(histo[b2]<=lmax) break;
330fwhm = (xbin[b2]-xbin[b1]);
331if(deb>1) printf(" fwhm = %f (%d,%d)\n",fwhm,b1,b2);
332
333/* find limites a frac*max */
334lmax = (int) ( (float) max * frac ) + 1;
335if( lmax < 1 )
336 {if(deb>1) printf(" lmax=%d\n",lmax); rc= -5.; goto END;}
337for(b1=imax;b1>0;b1--) if(histo[b1]<=lmax) break;
338for(b2=imax;b2<nbin;b2++) if(histo[b2]<=lmax) break;
339
340n=0;
341for(i=b1;i<=b2;i++) if( histo[i] > 0 ) n++;
342if(deb>1) printf(" limites bin %d %d (v=%d) bin non nuls=%d\n"
343 ,b1,b2,lmax,n);
344/* pas assez de bin, on en ajoute artificiellement */
345/* 8 est une borne large pour fiter gauss + parab (3+3) */
346if( n < 8 ) {
347 j1 = b1; j2 = b2;
348 for(i=1;i<nbin;i++) {
349 if ( b1-i >= 0 ) j1=b1-i;
350 if ( b2+i < nbin ) j2=b2+i;
351 if ( histo[j1] > 0 ) n++;
352 if ( histo[j2] > 0 ) n++;
353 if( n >= 6 ) break;
354 }
355 b1 = j1; b2 = j2;
356 if(deb>1) printf(" Re-update des limites bin %d %d (v=%d) bin non nuls=%d\n"
357 ,b1,b2,lmax,n);
358}
359if( n < 4 ) {rc= -6.; goto END;}
360ndata = n;
361
362/* fit dans les limites */
363npar = 3+deg+1; Set_FitFunDPol(deg);
364if(ndata<8) {npar=4; Set_FitFunDPol(0);}
365if(ndata<5) {npar=3; Set_FitFunDPol(-1);}
366X.fx = xbin;
367Y.fx = histo;
368EY.fx = ehisto;
369stochi2 = 0.0001;
370n = nbin;
371par[0] = max;
372par[1] = histo[imax-1]*xbin[imax-1]+histo[imax]*xbin[imax]+histo[imax+1]*xbin[imax+1];
373par[1] = par[1]/(histo[imax-1]+histo[imax]+histo[imax+1]);
374par[2] = (fwhm<1.) ? 1. : fwhm/2.36;
375par[3] = par[4] = par[5] = 0.;
376stepar[0] = par[0]/20.; stepar[1] = 1.; stepar[2] = par[2]/5.;
377stepar[3] = 1.; stepar[5] = 0.1;
378stepar[4] = 0.;
379for(i=0;i<6;i++) {minpar[i] = 1.; maxpar[i] = -1.;}
380minpar[2] = 0.1; maxpar[2] = 65530.;
381if(deb>1) printf(" choix de fit a %d parametres b1=%d,b2=%d,ndata=%d)\n"
382 ,npar,b1,b2,ndata);
383
384i = (deb>=3) ? deb-2 : 0 ;
385vp = FitFun(Fun,FITFUN_FLOAT,X,Y,EY,&n,b1,b2,par
386 ,epar,stepar,minpar,maxpar,npar,stochi2,50,i);
387
388if( vp<=0. && vp != -4. ) {rc= -7.; goto END;}
389rc = par[2];
390if(deb>1) printf(" sigma=%f (h=%f centrage=%f) xi2=%f (%f pour %d)\n"
391 ,rc,par[0],par[1],vp,vp/(n-npar),npar);
392
393END:
394free(histo); free(xbin); free(ehisto);
395return(rc);
396}
397
398/* ................................................................................. */
399/* Filtrage ( Convolution ) d' images */
400/* ................................................................................. */
401
402/* Nouvelle-Fonction */
403template <class T>
404//++
405Image<T> * FilterImage(Image<T> const * pim, Image<T> * pom, ImageR4 *matx)
406//
407// Filtrage ( Convolution ) d'une image `pim' par une matrice carree
408// representee par l'ImageR4 `matx'.
409// Le resultat se retrouve dans `pom'.
410// Si `pom' = NULL, elle est cree.
411// `pim' et `pom' doivent avoir la meme taille.
412//--
413{
414
415float vpix, *mxe;
416bool fgcr;
417int i,j,k,l,m,n;
418int ndim, nd2;
419
420fgcr = false;
421
422if ((pim == NULL) || (matx == NULL))
423 { printf(" FilterImage_Erreur : pim ou matx = NULL ! \n");
424 return(NULL); }
425if (matx->XSize() != matx->YSize())
426 { printf("FilterImage_Erreur: Filtrage par matrice carre uniquement (%d,%d)\n",
427 matx->XSize(), matx->YSize()); return(NULL); }
428ndim = matx->XSize();
429if ( (ndim < 3) || ((ndim%2) == 0) || (ndim > 25))
430 { printf("FilterImage_Erreur: Pb dimension matrice (%d) (Impair, entre 3 et 25)\n", ndim);
431 return(NULL); }
432if (pom == NULL)
433 { pom = new Image<T>(pim->XSize(), pim->YSize(), false); fgcr = true; }
434
435if ( (pim->XSize() != pom->XSize()) || (pim->YSize() != pom->YSize()) )
436 { printf("FilterImage_Erreur : Taille pim et pom differentes ! \n") ;
437 if (fgcr) delete pom; return(NULL) ;}
438
439nd2 = ndim/2;
440
441/* Traitement de la partie centrale */
442for(j=nd2; j<(pom->YSize()-nd2); j++)
443 {
444 for(i=nd2; i<(pom->XSize()-nd2); i++)
445 {
446 vpix = 0.; mxe = matx->ImagePtr();
447 for(l=-nd2; l<=nd2; l++)
448 for(k=-nd2; k<=nd2; k++)
449 { vpix += (float)((*pim)((i+k), (j+l))) * (*mxe);
450 mxe++; }
451 (*pom)(i, j) = (T)vpix;
452 }
453 }
454
455/* Traitement bord en Y en haut */
456for(j=0; j<nd2; j++)
457 {
458 for(i=nd2; i<(pom->XSize()-nd2); i++)
459 {
460 vpix = 0.; mxe = matx->ImagePtr();
461 for(l=-nd2; l<=nd2; l++)
462
463 for(k=-nd2; k<=nd2; k++)
464 {
465 m = i+k; n = j+l;
466 if (n < 0) n = 0;
467 vpix += (float)((*pim)(m, n)) * (*mxe);
468 mxe++;
469 }
470 (*pom)(i, j) = (T)vpix;
471 }
472 }
473
474/* Traitement bord en Y en bas */
475for(j=(pom->YSize()-nd2); j<pom->YSize(); j++)
476 {
477 for(i=nd2; i<(pom->XSize()-nd2); i++)
478 {
479 vpix = 0.; mxe = matx->ImagePtr();
480 for(l=-nd2; l<=nd2; l++)
481 for(k=-nd2; k<=nd2; k++)
482 {
483 m = i+k; n = j+l;
484 if (n >= pim->YSize()) n = (pim->YSize()-1);
485 vpix += (float)((*pim)(m, n)) * (*mxe);
486 mxe++;
487 }
488 (*pom)(i, j) = (T)vpix;
489 }
490 }
491
492/* traitement bords en X et coins */
493for(j=0; j<pom->YSize(); j++)
494 {
495 for(i=0; i<nd2; i++)
496 {
497 vpix = 0.; mxe = matx->ImagePtr();
498 for(l=-nd2; l<=nd2; l++)
499 for(k=-nd2; k<=nd2; k++)
500 {
501 m = i+k; n = j+l;
502 if (m < 0) m = 0;
503 if (n < 0) n = 0;
504 if (n >= pim->YSize()) n = (pim->YSize()-1);
505 vpix += (float)((*pim)(m, n)) * (*mxe);
506 mxe++;
507 }
508 (*pom)(i, j) = (T)vpix;
509 }
510 for(i=(pom->XSize()-nd2); i<pom->XSize(); i++)
511 {
512 vpix = 0.; mxe = matx->ImagePtr();
513 for(l=-nd2; l<=nd2; l++)
514 for(k=-nd2; k<=nd2; k++)
515 {
516 m = i+k; n = j+l;
517 if (m >= pim->XSize()) m = (pim->XSize()-1);
518 if (n < 0) n = 0;
519 if (n >= pim->YSize()) n = (pim->YSize()-1);
520 vpix += (float)((*pim)(m, n)) * (*mxe);
521 mxe++;
522 }
523 (*pom)(i, j) = (T)vpix;
524 }
525 }
526
527
528return(pom);
529
530}
531
532/* Nouvelle-Fonction */
533//++
534ImageR4* NoiseFiltImage(ImageR4& img, ImageR4& filtre, DynCCD& dynccd)
535// Calcule une image de bruit, a partir de "img" et de "dynccd"
536// prenant en compte le filtre de convolution "filtre".
537// L'image de bruit renvoyee contient les correlations de bruit entre
538// pixels, du au filtrage (convolution) de l'image "img" par le
539// filtre "filtre".
540//--
541{
542 ImageR4* ImgBrt = NoiseImage(&img,&dynccd);
543
544 int i,j;
545 for (i=0; i<ImgBrt->XSize(); i++)
546 for (j=0; j<ImgBrt->YSize(); j++) (*ImgBrt)(i,j) =(*ImgBrt)(i,j) *(*ImgBrt)(i,j);
547
548 ImageR4 Filter2(filtre.XSize(), filtre.YSize()) ;
549 for (i=0; i<filtre.XSize(); i++)
550 for (j=0; j<filtre.YSize(); j++) (Filter2)(i,j) = filtre(i,j) * filtre(i,j);
551
552 ImageR4* out2;
553 out2 = new ImageR4(ImgBrt->XSize(), ImgBrt->YSize());
554 FilterImage(ImgBrt, out2, &Filter2);
555 for (i=0; i<out2->XSize(); i++)
556 for (j=0; j<out2->YSize(); j++) (*out2)(i,j) =sqrt((*out2)(i,j));
557 delete ImgBrt;
558 return(out2);
559}
560
561
562//////////////////////////////////////////////////////////////////////////////
563/* Nouvelle-Fonction CMV 14/10/97 */
564template <class T>
565//++
566Image<T> * FilterStat(Image<T> const * pim, Image<T> * pom
567 ,int N,int M,float lowbad,float highbad)
568//
569// Filtrage median d'une image `(*pim)'.
570// Le resultat se retrouve dans `(*pom)'.
571// Si pom = NULL l'image de sortie est cree automatiquement.
572// pim et pom doivent avoir la meme taille.
573// Le masque de calcul du filtre median est un rectangle `N'x`M'
574// La dynamique de l'image est `lowbad,highbad'.
575// Si lowbad>=highbad dynamique infinie.
576//--
577{
578if(pim == NULL)
579 {printf(" FilterStat_Erreur : pim = NULL ! \n");
580 return(NULL);}
581if( N<=0 || M<=0 )
582 {printf("FilterStat_Erreur: Filtrage par matrice (%d,%d)\n",N,M);
583 return(NULL);}
584bool fgcr = false;
585if(pom == NULL)
586 {pom = new Image<T>(pim->XSize(), pim->YSize(), false); fgcr = true;}
587if( (pim->XSize() != pom->XSize()) || (pim->YSize() != pom->YSize()) )
588 { printf("FilterImage_Erreur : Taille pim et pom differentes ! \n");
589 if(fgcr) delete pom; return(NULL);}
590
591T low = (T) lowbad;
592T high = (T) highbad;
593if(high>low) fgcr = true; else fgcr = false;
594
595N /= 2; M /= 2;
596double * tab = new double[(2*N+1)*(2*M+1)];
597
598for(int j=0;j<pom->YSize();j++) {
599 for(int i=0;i<pom->XSize();i++) {
600 int n=0;
601 for(int jj=j-M;jj<=j+M;jj++) {
602 if( jj<0 || jj >= pom->YSize() ) continue;
603 for(int ii=i-N;ii<=i+N;ii++) {
604 if( ii<0 || ii >= pom->XSize() ) continue;
605 if( fgcr && ( (*pim)(ii,jj)<=low || high<=(*pim)(ii,jj) ) ) continue;
606 tab[n] = (double) (*pim)(ii,jj);
607 n++;
608 }
609 }
610 if(n==0) (*pom)(i,j) = 0;
611 else if(n==1) (*pom)(i,j) = (T) tab[0];
612 else {
613 HeapSort(n,tab);
614 if(n%2==1) (*pom)(i,j) = (T) tab[(n-1)/2];
615 else (*pom)(i,j) = (T) (tab[n/2]+tab[n/2-1])/2.;
616 }
617 }
618}
619
620delete [] tab;
621return(pom);
622}
623
624//////////////////////////////////////////////////////////////////////////////
625/* Nouvelle-Fonction CMV 7/8/98 */
626template <class T>
627//++
628Image<T> * CSplineImage(Image<T> const & f1,double pxsz,double x0,double y0,int natural)
629//
630// Fonction retournant un pointeur sur une image qui contient
631// le Spline 3 a 2 dimensions de l'images f1 (NULL est retourne si probleme).
632// `pxsz' est la taille du pixel de l'image retournee en considerant
633// que l'image `f1' en entree a une taille de pixel egale a 1.
634// `x0,y0' est le decalage eventuel a appliquer de l'image origine a
635// l'image retournee, ce decalage etant mesure en unites de pixels
636// de l'image origine (1).
637// `natural' a la signification decrite dans la classe `CSpline2'.
638// Attention, la taille de l'image retournee sera `nx/pxsz' et `ny/pxsz'.
639// Pour eviter les problemes d'extrapolation sur les bords de l'image
640// d'origine, il leurs est assigne la valeur du centre du pixel auquel ils
641// appartiennent.
642//--
643{
644//// open for result
645int nx = f1.XSize();
646int ny = f1.YSize();
647float dum;
648dum = nx/pxsz;
649int dx = (int) (dum+0.5);
650dum = ny/pxsz;
651int dy = (int) (dum+0.5);
652cout<<"Taille image en entree: "<<nx<<" "<<ny<<endl
653 <<"Taille image en sortie: "<<dx<<" "<<dy<<" (pix size="<<pxsz<<")"<<endl
654 <<"Decalage en x="<<x0<<" en y="<<y0<<endl;
655if(nx<=1 || ny<=1) return NULL;
656if(dx<=0 || dy<=0) return NULL;
657
658Image<T> *f2 = NULL;
659TRY {
660 f2 = new Image<T>(dx,dy);
661} CATCHALL {
662 cout<<"Impossible d'ouvrir le fichier resultat"<<endl;
663 return NULL;
664} ENDTRY
665
666//// Spline 3 2D ... nx+2 pour tenir compte des effets de bords !
667double *x = new double[nx+2];
668double *y = new double[ny+2];
669double *z = new double[(nx+2)*(ny+2)];
670{
671int k;
672// Centre pixel (0,0) en (0.5,0.5) et on ajoute les coordonnees
673// des bords pour extrapolation sur tout l'intervalle.
674// Les coordonnees vont de [0,nx] et [0,ny] (valeurs extremes comprises).
675x[0] = 0.; x[nx+1] = nx;
676y[0] = 0.; y[ny+1] = ny;
677for(k=1;k<=nx;k++) x[k] = k-0.5;
678for(k=1;k<=ny;k++) y[k] = k-0.5;
679// Remplissage des valeurs avec effet de bord [0] et [nx+1],[ny+1].
680k = 0;
681int ii,jj;
682for(int j=0;j<ny+2;j++) {
683 jj = j-1;
684 if(jj<0) jj=0; else if(jj>=ny) jj=ny-1;
685 for(int i=0;i<nx+2;i++) {
686 ii = i-1;
687 if(ii<0) ii=0; else if(ii>=nx) ii=nx-1;
688 z[k] = (double) f1(ii,jj);
689 k++;
690 }
691}
692}
693
694CSpline2 spl(nx+2,x,ny+2,y,z,natural,true);
695spl.ComputeCSpline();
696
697*f2 = 0.;
698{
699// Attention au decalage possible !
700// Les valeurs vont de [0.5*pxsz,(dx-0.5)*pxsz] idem pour y...
701// a corriger eventuellement du decalage x0,y0.
702double xmin = -pxsz/10000., xmax = nx+pxsz/10000.;
703double ymin = -pxsz/10000., ymax = ny+pxsz/10000.;
704for(int j=0;j<dy;j++) {
705 double yy = (j+0.5)*pxsz - y0;
706 if( yy<=ymin || yy>=ymax ) continue;
707 for(int i=0;i<dx;i++) {
708 double xx = (i+0.5)*pxsz - x0;
709 if( xx<=xmin || xx>=xmax ) continue;
710 (*f2)(i,j) = (T) spl.CSplineInt(xx,yy);
711}}}
712
713//// nettoyage
714delete [] x;
715delete [] y;
716delete [] z;
717
718return f2;
719}
720//////////////////////////////////////////////////////////////////////////////
721
722
723
724/* Nouvelle-Fonction */
725int FondSigmaCiel(RzImage *pim, float low, float high,int pvsz
726 ,float nbsig1,float nbsig2,float nbsig3
727 ,float binsg,float frac,int modu,int deb)
728{
729int rc;
730
731 switch (pim->PixelType())
732 {
733 case kuint_2:
734 {
735 ImageU2 imt((RzImage&)(*pim));
736 rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
737 pim->fond = imt.fond;
738 pim->sigmaFond = imt.sigmaFond;
739 break;
740 }
741 case kint_2:
742 {
743 ImageI2 imt((RzImage&)(*pim));
744 rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
745 pim->fond = imt.fond;
746 pim->sigmaFond = imt.sigmaFond;
747 break;
748 }
749 case kint_4:
750 {
751 ImageI4 imt((RzImage&)(*pim));
752 rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
753 pim->fond = imt.fond;
754 pim->sigmaFond = imt.sigmaFond;
755 break;
756 }
757 case kr_4:
758 {
759 ImageR4 imt((RzImage&)(*pim));
760 rc = imt.FondSigmaCiel(low, high, pvsz, nbsig1, nbsig2, nbsig3, binsg, frac, modu, deb);
761 pim->fond = imt.fond;
762 pim->sigmaFond = imt.sigmaFond;
763 break;
764 }
765 default:
766 {
767 cout<<"FindSigmaCiel: type d'image non-prevu "<<(int) pim->PixelType()<<endl;
768 exit(-1);
769 }
770 }
771
772return(rc);
773}
774
775/* Nouvelle-Fonction */
776int CalFondImage(RzImage *pim, double low, double high, int action)
777
778/* ATTENTION: Adaptation a partir d'une version C */
779/* CalculFond() de imgstat.c */
780
781/* Cette fonction calcule le fond et le sigma du ciel a la CMV */
782/* CMV / E.Aubourg / Reza */
783/* action Bit 0 CheckDyn() */
784/* Bit 1 ImgFondCielHisto() */
785/* Bit 2 ImgSigmaCiel */
786/* Bit 3 FondSig2Ciel() */
787/* Bit 1+2+3 (=14) Fond_Ciel() + Sigma_Ciel() Si OK */
788/* FondSig2Ciel() sinon */
789/* action = 0 ou >= 15 CheckDyn + Bit 1+2+3 (=14) */
790/* Rc=0 OK Erreur Sinon */
791/* Sortie : Mise a jour de la struct RzImage */
792
793{
794int n,ok,rc,fga;
795float fnd, sg;
796
797int nbsig, vitesse, nbin;
798float minaff, maxaff, minpix, maxpix;
799float fracmax,fracnul,fracsat;
800char *meth;
801
802#define DBG 0
803
804if (action >= 15) action = 15;
805if (action <= 0) action = 15;
806fga = action;
807fnd = sg = -1.;
808ok = 1;
809
810if (action & 1)
811 { pim->CheckDyn(low, high);
812 action--; }
813
814n = pim->XSize()*pim->YSize();
815
816if (action & 6)
817 {
818 int npu = (n/4 < 5000) ? 5000 : n/4;
819 meth = (char *)"Fond/Sig";
820 switch (pim->PixelType())
821 {
822 case kuint_2:
823 {
824 ImageU2 imt((RzImage&)(*pim));
825 if (action & 2)
826 fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
827 if (action & 4)
828 { sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2,
829 (float)low, (float)high, DBG);
830 if (sg < -1.e-3 ) ok = 0; }
831 break;
832 }
833 case kint_2:
834 {
835 ImageI2 imt((RzImage&)(*pim));
836 if (action & 2)
837 fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
838 if (action & 4)
839 { sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2,
840 (float)low, (float)high, DBG);
841 if (sg < -1.e-3 ) ok = 0; }
842 break;
843 }
844 case kint_4:
845 {
846 ImageI4 imt((RzImage&)(*pim));
847 if (action & 2)
848 fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
849 if (action & 4)
850 { sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2,
851 (float)low, (float)high, DBG);
852 if (sg < -1.e-3 ) ok = 0; }
853 break;
854 }
855 case kr_4:
856 {
857 ImageR4 imt((RzImage&)(*pim));
858 if (action & 2)
859 fnd = ImgFondCielHisto(imt, 1.0f, 2, 0.5f, (float)low, (float)high, DBG);
860 if (action & 4)
861 { sg = ImgSigmaCiel(imt, npu, 1, 200, 0.1f, 2,
862 (float)low,(float) high, DBG);
863 if (sg < -1.e-3 ) ok = 0; }
864 break;
865 }
866
867 case kuint_4:
868 case kr_8:
869 cerr << "CalFondImage(RzImage *, ...) Unsupported image type (kuint_4/kr_8) "
870 << (int) pim->PixelType() << "\n" ;
871 break;
872 default:
873 cerr << "CalFondImage(RzImage *, ...) Unknown image type !! " << (int) pim->PixelType() << "\n" ;
874 break;
875 }
876
877 }
878
879if ( (action == 8) || (!ok && (action & 8)) )
880 {
881 fracmax = 0.85; nbsig = 1; vitesse = 3;
882 minpix = low; maxpix = high; nbin = 10000;
883 rc = RzFondSig2Ciel(pim, vitesse, fracmax, nbsig, nbin,
884 &minpix, &maxpix, &fnd, &sg,
885 &minaff, &maxaff, &fracnul, &fracsat);
886 if (rc != 0) ok = 0;
887 meth = (char *)"FondSig2";
888 }
889
890if (ok > 0) /* Mise a jour */
891 {
892 if (fnd > -1.) pim->fond = fnd;
893 if (sg > -1. ) pim->sigmaFond = sg;
894 printf("CalculFond_Info: Meth=%s (%d) Fond/Sig= %5g %4g Nul/Sat=%d %d\n",
895 meth, fga, fnd, sg, pim->nbNul, pim->nbSat);
896 return(0);
897 }
898else
899 {
900 printf("CalculFond_Erreur: Action=%d NbNul/Sat=%d %d\n",fga,
901 pim->nbNul, pim->nbSat);
902 return(1);
903 }
904
905}
906
907
908
909/* Nouvelle-Fonction */
910
911int RzFondSig2Ciel(RzImage * pim, int vitesse, float FracMax, int NbSigmas, int NBin,
912 float *MinPix, float *MaxPix, float *Fond, float *SigFond,
913 float *MinAff, float *MaxAff, float *FracNul, float *FracSature)
914
915/* ATTENTION : Cette fonction a ete mis en C++, pour traiter des */
916/* RzImage a partir de la fonction FondSig2Ciel() */
917/* de imgstat.c */
918/* R. Ansari 04/95 */
919
920/* Cette fonction calcule les valeurs suivantes pour une image */
921/* Arguments d'entree : */
922/* - pim : Pointeur RzImage */
923/* - vitesse : 1 pixel / vitesse pris en compte */
924/* - FracMax : Parametre de calcul (voir ci-dessous) */
925/* - NbSigmas: MinAff = Fond + NbSigmas*SigFond */
926/* **Note : Si NbSigmas < 0 , on prend : */
927/* MinAff = Fond + NbSigmas*SigFond */
928/* et MaxAff = Fond - NbSigmas*SigFond */
929/* Dans ce cas la gamme d'affichage fait ressortir le fond. */
930/* - NBin : Nb de bins de l'histo */
931/* - La valeur mini des pixels et la valeur de saturation */
932/* (ValNul et ValSat) doivent etre fournies ds Min/MaxPix */
933/* sorties : */
934/* - MinPix,MaxPix = Valeurs Min et Max des pixels (<ValSat) */
935/* - Fond = Pic de l'histo des pixels */
936/* - SigFond = sigma du fond calcule par la relation */
937/* NbPixel(Fond+SigFond) = Pic histo * 0.606 (=Exp(-0.5)) */
938/* - MinAff = Bruit + NbSigmas * Sigma bruit */
939/* - MaxAff est calcule tel que la relation suivante soit verifie */
940/* NbPixel(MinAff<Pixel<MaxAff) = Nb de pixel entre Min-MaxAff */
941/* = FracMax * NbPixel(MinAff<Pixel<ValSat */
942/* - FracNul = Nb pixel < ValNul / Nb total pixel */
943/* - FracSat = Nb pixel satures (>ValSat) / Nb total pixel */
944
945/* return code : 0 OK , 1 Erreur calloc */
946/* 2 : Pic histos non trouve */
947/* 3,4 : Probleme de Min/MaxAff avec ValSat */
948/* 5 : Gamme dynamique minimale non respectee */
949/* (MaxAff - MinAff > MINDYNAMIC */
950
951{
952
953#define MINNBIN 100
954#define MAXNBIN 32765
955
956int nbin;
957
958float minhis = 0., maxhis = 65530. , binwidth;
959int xdim, ydim;
960int_4 under,over,tothis;
961int_4 *phis;
962float ValNul,ValSat;
963int_4 *q;
964
965float val;
966int i;
967float min, max;
968int_4 nbpix,nbpix2, imin, imax;
969int bin,binmax;
970int sigl,sigr;
971float fracfond,fracsignal;
972
973float MINDYNAMIC = 5; /* Gamme dynamique minimale (MaxAff-MinAff)*/
974
975#define FRACPIC 0.6065 /* Valeur de exp(-0.5) */
976
977xdim = pim->XSize();
978ydim = pim->YSize();
979
980ValNul = *MinPix;
981ValSat = *MaxPix;
982if ((ValSat < ValNul+0.01)) ValSat = ValNul+0.01;
983
984nbin = NBin;
985if (nbin < MINNBIN) nbin = MINNBIN;
986if (nbin > MAXNBIN) nbin = MAXNBIN;
987
988if (vitesse < 1) vitesse = 1;
989if (vitesse > xdim*ydim/100) vitesse = xdim*ydim/100;
990
991
992/* Je modifie min et max histo et nbin en fonction de valnul et valsat */
993/* je garde binwidth constant */
994
995binwidth = (ValSat-ValNul) / nbin;
996minhis = ValNul;
997maxhis = minhis + nbin*binwidth;
998
999
1000(*MinAff) = minhis;
1001(*MaxAff) = maxhis;
1002*MinPix = 0.;
1003*MaxPix = -1.;
1004*Fond = -1.;
1005*SigFond = -1.;
1006*FracNul = 0.0;
1007*FracSature = 0.0;
1008
1009/* allocation dynamique de memoire */
1010phis = new int_4[nbin];
1011
1012/* mise a zero de l'histogramme */
1013for ( q=phis ; q<phis+nbin ; *q++ = 0 ) ;
1014under = over = tothis = 0;
1015
1016
1017
1018nbpix2 = nbpix = 0;
1019min = ValSat+1.;
1020max = ValNul-1;
1021
1022/* Boucle sur les pixels */
1023/* Pour aller plus vite, on prend un pixel sur 7 */
1024for ( i = 0 ; i < xdim*ydim ; i += vitesse )
1025 {
1026 val = pim->IValue(i);
1027
1028 /* Calcul PixMax/Min et NbPix sature */
1029 if (val > ValSat) nbpix++;
1030 else {
1031 if(val < ValNul) nbpix2++;
1032 else {
1033 if (val > max) max = val;
1034 else if (val < min) min = val;
1035 }
1036 }
1037
1038 /* remplissage de l'histogramme */
1039 bin = (int)((val - minhis) / binwidth) ;
1040 if ( bin < 0 ) under++;
1041 else
1042 {
1043 if (bin < nbin) { (*(phis+bin))++; tothis++; }
1044 else over++;
1045 }
1046 } /* Fin de for : Boucle sur les pixels */
1047
1048(*MinPix) = min;
1049(*MaxPix) = max;
1050/* Ne pas oublier le facteur vitesse pour le calcul FracSatue */
1051(*FracNul) = (float)(nbpix2 * vitesse) / (float)(xdim*ydim);
1052(*FracSature) = (float)(nbpix * vitesse) / (float)(xdim*ydim);
1053
1054/* On verifie qu'il y a eu une entree ds l'histo */
1055if (tothis <= 0)
1056 { delete[] phis; return(2); }
1057
1058/* recherche du maximum */
1059binmax = -1 ; imax = -1 ; q = phis;
1060for ( i = 0 ; i < nbin ; i++, q++ )
1061 if ( *(q) > imax ) {imax = *(q) ; binmax = i ;}
1062
1063if ( (binmax < 0) || (binmax >= nbin) )
1064 { delete[] phis; return(2); }
1065
1066*Fond = minhis + binmax * binwidth;
1067
1068/* On va chercher le sigma du fond */
1069imin = (int_4) ( (float)(imax) * FRACPIC ) ;
1070
1071/*
1072printf("StatsDebug : min/maxhis %g %g (nbin= %d) w=%g\n",
1073 minhis,maxhis,nbin,binwidth);
1074printf(".binmax= %d imax(=pic) = %d , imin= %d \n ",binmax,imax,imin);
1075*/
1076
1077nbpix = 0;
1078
1079sigr = 0;
1080q = phis+binmax;
1081while ( ((*q) > imin) && (sigr < (nbin-binmax)) )
1082 { sigr++; nbpix += (*q) ; q++; }
1083
1084sigl = 0;
1085q = phis+binmax;
1086while ( ((*q) > imin) && (sigl < binmax) )
1087 { sigl++; nbpix += (*q) ; q--; }
1088nbpix -= (*(phis+binmax));
1089
1090(*SigFond) = (float)(sigl+sigr)*binwidth/2.;
1091fracfond = (float)(nbpix) / (float)(tothis);
1092
1093(*MinAff) = minhis+(binmax-sigl)*binwidth;
1094(*MaxAff) = minhis+(binmax+sigr)*binwidth;
1095
1096
1097/* Check sur NbSigmas : */
1098if (NbSigmas == 0) NbSigmas = 1;
1099if (NbSigmas < 0) /* On considere le cas ou NbSigmas est negatif */
1100 {
1101 if ( (bin = binmax+NbSigmas*sigr) < 0 ) bin = 0;
1102 (*MinAff) = minhis + (float)bin*binwidth;
1103 if ( (bin = binmax-NbSigmas*sigr) >= nbin ) bin = nbin-1;
1104 (*MaxAff) = minhis + (float)bin*binwidth;
1105 goto sortie;
1106 }
1107
1108if ( (bin = binmax+NbSigmas*sigl) >= nbin )
1109 { delete[] phis; return(4); }
1110
1111(*MinAff) = minhis + bin*binwidth;
1112
1113/* Calcul du nb de pixel entre MinAff et ValSat */
1114nbpix = 0;
1115for( q = phis+bin ; q < phis+nbin ; q++) nbpix += (*q);
1116
1117fracsignal = (float)(nbpix) / (float)(tothis+over+under);
1118nbpix = (int_4)((float)nbpix * FracMax);
1119nbpix2 = 0; q = phis+bin; i = bin;
1120while( (nbpix2 < nbpix) && (i < nbin) )
1121 { nbpix2 += (*q); i++; q++; }
1122
1123(*MaxAff) = minhis + (float)i*binwidth;
1124
1125
1126sortie: /* On sort en liberant l'espace */
1127
1128/*
1129printf(" Stats..Debug : Fond,sigl,sigr,SigFond = %g,%d,%d,%g \n",
1130(*Fond),sigl,sigr,(*SigFond));
1131 printf(" Min/MaxPix= %g %g , Min/MaxAff= %g %g \n",(*MinPix),
1132(*MaxPix),(*MinAff),(*MaxAff));
1133printf(" FracFond,FracSignal %g %g FracNul,FracSat %g %g\n",
1134 fracfond,fracsignal,(*FracNul),(*FracSature));
1135*/
1136
1137delete[] phis;
1138if ( (*MaxAff) < ((*MinAff)+MINDYNAMIC) ) return(5);
1139return(0);
1140
1141}
1142
1143
1144// ********** INSTANCES
1145
1146#if defined(__xlC) || defined(__aCC__)
1147void instancetempimageop(int n)
1148{
1149/* Cette fonction sert uniquement a forcer le compilo a instancier les
1150 classes/fonctions template */
1151
1152int m = n-50;
1153ImageU2 iu2(n,n), *tiu2[5], xiu2(n,n);
1154ImageI2 ii2(n,n), *tii2[5], xii2(n,n);
1155ImageI4 ii4(n,n), *tii4[5], xii4(n,n);
1156ImageR4 ir4(n,n), *tir4[5], xir4(n,n);
1157
1158ImageR4 matx(5,5);
1159
1160RzCmvSigmaCiel(iu2.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
1161RzCmvSigmaCiel(ii2.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
1162RzCmvSigmaCiel(ii4.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
1163RzCmvSigmaCiel(ir4.ImagePtr(), n, n, 100, 50, 50, (float)0.5, 2, (float)0., (float)65000., 0);
1164
1165ImgFondCielHisto(iu2);
1166ImgFondCielHisto(ii2);
1167ImgFondCielHisto(ii4);
1168ImgFondCielHisto(ir4);
1169
1170ImgSigmaCiel(iu2, 100, 1, 100);
1171ImgSigmaCiel(ii2, 100, 1, 100);
1172ImgSigmaCiel(ii4, 100, 1, 100);
1173ImgSigmaCiel(ir4, 100, 1, 100);
1174
1175CleanImages(5, tiu2);
1176CleanImages(5, tii2);
1177CleanImages(5, tii4);
1178CleanImages(5, tir4);
1179
1180
1181FilterImage(&iu2, &xiu2, &matx);
1182FilterImage(&ii2, &xii2, &matx);
1183FilterImage(&ii4, &xii4, &matx);
1184FilterImage(&ir4, &xir4, &matx);
1185
1186FilterStat(&iu2, &xiu2,(int)3,(int)3,(float)1.,(float)-1.);
1187FilterStat(&ii2, &xii2,(int)3,(int)3,(float)1.,(float)-1.);
1188FilterStat(&ii4, &xii4,(int)3,(int)3,(float)1.,(float)-1.);
1189FilterStat(&ir4, &xir4,(int)3,(int)3,(float)1.,(float)-1.);
1190
1191CSplineImage(iu2,(double)1,(double)0,(double)0,(int)0);
1192CSplineImage(ii2,(double)1,(double)0,(double)0,(int)0);
1193CSplineImage(ii4,(double)1,(double)0,(double)0,(int)0);
1194CSplineImage(ir4,(double)1,(double)0,(double)0,(int)0);
1195
1196return;
1197}
1198
1199#endif
1200
1201
1202#ifdef __CXX_PRAGMA_TEMPLATES__
1203
1204#pragma define_template RzCmvSigmaCiel<uint_2>
1205#pragma define_template RzCmvSigmaCiel<int_2>
1206#pragma define_template RzCmvSigmaCiel<int_4>
1207#pragma define_template RzCmvSigmaCiel<r_4>
1208
1209#pragma define_template ImgFondCielHisto<uint_2>
1210#pragma define_template ImgFondCielHisto<int_2>
1211#pragma define_template ImgFondCielHisto<int_4>
1212#pragma define_template ImgFondCielHisto<r_4>
1213
1214#pragma define_template ImgSigmaCiel<uint_2>
1215#pragma define_template ImgSigmaCiel<int_2>
1216#pragma define_template ImgSigmaCiel<int_4>
1217#pragma define_template ImgSigmaCiel<r_4>
1218
1219#pragma define_template CleanImages<uint_2>
1220#pragma define_template CleanImages<int_2>
1221#pragma define_template CleanImages<int_4>
1222#pragma define_template CleanImages<r_4>
1223
1224
1225#pragma define_template FilterImage<uint_2>
1226#pragma define_template FilterImage<int_2>
1227#pragma define_template FilterImage<int_4>
1228#pragma define_template FilterImage<r_4>
1229
1230#pragma define_template FilterStat<uint_2>
1231#pragma define_template FilterStat<int_2>
1232#pragma define_template FilterStat<int_4>
1233#pragma define_template FilterStat<r_4>
1234
1235#pragma define_template CSplineImage<uint_2>
1236#pragma define_template CSplineImage<int_2>
1237#pragma define_template CSplineImage<int_4>
1238#pragma define_template CSplineImage<r_4>
1239
1240#endif
1241
1242
1243#if ( defined(ANSI_TEMPLATES) && !defined(__aCC__) )
1244template float RzCmvSigmaCiel<uint_2 const>(uint_2 const*, int, int, int, int, int, float, int, float, float, int);
1245template float RzCmvSigmaCiel< int_2 const>( int_2 const*, int, int, int, int, int, float, int, float, float, int);
1246template float RzCmvSigmaCiel< int_4 const>( int_4 const*, int, int, int, int, int, float, int, float, float, int);
1247template float RzCmvSigmaCiel< r_4 const>( r_4 const*, int, int, int, int, int, float, int, float, float, int);
1248template float ImgFondCielHisto<uint_2>(Image<uint_2> const&, float, int, float, float , float , int);
1249template float ImgFondCielHisto< int_2>(Image< int_2> const&, float, int, float, float , float , int);
1250template float ImgFondCielHisto< int_4>(Image< int_4> const&, float, int, float, float , float , int);
1251template float ImgFondCielHisto< r_4>(Image< r_4> const&, float, int, float, float , float , int);
1252template float ImgSigmaCiel<uint_2>(Image<uint_2> const& img, int, int, int, float, int, float , float, int);
1253template float ImgSigmaCiel< int_2>(Image< int_2> const& img, int, int, int, float, int, float , float, int);
1254template float ImgSigmaCiel< int_4>(Image< int_4> const& img, int, int, int, float, int, float , float, int);
1255template float ImgSigmaCiel< r_4>(Image< r_4> const& img, int, int, int, float, int, float , float, int);
1256template void CleanImages<uint_2>(int, Image<uint_2>**, double, double, double);
1257template void CleanImages< int_2>(int, Image< int_2>**, double, double, double);
1258template void CleanImages< int_4>(int, Image< int_4>**, double, double, double);
1259template void CleanImages< r_4>(int, Image< r_4>**, double, double, double);
1260
1261
1262template Image<uint_2>* FilterImage<uint_2>(Image<uint_2> const *, Image<uint_2>*, ImageR4 *);
1263template Image< int_2>* FilterImage< int_2>(Image< int_2> const *, Image< int_2>*, ImageR4 *);
1264template Image< int_4>* FilterImage< int_4>(Image< int_4> const *, Image< int_4>*, ImageR4 *);
1265template Image< r_4>* FilterImage< r_4>(Image< r_4> const *, Image< r_4>*, ImageR4 *);
1266
1267template Image<uint_2>* FilterStat<uint_2>(Image<uint_2> const *, Image<uint_2>*,int,int,float,float);
1268template Image< int_2>* FilterStat< int_2>(Image< int_2> const *, Image< int_2>*,int,int,float,float);
1269template Image< int_4>* FilterStat< int_4>(Image< int_4> const *, Image< int_4>*,int,int,float,float);
1270template Image< r_4>* FilterStat< r_4>(Image< r_4> const *, Image< r_4>*,int,int,float,float);
1271
1272template Image<uint_2>* CSplineImage<uint_2>(Image<uint_2> const &,double,double,double,int);
1273template Image< int_2>* CSplineImage< int_2>(Image< int_2> const &,double,double,double,int);
1274template Image< int_4>* CSplineImage< int_4>(Image< int_4> const &,double,double,double,int);
1275template Image< r_4>* CSplineImage< r_4>(Image< r_4> const &,double,double,double,int);
1276
1277#endif
1278
1279
1280#ifdef GNU_TEMPLATES
1281template float RzCmvSigmaCiel(uint_2 const*, int, int, int, int, int, float, int, float, float, int);
1282template float RzCmvSigmaCiel( int_2 const*, int, int, int, int, int, float, int, float, float, int);
1283template float RzCmvSigmaCiel( int_4 const*, int, int, int, int, int, float, int, float, float, int);
1284template float RzCmvSigmaCiel( r_4 const*, int, int, int, int, int, float, int, float, float, int);
1285template float ImgFondCielHisto(Image<uint_2> const&, float, int, float, float , float , int);
1286template float ImgFondCielHisto(Image< int_2> const&, float, int, float, float , float , int);
1287template float ImgFondCielHisto(Image< int_4> const&, float, int, float, float , float , int);
1288template float ImgFondCielHisto(Image< r_4> const&, float, int, float, float , float , int);
1289template float ImgSigmaCiel(Image<uint_2> const& img, int, int, int, float, int, float , float, int);
1290template float ImgSigmaCiel(Image< int_2> const& img, int, int, int, float, int, float , float, int);
1291template float ImgSigmaCiel(Image< int_4> const& img, int, int, int, float, int, float , float, int);
1292template float ImgSigmaCiel(Image< r_4> const& img, int, int, int, float, int, float , float, int);
1293template void CleanImages(int, Image<uint_2>**, double, double, double);
1294template void CleanImages(int, Image< int_2>**, double, double, double);
1295template void CleanImages(int, Image< int_4>**, double, double, double);
1296template void CleanImages(int, Image< r_4>**, double, double, double);
1297
1298
1299template Image<uint_2>* FilterImage(Image<uint_2> const *, Image<uint_2>*, ImageR4 *);
1300template Image< int_2>* FilterImage(Image< int_2> const *, Image< int_2>*, ImageR4 *);
1301template Image< int_4>* FilterImage(Image< int_4> const *, Image< int_4>*, ImageR4 *);
1302template Image< r_4>* FilterImage(Image< r_4> const *, Image< r_4>*, ImageR4 *);
1303
1304template Image<uint_2>* FilterStat(Image<uint_2> const *, Image<uint_2>*,int,int,float,float);
1305template Image< int_2>* FilterStat(Image< int_2> const *, Image< int_2>*,int,int,float,float);
1306template Image< int_4>* FilterStat(Image< int_4> const *, Image< int_4>*,int,int,float,float);
1307template Image< r_4>* FilterStat(Image< r_4> const *, Image< r_4>*,int,int,float,float);
1308
1309template Image<uint_2>* CSplineImage(Image<uint_2> const &,double,double,double,int);
1310template Image< int_2>* CSplineImage(Image< int_2> const &,double,double,double,int);
1311template Image< int_4>* CSplineImage(Image< int_4> const &,double,double,double,int);
1312template Image< r_4>* CSplineImage(Image< r_4> const &,double,double,double,int);
1313
1314#endif
Note: See TracBrowser for help on using the repository browser.