source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/imageop.cc@ 2830

Last change on this file since 2830 was 658, checked in by ansari, 26 years ago

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