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

Last change on this file since 3366 was 2615, checked in by cmv, 21 years ago

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