00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include <locale.h>
00029 #include <libintl.h>
00030
00031
00032 #include "mem2d.h"
00033 #include "glbvars.h"
00034 #include "clip.h"
00035
00036 #define UTF8_DEGREE "\302\260"
00037
00038
00039 #define SET_PROGRESS(P) { gapp->SetProgress((gfloat)(P)); while (gtk_events_pending()) gtk_main_iteration(); }
00040
00041
00042 #define SAVECONVOLSRC
00043
00044 char *ZD_name[] = { "I", "Byte", "Short", "Long", "ULong", "LLong",
00045 "Float", "Double", "Complex", "RGBA", "Event"
00046 };
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 ZData::ZData(int Nx, int Ny, int Nv){
00057 XSM_DEBUG (DBG_L6, "ZData::ZData");
00058 nx=Nx; ny=Ny; nv=Nv; vlayer=0;
00059 Li = new LineInfo[ny];
00060 Xlookup=new double[nx];
00061 Ylookup=new double[ny];
00062 Vlookup=new double[nv];
00063 memset (Xlookup, 0, Nx);
00064 memset (Ylookup, 0, Ny);
00065 memset (Vlookup, 0, Nv);
00066 }
00067
00068 ZData::~ZData(){
00069 XSM_DEBUG (DBG_L6, "ZData::~ZData");
00070 delete [] Ylookup;
00071 delete [] Xlookup;
00072 delete [] Vlookup;
00073 delete [] Li;
00074 }
00075
00076 int ZData::ZResize(int Nx, int Ny, int Nv){
00077 XSM_DEBUG (DBG_L6, "ZData: ZResize, delete Li Li=" << Li << " " << Li[0].IsValid());
00078 double *xl = Xlookup;
00079 double *yl = Ylookup;
00080 double *vl = Vlookup;
00081
00082 if (Nv == 0) Nv=nv;
00083
00084 if(Li) delete [] Li;
00085 Li=new LineInfo[Ny];
00086 Xlookup=new double[Nx];
00087 Ylookup=new double[Ny];
00088 Vlookup=new double[Nv];
00089 memset (Xlookup, 0, Nx);
00090 memset (Ylookup, 0, Ny);
00091 memset (Vlookup, 0, Nv);
00092 if(xl){
00093 memcpy(Xlookup, xl, (Nx > nx ? nx:Nx)*sizeof(double));
00094 delete [] xl;
00095 }
00096 if(yl){
00097 memcpy(Ylookup, yl, (Ny > ny ? ny:Ny)*sizeof(double));
00098 delete [] yl;
00099 }
00100 if(vl){
00101 memcpy(Vlookup, vl, (Nv > nv ? nv:Nv)*sizeof(double));
00102 delete [] vl;
00103 }
00104 nx=Nx; ny=Ny; nv=Nv; vlayer=0; cp_ix0=0; cp_num=nx;
00105 XSM_DEBUG (DBG_L6, "ZData: ZResize OK");
00106
00107 return 0;
00108 }
00109
00110 void ZData::ZPutDataSetDest(int ix0, int num)
00111 {
00112 XSM_DEBUG (DBG_L6, "ZData:ZPutDataSetDest");
00113 if((ix0 >= 0) && (num > 0) && (ix0+num) <= nx) {
00114 cp_ix0=ix0;
00115 cp_num=num;
00116 XSM_DEBUG (DBG_L6, "ZData:ZPutDataSetDest => " << cp_ix0 << " #=" << cp_num);
00117 }
00118 };
00119
00120
00121
00122
00123
00124 template <class ZTYP>
00125 TZData<ZTYP>::TZData(int Nx, int Ny, int Nv)
00126 :ZData(Nx, Ny, Nv){ XSM_DEBUG (DBG_L6, "TZData::TZData"); Zdat=NULL; TNew(); }
00127
00128 template <class ZTYP>
00129 TZData<ZTYP>::~TZData(){ XSM_DEBUG (DBG_L6, "TZData::~TZData"); TDel(); }
00130
00131 template <class ZTYP>
00132 int TZData<ZTYP>::Resize(int Nx, int Ny, int Nv){
00133 XSM_DEBUG (DBG_L6, "TZData: Resize(" << Nx << " " << Ny<< ")");
00134 if(Nx != nx || Ny != ny || Nv != nv){
00135 if(Nx == nx && Ny < ny && Nv == nv){
00136 while(ny-- > Ny){
00137 for(int v=0; v<nv; ++v){
00138 delete [] Zdat[ny*nv+v];
00139 Zdat[ny*nv+v] = NULL;
00140 }
00141 }
00142 ny++;
00143 ZResize(Nx, Ny, Nv);
00144 }else{
00145 if(!Zdat){
00146 XSM_DEBUG (DBG_L6, "TZData: Resize Del old");
00147 TDel();
00148 XSM_DEBUG (DBG_L6, "TZData: Resize set size, Info");
00149 ZResize(Nx, Ny, Nv);
00150 XSM_DEBUG (DBG_L6, "TZData: Resize New");
00151 TNew();
00152 }else{
00153 XSM_DEBUG (DBG_L6, "TZData: real Resize, copy old ");
00154 ZTYP** Zdold = Zdat;
00155 while(ny-- > Ny){
00156 for(int v=0; v<nv; ++v){
00157 delete [] Zdold[ny*nv+v];
00158 Zdold[ny*nv+v] = NULL;
00159 }
00160 }
00161 ny++;
00162 Zdat = new ZTYP*[Ny*Nv];
00163 for(int i=0; i<Ny; i++){
00164 for(int v=0; v<Nv; ++v){
00165 Zdat[i*Nv+v] = new ZTYP[Nx];
00166 if(!Zdat[i*Nv+v]){
00167 Ny=i-1;
00168 XSM_DEBUG_ERROR (DBG_L1, "No Mem in Mem2d:New" );
00169 return 0;
00170 }
00171 if(Nx > nx || i >= ny || v >= nv)
00172 memset(Zdat[i*Nv+v], 0, sizeof(ZTYP)*Nx);
00173 if(i < ny && v < nv){
00174 memcpy(Zdat[i*Nv+v], Zdold[i*nv+v], sizeof(ZTYP)*(Nx < nx ? Nx : nx));
00175 delete [] Zdold[i*nv+v];
00176 }
00177 }
00178 }
00179 delete [] Zdold;
00180
00181 ZResize(Nx, Ny, Nv);
00182 }
00183 }
00184 return 1;
00185 }
00186 return 0;
00187 }
00188
00189
00190 template <class ZTYP>
00191 double TZData<ZTYP>::Z(double vx, double vy){
00192 return 0.;
00193 }
00194
00195 template <class ZTYP>
00196 double TZData<ZTYP>::Z(double z, double vx, double vy){
00197 return 0.;
00198 }
00199
00200
00201 template <class ZTYP>
00202 void TZData<ZTYP>::ZPutDataLine(int y, void *src, int mode){
00203 switch(mode){
00204 case MEM_SET:
00205 memcpy((void*)(Zdat[y*nv+vlayer]+cp_ix0), src, sizeof(ZTYP)*cp_num);
00206 break;
00207 case MEM_ADDTO:
00208 for(int i=cp_ix0; i<(cp_ix0+cp_num); i++)
00209 Zdat[y*nv+vlayer][i] += *((ZTYP*)src+i);
00210 break;
00211 case MEM_AVG:
00212 for(int i=cp_ix0; i<(cp_ix0+cp_num); i++){
00213 Zdat[y*nv+vlayer][i] += *((ZTYP*)src+i);
00214 Zdat[y*nv+vlayer][i] /= 2;
00215 }
00216 break;
00217 }
00218 Li[y].invalidate();
00219 }
00220
00221 template <class ZTYP>
00222 int TZData<ZTYP>::CopyFrom(ZData *src, int x, int y, int tox, int toy, int nx, int ny){
00223 for(int i=0; i<ny; i++){
00224 memcpy((void*)&Zdat[(toy++)*nv+vlayer][tox], src->GetPtr(x,y++), nx*sizeof(ZTYP));
00225 Li[i].invalidate();
00226 }
00227 return 0;
00228 }
00229
00230 template <class ZTYP>
00231 void TZData<ZTYP>::TNew( void ){
00232 XSM_DEBUG (DBG_L6, "TZData: TNew: (" << nx << "," << ny << ")");
00233 Zdat = new ZTYP*[ny*nv];
00234 for(int i=0; i<ny*nv; i++){
00235 Zdat[i] = new ZTYP[nx];
00236 if(!Zdat[i]){
00237 ny=i-1;
00238 XSM_DEBUG_ERROR (DBG_L1, "No Mem in Mem2d:New" );
00239 return;
00240 }
00241 memset(Zdat[i], 0, sizeof(ZTYP)*nx);
00242 }
00243 }
00244
00245 template <class ZTYP>
00246 void TZData<ZTYP>::TDel( void ){
00247 if(!Zdat) return;
00248 XSM_DEBUG (DBG_L6, "TZData: TDel ny=" << ny);
00249 for(int i=0; i<ny*nv; i++){
00250 delete [] Zdat[i];
00251 }
00252 delete [] Zdat;
00253 Zdat = NULL;
00254 XSM_DEBUG (DBG_L6, "TZData: TDel OK. Zdat=" << Zdat);
00255 }
00256
00257
00258 template <class ZTYP>
00259 void TZData<ZTYP>::set_all_Z (double z, int v, int x0, int y0, int xs, int ys){
00260 ZTYP zval = (ZTYP)z;
00261 int vs = 0;
00262 int ve = nv-1;
00263 if (v >= 0 && v < nv)
00264 vs=ve=v;
00265 if (!xs) xs = nx;
00266 if (!ys) ys = ny;
00267
00268 for (int i=x0; i < (x0+xs); Zdat[y0*nv+vs][i++] = zval);
00269
00270 for (v = vs; v<=ve; ++v)
00271 for (int i=y0+1; i < (y0+ys); ++i)
00272 memcpy (&Zdat[i*nv+v][x0], &Zdat[y0*nv+vs][x0], xs*sizeof(ZTYP));
00273 }
00274
00275
00276
00277 template <class ZTYP>
00278 void TZData<ZTYP>::NcPut(NcVar *ncfield){
00279 gdouble yy, yy2;
00280 yy = yy2 = 0.;
00281 for(int y=0; y<ny; y++){
00282 for(int v=0; v<nv; ++v){
00283 ncfield->set_cur(0,v,y);
00284 ncfield->put(Zdat[y*nv+v], 1,1, 1,nx);
00285 }
00286 yy = (double)y/(double)ny;
00287 if (ceil(yy) - ceil(yy2) > 1.)
00288 gapp->progress_info_set_bar_fraction (yy2=yy, 2);
00289 }
00290 }
00291
00292 template <class ZTYP>
00293 void TZData<ZTYP>::NcGet(NcVar *ncfield){
00294 gdouble yy, yy2;
00295 yy = yy2 = 0.;
00296 for(int y=0; y<ny; y++){
00297 for(int v=0; v<nv; ++v){
00298 ncfield->set_cur(0,v,y);
00299 ncfield->get(Zdat[y*nv+v], 1,1, 1,nx);
00300 Li[y].invalidate();
00301 }
00302 yy = (double)y/(double)ny;
00303 if (ceil(yy) - ceil(yy2) > 1.)
00304 gapp->progress_info_set_bar_fraction (yy2=yy, 2);
00305 }
00306 }
00307
00308
00309
00310
00311
00312
00313
00314
00315 Mem2d::Mem2d(int Nx, int Ny, ZD_TYPE type){
00316 data = NULL;
00317 scan_event_list = NULL;
00318 SetPlaneMxyz ();
00319 Mnew (Nx, Ny, 1, type);
00320 }
00321
00322 Mem2d::Mem2d(int Nx, int Ny, int Nv, ZD_TYPE type){
00323 data=NULL;
00324 scan_event_list = NULL;
00325 SetPlaneMxyz ();
00326 Mnew (Nx, Ny, Nv, type);
00327 }
00328
00329
00330 Mem2d::Mem2d(Mem2d *m, MEM2D_CREATE_MODE Option){
00331 data=NULL;
00332 scan_event_list = NULL;
00333 GetPlaneMxyz (plane_mx, plane_my, plane_z0);
00334 Mnew(m->GetNx(), m->GetNy(), m->GetNv(), m->GetTyp());
00335 switch(Option){
00336 case M2D_ZERO: break;
00337 case M2D_COPY:
00338 m->data->StoreLayer();
00339 for(int v=0; v<m->GetNv(); ++v){
00340 m->SetLayer(v);
00341 ConvertFrom(m, 0,0, 0,0, GetNx(),GetNy());
00342 }
00343 m->data->RestoreLayer();
00344 break;
00345 }
00346 }
00347
00348
00349 Mem2d::Mem2d(Mem2d *m, int Nx, int Ny, int Nv){
00350 data=NULL;
00351 scan_event_list = NULL;
00352 GetPlaneMxyz (plane_mx, plane_my, plane_z0);
00353 Mnew(Nx, Ny, Nv, m->GetTyp());
00354 }
00355
00356 Mem2d::~Mem2d(){
00357 XSM_DEBUG (DBG_L6, "Mem2d::~Mem2d");
00358 RemoveScanEvents ();
00359 Mdelete();
00360 }
00361
00362 void Mem2d::evl_remove(gpointer entry, gpointer from){
00363 delete (ScanEvent*)entry;
00364 }
00365
00366 void Mem2d::RemoveScanEvents (){
00367 g_slist_foreach(scan_event_list, (GFunc) Mem2d::evl_remove, scan_event_list);
00368 g_slist_free (scan_event_list);
00369 scan_event_list = NULL;
00370 }
00371
00372 gint compare_events_distance (ScanEvent *a, ScanEvent *b, double *xy){
00373 double da = a->distance (xy);
00374 double db = b->distance (xy);
00375 if (da < db) return -1;
00376 if (da > db) return 1;
00377 return 0;
00378 }
00379
00380 GSList* Mem2d::ReportScanEvents (GFunc report_obj_func, gpointer gp, double *xy, double distance, int number){
00381 if (!scan_event_list)
00382 return NULL;
00383 if (!xy && report_obj_func)
00384 g_slist_foreach(scan_event_list, (GFunc)report_obj_func, gp);
00385 else{
00386 GSList *ev_sort = g_slist_copy (scan_event_list);
00387 ev_sort = g_slist_sort_with_data (ev_sort, GCompareDataFunc (compare_events_distance), (gpointer)xy);
00388 if (!report_obj_func)
00389 return ev_sort;
00390
00391 GSList *ev = ev_sort;
00392 int i=0;
00393 while (ev){
00394 if (++i > number || ((ScanEvent*)(ev->data))->distance (xy) > distance)
00395 break;
00396
00397 (*report_obj_func)(ev->data, gp);
00398 ev = g_slist_next (ev);
00399 }
00400 g_slist_free (ev_sort);
00401 return NULL;
00402 }
00403 return NULL;
00404 }
00405
00406 int Mem2d::WriteScanEvents (NcFile *ncf){
00407 GSList *sel;
00408 int p_index;
00409 int u_index;
00410 int p_dim_sets=0;
00411 int p_dim_samples=0;
00412 NcVar* evdata_var=NULL;
00413 NcVar* evcoords_var=NULL;
00414
00415 class Event{
00416 public:
00417 Event (Event *e, ScanEvent *s, ProbeEntry *p, int _j, int _n){
00418 se=s; pe=p; n=_n; j=_j;
00419 if (e)
00420 same_count = same (*e) ? e->same_count+1 : 1;
00421 else
00422 same_count = 1;
00423 };
00424 ~Event (){};
00425 int same (Event &e) { return n == e.n && j == e.j ? TRUE:FALSE; };
00426 ScanEvent *se;
00427 ProbeEntry *pe;
00428 int same_count;
00429 int n;
00430 int j;
00431 };
00432 GSList *ProbeEventsList=NULL;
00433 Event *e=NULL;
00434
00435 u_index=0;
00436
00437 sel = scan_event_list;
00438 while (sel){
00439 ScanEvent *se = (ScanEvent*)sel->data;
00440 GSList *eel = se->event_list;
00441 while (eel){
00442 EventEntry *ee = (EventEntry*)eel->data;
00443 switch (ee->description_id ()){
00444 case 'P':
00445 {
00446 ProbeEntry *pe = (ProbeEntry*)eel->data;
00447 p_dim_sets = pe->get_num_sets ();
00448 p_dim_samples = pe->get_chunk_size ();
00449 ProbeEventsList = g_slist_prepend (ProbeEventsList, e=new Event (e, se, pe, p_dim_sets, p_dim_samples));
00450
00451 gapp->progress_info_set_bar_pulse (2);
00452 }
00453 break;
00454 case 'U':
00455 {
00456 UserEntry *ue = (UserEntry*)eel->data;
00457 ue->store_event_to_nc (ncf, ++u_index, se);
00458 }
00459 break;
00460 }
00461 eel = g_slist_next (eel);
00462 }
00463
00464 sel = g_slist_next (sel);
00465 }
00466
00467 if (!ProbeEventsList)
00468 return 0;
00469
00470 gdouble sc;
00471 p_index=0;
00472 GSList *pel=ProbeEventsList;
00473 Event *prev = e = (Event*) pel->data;
00474 e->pe->write_nc_variable_set (ncf, ++p_index, evdata_var, evcoords_var, e->same_count);
00475 sc = (gdouble)e->same_count;
00476
00477 while (pel){
00478 e = (Event*) pel->data;
00479
00480
00481 if (! e->same(*prev)){
00482 e->pe->write_nc_variable_set (ncf, ++p_index, evdata_var, evcoords_var, e->same_count);
00483 sc = (gdouble)e->same_count;
00484 }
00485
00486
00487 e->pe->write_nc_data (evdata_var, evcoords_var, e->se, e->same_count-1);
00488
00489 prev=e;
00490 pel = g_slist_next (pel);
00491 gapp->progress_info_set_bar_fraction (1.-(gdouble)e->same_count/sc, 2);
00492 }
00493
00494
00495
00496 return 0;
00497 }
00498
00499 int Mem2d::LoadScanEvents (NcFile *ncf){
00500 ScanEvent *se=NULL;
00501 ProbeEntry *pe=NULL;
00502 NcVar* evdata_var=NULL;
00503 NcVar* evcoords_var=NULL;
00504
00505
00506 for (int p_index=1; ; ++p_index){
00507 int limit=1;
00508 pe = NULL;
00509 se = NULL;
00510 for (int count=0; count<limit; ++count){
00511 pe = new ProbeEntry ("Probe", ncf, p_index, evdata_var, evcoords_var, se, count, limit, pe);
00512 if (se){
00513 se->add_event (pe);
00514 AttachScanEvent (se);
00515 gapp->progress_info_set_bar_fraction ((double)count/(gdouble)limit, 2);
00516 } else break;
00517 }
00518 if (!se)
00519 break;
00520 }
00521
00522 if (pe)
00523 delete pe;
00524
00525
00526
00527 UserEntry *ue=NULL;
00528 for (int u_index=1; ; ++u_index){
00529 ue = NULL;
00530 se = NULL;
00531 ue = new UserEntry ("User", ncf, u_index, se);
00532 if (se){
00533 se->add_event (ue);
00534 AttachScanEvent (se);
00535 } else
00536 break;
00537 }
00538 if (ue)
00539 delete ue;
00540
00541 return 0;
00542 }
00543
00544 void Mem2d::AttachScanEvent (ScanEvent *se){
00545 scan_event_list = g_slist_prepend(scan_event_list, se);
00546 }
00547
00548
00549
00550 void Mem2d::Mnew(int Nx, int Ny, int Nv, ZD_TYPE type){
00551 Init();
00552 zdtyp=type;
00553 XSM_DEBUG (DBG_L6, "Mem2d::Mnew");
00554 if(data) Mdelete();
00555 switch(zdtyp){
00556 case ZD_BYTE: data = new TZData<unsigned char>(Nx, Ny, Nv); break;
00557 case ZD_SHORT: data = new TZData<SHT>(Nx, Ny, Nv); break;
00558 case ZD_LONG: data = new TZData<long>(Nx, Ny, Nv); break;
00559
00560
00561 case ZD_FLOAT: data = new TZData<float>(Nx, Ny, Nv); break;
00562 case ZD_DOUBLE: data = new TZData<double>(Nx, Ny, Nv); break;
00563 case ZD_COMPLEX: data = new TZData<double>(Nx, Ny, 3); break;
00564 case ZD_RGBA: data = new TZData<short>(Nx, Ny, 4); break;
00565 case ZD_EVENT: data = new TZData<unsigned char>(Nx, Ny, Nv); break;
00566 default:
00567 XSM_DEBUG_ERROR (DBG_L1, "Wrong type, using SHT" );
00568 data = new TZData<SHT>(Nx, Ny, Nv);
00569 break;
00570 }
00571 }
00572
00573
00574 void Mem2d::Mdelete(){
00575 XSM_DEBUG (DBG_L6, "Mem2d::Mdelete: delete data");
00576 delete data;
00577 data=NULL;
00578
00579 if (Zbins) delete [] Zbins;
00580 Zbins = NULL;
00581 if (Zibins) delete [] Zibins;
00582 Zibins = NULL;
00583 }
00584
00585 void Mem2d::Init(){
00586 SetDataPktMode(SCAN_V_DIRECT);
00587 Mod = MEM_SET;
00588 SetDataRange(0, 64);
00589 SetDataSkl(0.1, 32.);
00590 data_valid=FALSE;
00591 Zbins = NULL;
00592 Zibins = NULL;
00593 Zbin_num = 0;
00594 }
00595
00596 const char * Mem2d::GetEname(){
00597 strcpy(MemElemName, ZD_name[GetTyp()]);
00598 sprintf(MemElemName+strlen(MemElemName),"[%d]", (int)GetEsz());
00599 return((const char*)MemElemName);
00600 }
00601
00602
00603 void Mem2d::Modus(int mod){
00604 Mod = mod;
00605 }
00606
00607
00608 int Mem2d::Resize(int Nx, int Ny, ZD_TYPE type){
00609 XSM_DEBUG (DBG_L6, "Mem2d::Resize");
00610 if((zdtyp == type || type == ZD_IDENT) && data)
00611 return(data->Resize(Nx, Ny));
00612 if(type == ZD_IDENT)
00613 return -1;
00614 else
00615 Mnew(Nx, Ny, GetNv (), type);
00616 return 0;
00617 }
00618
00619 int Mem2d::Resize(int Nx, int Ny, int Nv, ZD_TYPE type){
00620 XSM_DEBUG (DBG_L6, "Mem2d::Resize with multiple layers");
00621 if((zdtyp == type || type == ZD_IDENT) && data)
00622 return(data->Resize(Nx, Ny, Nv));
00623 if(type == ZD_IDENT)
00624 return -1;
00625 else
00626 Mnew(Nx, Ny, Nv, type);
00627 return 0;
00628 }
00629
00630
00631 inline double Mem2d::GetDataPktLineReg(int x, int y){
00632 if(!data->Li[y].valid)
00633 CalcLinRegress(y,y);
00634 return(data->Z(x,y) - data->Li[y].getY(x));
00635 }
00636
00637
00638 inline double Mem2d::GetDataPktHorizont(int x, int y){
00639 if(!data->Li[y].valid)
00640 CalcLinRegress(y,y);
00641 return(data->Z(x,y) - data->Li[y].getB());
00642 }
00643
00644
00645
00646 #define KLEN 16
00647 double Mem2d::GetDataPktDiff(int x, int y){
00648 const double Lproz = 0.92;
00649 const double Rproz = 0.08;
00650 double v[KLEN];
00651 v[0] = GetDataPkt (x,y);
00652 for (int i=1; i<KLEN; ++i)
00653 if (x+i < GetNx())
00654 v[i] = Lproz * v[i-1] + Rproz * GetDataPkt (x+i,y);
00655 else break;
00656 for (int i=KLEN-2; i>=0; --i)
00657 if (x+i+1 < GetNx())
00658 v[i] = Lproz * v[i+1] + Rproz * v[i];
00659
00660 return GetDataPkt (x,y) - v[0];
00661 }
00662
00663 inline double Mem2d_LogLimit(double x){ return (x>0. ? x:0.); }
00664
00665 double Mem2d::GetDataPktLog(int x, int y){
00666 return logFac * log (1. + Mem2d_LogLimit (((GetDataPkt (x,y)) - Zmin)));
00667 }
00668
00669 double Mem2d::GetDataPktPlaneSub(int x, int y){
00670 return GetDataPkt (x,y) + plane_z0 + plane_mx*x + plane_my*y;
00671 }
00672
00673 double Mem2d::GetDataPktInterpol(double x, double y){
00674 double f1,f2,f3,f4;
00675 int x1,x2,y1,y2;
00676 x1 = (int)x;
00677 x2 = x1 + 1;
00678 y1 = (int)y;
00679 y2 = y1 + 1;
00680 f1 = ((double)x2-x) *(y-(double)y1);
00681 f2 = (x - (double)x1)*(y-(double)y1);
00682 f3 = (x - (double)x1)*((double)y2-y);
00683 f4 = ((double)x2-x) *((double)y2-y);
00684
00685 y2 = MIN (y2, GetNx()-1);
00686 x2 = MIN (x2, GetNx()-1);
00687 if(x1 < 0 || x1 >= GetNx() || y1 < 0 || y1 >= GetNy() ||
00688 x2 < 0 || x2 >= GetNx() || y2 < 0 || y2 >= GetNy())
00689 return 0.;
00690 else
00691 return( f4 * (double)GetDataPkt(x1,y1) + f3 * (double)GetDataPkt(x2,y1)
00692 + f1 * (double)GetDataPkt(x1,y2) + f2 * (double)GetDataPkt(x2,y2));
00693 }
00694
00695
00696 #define MEM2D_LINEAR_SCALE(Z) ((Z)*Zcontrast+Zbright)
00697 #define MEM2D_PERIODIC_SCALE(Z) (fmod((Z)*Zcontrast+Zbright, (double)ZVrange))
00698 #define MEM2D_RANGE_CHECK(X) ((val=(X)) > (double)ZVmax ? ZVmax : val < (double)ZVmin ? ZVmin : (ZVIEW_TYPE)val)
00699
00700 ZVIEW_TYPE Mem2d::ZQuick(int &x, int &y){
00701 double val;
00702 return MEM2D_RANGE_CHECK (MEM2D_LINEAR_SCALE (GetDataPktLineReg (x,y)));
00703 }
00704
00705 ZVIEW_TYPE Mem2d::ZDirect(int &x, int &y){
00706 double val;
00707 return MEM2D_RANGE_CHECK (MEM2D_LINEAR_SCALE (GetDataPkt (x,y)));
00708 }
00709
00710 ZVIEW_TYPE Mem2d::ZPlaneSub(int &x, int &y){
00711 double val;
00712 return MEM2D_RANGE_CHECK (MEM2D_LINEAR_SCALE (GetDataPktPlaneSub (x,y)));
00713 }
00714
00715
00716 ZVIEW_TYPE Mem2d::ZLog(int &x, int &y){
00717 double val;
00718 return MEM2D_RANGE_CHECK (MEM2D_LINEAR_SCALE (GetDataPktLog (x,y)));
00719 }
00720
00721 ZVIEW_TYPE Mem2d::ZPeriodic(int &x, int &y){
00722 double val;
00723 return MEM2D_RANGE_CHECK (MEM2D_PERIODIC_SCALE (GetDataPkt (x,y)));
00724 }
00725
00726 ZVIEW_TYPE Mem2d::ZDifferential(int &x, int &y){
00727 double val;
00728 return MEM2D_RANGE_CHECK (MEM2D_LINEAR_SCALE (GetDataPktDiff (x,y)));
00729 }
00730
00731 ZVIEW_TYPE Mem2d::ZHorizontal(int &x, int &y){
00732 double val;
00733 return MEM2D_RANGE_CHECK (MEM2D_LINEAR_SCALE (GetDataPktHorizont (x,y)));
00734 }
00735
00736 ZVIEW_TYPE Mem2d::GetDataVMode(int x, int y){
00737 return((ZVIEW_TYPE)(this->*ZVFkt)(x,y));
00738 }
00739
00740 double Mem2d::GetDataMode(int x, int y){
00741 return((double)(this->*ZFkt)(x,y));
00742 }
00743
00744 void Mem2d::SetDataRange(ZVIEW_TYPE Min, ZVIEW_TYPE Max){
00745 ZVmin=Min;
00746 ZVmax=Max;
00747 ZVrange=ZVmax-ZVmin;
00748 logFac = (double)Zrange/log(1.+fabs(Zrange));
00749 XSM_DEBUG(DBG_L5, "Mem2d::SetDataRange ZMin=" << ZVmin << " ZVMax=" << ZVmax << " => ZVRange=" << ZVrange);
00750 }
00751
00752 void Mem2d::SetDataVRangeZ(double VRangeZ, double VOffsetZ, double dz){
00753 double hi, lo, avg, zvr;
00754 GetZHiLo (&hi, &lo);
00755 avg = lo + (hi-lo)/2. + VOffsetZ/dz;
00756 zvr = VRangeZ/dz/2.;
00757 SetHiLo (avg+zvr, avg-zvr);
00758 }
00759
00760 void Mem2d::SetDataPktMode(int mode){
00761 switch (mode){
00762 case SCAN_V_QUICK: SetDataFkt (&Mem2d::ZQuick, &Mem2d::GetDataPktLineReg); break;
00763 case SCAN_V_DIRECT: SetDataFkt (&Mem2d::ZDirect, &Mem2d::GetDataPkt); break;
00764 case SCAN_V_PLANESUB: SetDataFkt (&Mem2d::ZPlaneSub, &Mem2d::GetDataPktPlaneSub); break;
00765 case SCAN_V_LOG: SetDataFkt (&Mem2d::ZLog, &Mem2d::GetDataPktLog); break;
00766 case SCAN_V_PERIODIC: SetDataFkt (&Mem2d::ZPeriodic, &Mem2d::GetDataPkt); break;
00767 case SCAN_V_DIFFERENTIAL: SetDataFkt (&Mem2d::ZDifferential, &Mem2d::GetDataPktDiff); break;
00768 case SCAN_V_HORIZONTAL: SetDataFkt (&Mem2d::ZHorizontal, &Mem2d::GetDataPktHorizont); break;
00769 }
00770 XSM_DEBUG (DBG_L5, "Mem2d::SetDataPktMode to " << mode);
00771 }
00772
00773 void Mem2d::SetDataFkt(ZVIEW_TYPE (Mem2d::*newZVFkt)(int &, int &),
00774 double (Mem2d::*newZFkt)(int, int)){
00775 ZFkt = newZFkt;
00776 ZVFkt = newZVFkt;
00777 }
00778
00779 void Mem2d::GetZHiLo(double *hi, double *lo){
00780 *hi=Zmax; *lo=Zmin;
00781 }
00782
00783 double Mem2d::GetZRange(){
00784 return Zrange;
00785 }
00786
00787 void Mem2d::SetHiLo(double hi, double lo){
00788 Zmin=lo;
00789 Zmax=hi;
00790 Zrange=Zmax-Zmin;
00791 }
00792
00793 void Mem2d::AutoDataSkl(double *contrast, double *bright){
00794
00795 if(fabs(Zrange) > 0.){
00796 Zcontrast = (double)ZVrange/Zrange;
00797 Zbright = -(double)ZVrange*Zmin/Zrange;
00798 if (contrast && bright){
00799 *contrast=Zcontrast; *bright=Zbright;
00800 }
00801 XSM_DEBUG (DBG_L6, "AutoSkl: Zmax=" << Zmax << " Zmin=" << Zmin << " ZRange=" << Zrange << " Contrast=" << Zcontrast << " Bright=" << Zbright);
00802 }
00803 }
00804
00805 void Mem2d::CalcLinRegress(int yfirst, int ylast)
00806
00807
00808
00809 {
00810 int y;
00811 double sumi, sumy, sumiy, sumyq, sumiq;
00812 int i, istep;
00813
00814 int Xmax, Xmin;
00815 double a, b, n, imean, ymean;
00816
00817 for(y=yfirst; y<=ylast; y++){
00818 sumi = sumiq = sumy = sumyq = sumiy = 0.0;
00819 Xmax = data->GetNx()-data->GetNx()/10;
00820 Xmin = data->GetNx()/10;
00821 istep = (long) MAX(1,(Xmax-Xmin)/30);
00822 for (i=Xmin; i < Xmax; i+=istep) {
00823 sumi += (double)i;
00824 sumiq += (double)i*(double)i;
00825 a = data->Z(i,y);
00826 sumy += a;
00827 sumyq += a * a;
00828 sumiy += a * (double)i;
00829 }
00830 n = (double)((Xmax-Xmin)/istep);
00831 if ((Xmax-Xmin) % istep != 0) n += 1.0;
00832 imean = sumi / n;
00833 ymean = sumy / n;
00834 a = (sumiy- n * imean * ymean ) / (sumiq - n * imean * imean);
00835 b = ymean - a * imean;
00836 data->Li[y].set(a,b);
00837 }
00838 }
00839
00840 int Mem2d::GetDataLineFrom(Point2D *start, Point2D *end, Mem2d *Mob, SCAN_DATA *sdata, GETLINEORGMODE orgmode){
00841 int i0=0;
00842 int i, LineLen;
00843 double dx, dy;
00844 double x,y;
00845
00846 XSM_DEBUG (DBG_L6, "--- Clipping ---");
00847 XSM_DEBUG (DBG_L6, "Line ("
00848 << start->x << ", " << start->y << ")-(" << end->x << ", " << end->y << ")"
00849 " ClipBox: (0,0)-(" << (Mob->GetNx()-1) << ", " << (Mob->GetNx()-1) << ")");
00850
00851
00852 cohen_sutherland_line_clip_i (&(start->x), &(start->y), &(end->x), &(end->y),
00853 0, Mob->GetNx()-1,
00854 0, Mob->GetNy()-1);
00855
00856 XSM_DEBUG (DBG_L6, "Result: (" << start->x << ", " << start->y << ")-(" << end->x << ", " << end->y << ")");
00857
00858
00859 dx = (double)(end->x - start->x);
00860 dy = (double)(end->y - start->y);
00861
00862
00863 LineLen = (int)sqrt(dx*dx+dy*dy);
00864
00865 Resize(LineLen,1);
00866
00867
00868 dx /= (double)(LineLen-1);
00869 dy /= (double)(LineLen-1);
00870
00871 for (i=0; i < LineLen; ++i)
00872 PutDataPkt(Mob->GetDataPktInterpol((double)start->x + (double)i * dx,
00873 (double)start->y + (double)i * dy),
00874 i, i0);
00875
00876 double xy[4];
00877 xy[0] = Mob->data->GetXLookup((int)start->x);
00878 xy[1] = Mob->data->GetYLookup((int)start->y);
00879 xy[2] = Mob->data->GetXLookup((int)end->x);
00880 xy[3] = Mob->data->GetYLookup((int)end->y);
00881 x = data->GetXLookup((int)(start->x+dx/2));
00882 y = data->GetYLookup((int)(start->y+dy/2));
00883 sdata->s.x0 = x;
00884 sdata->s.y0 = y;
00885 sdata->s.dx *= dx;
00886 sdata->s.dy *= dy;
00887
00888 sdata->s.dl = sqrt((xy[2]-xy[0])*(xy[2]-xy[0]) + (xy[3]-xy[1])*(xy[3]-xy[1]));
00889 sdata->s.rx = sdata->s.dl;
00890 sdata->s.ry = 0.;
00891 sdata->s.nx = LineLen;
00892 sdata->s.ny = 1;
00893
00894 switch(orgmode){
00895 case GLORG_CENTER: data->MkXLookup(-sdata->s.rx/2., sdata->s.rx/2); break;
00896 case GLORG_PRJX: data->MkXLookup(xy[0], xy[2]); break;
00897 case GLORG_PRJY: data->MkXLookup(xy[1], xy[3]); break;
00898 case GLORG_ZERO:
00899 default: data->MkXLookup(0., sdata->s.rx); break;
00900 }
00901
00902 if(fabs(dx) > 1e-5)
00903 sdata->s.alpha = 180./M_PI*atan(-dy/dx);
00904 else
00905 sdata->s.alpha = 90.*(-dy>0.?1.:-1.);
00906 if(dx < 0.)
00907 sdata->s.alpha = 180.+sdata->s.alpha;
00908 else
00909 if(dy>0.)
00910 sdata->s.alpha = 360.+sdata->s.alpha;
00911 return 0;
00912 }
00913
00914 int Mem2d::GetLayerDataLineFrom(Point2D *start, Mem2d *Mob, SCAN_DATA *sdata){
00915 Resize (Mob->GetNv(), 1);
00916
00917 for (int i=0; i < Mob->GetNv(); ++i)
00918 PutDataPkt (Mob->GetDataPkt((int)start->x,
00919 (int)start->y,
00920 i),
00921 i, 0);
00922
00923 sdata->s.nx = Mob->GetNv();
00924 sdata->s.ny = 1;
00925
00926 sdata->s.dl = Mob->data->GetVLookup (Mob->GetNv()-1) - Mob->data->GetVLookup (0);
00927
00928 sdata->s.rx = sdata->s.dl;
00929 sdata->s.ry = 0.;
00930
00931 sdata->s.x0 = Mob->data->GetXLookup (0);
00932 sdata->s.y0 = 0.;
00933 sdata->s.dx = sdata->s.rx/(sdata->s.nx-1);
00934 sdata->s.dy = 0.;
00935
00936 data->MkXLookup(Mob->data->GetVLookup (0), Mob->data->GetVLookup (Mob->GetNv()-1));
00937
00938 sdata->s.alpha = 0.;
00939
00940 if (sdata->Vunit)
00941 sdata->SetXUnit (sdata->Vunit);
00942
00943 return 0;
00944 }
00945
00946
00947 int Mem2d::GetArcDataLineFrom(Point2D *center, Point2D *radius, Mem2d *Mob, SCAN_DATA *sdata){
00948 double dx, dy, r, circumference, dphi;
00949 int n;
00950
00951
00952 dx = (double)(center->x - radius->x);
00953 dy = (double)(center->y - radius->y);
00954 r = sqrt(dx*dx+dy*dy);
00955 circumference = 2.*r*M_PI;
00956 n = (int)circumference;
00957 dphi = 2.*M_PI/(double)n;
00958
00959 Resize (n, 1);
00960
00961 for (int i=0; i < n; ++i)
00962 PutDataPkt (Mob->GetDataPktInterpol (center->x + r*cos((double)i*dphi),
00963 center->y + r*sin((double)i*dphi)
00964 ),
00965 i, 0);
00966
00967 sdata->s.nx = n;
00968 sdata->s.ny = 1;
00969
00970 sdata->s.dl = 360.;
00971
00972 sdata->s.rx = sdata->s.dl;
00973 sdata->s.ry = r*sdata->s.dx;
00974
00975 sdata->s.x0 = center->x;
00976 sdata->s.y0 = center->y;
00977 sdata->s.dy = r*sdata->s.dx;
00978 sdata->s.dx = 180.*dphi/M_PI;
00979
00980 data->MkXLookup(0., 360.);
00981
00982 sdata->s.alpha = 0.;
00983
00984 UnitObj *ArcUnit = new UnitObj (UTF8_DEGREE, "\260");
00985 sdata->SetXUnit (ArcUnit);
00986 delete ArcUnit;
00987
00988 return 0;
00989 }
00990
00991 void Mem2d::DataRead(std::ifstream &f, int ydir){
00992 int y;
00993 if (ydir>0){
00994 ydir = 1;
00995 y = 0;
00996 } else {
00997 ydir = -1;
00998 y = data->GetNy()-1;
00999 }
01000 for(int x = 0, l = 0; l<data->GetNy(); l++, y+=ydir){
01001 if(f.good()){
01002 if (WORDS_BIGENDIAN){
01003
01004 int bpz = data->Zsize();
01005 char *tmp = new char[bpz*(data->GetNx()+1)];
01006 f.read(tmp+bpz, bpz*data->GetNx());
01007 for(int i=0; i<bpz*data->GetNx(); i+=bpz)
01008 for(int j=0; j<bpz; ++j)
01009 tmp[i+j] = tmp[i+2*bpz-j-1];
01010 memcpy(data->GetPtr(x,y), tmp, data->Zsize()*data->GetNx());
01011 delete[] tmp;
01012 }else{
01013 f.read((char*)data->GetPtr(x,y), data->Zsize()*data->GetNx());
01014 }
01015 }
01016 else
01017 XSM_DEBUG (DBG_L6, "File ends unexpectetly !");
01018 data->Li[y].invalidate();
01019 }
01020 data_valid=1;
01021 }
01022
01023
01024
01025
01026 void Mem2d::DataD2DRead(std::ifstream &f, double gate){
01027 int i,j;
01028 double fac, val;
01029
01030
01031 short *buf, *ptr;
01032 ptr = buf = g_new( short, (data->GetNx()+1)*(data->GetNy()+1) );
01033 f.read((char*)buf, sizeof(short)*(data->GetNx()+1)*data->GetNy());
01034
01035 fac = (gate*1e-3)/2500.;
01036 XSM_DEBUG (DBG_L6, "D2DRead factor:" << fac);
01037 for(i=0; i<data->GetNy(); i++){
01038 for(j=0; j<data->GetNx(); j++,buf++){
01039
01040 if (WORDS_BIGENDIAN)
01041 val = (double) (GINT16_FROM_LE (*buf)) + 32768.;
01042 else
01043 val = (double) (*buf) + 32768.;
01044
01045 data->Z(val*val*fac,j,i);
01046 }
01047 buf++;
01048 data->Li[i].invalidate();
01049 }
01050 g_free( ptr );
01051 data_valid=1;
01052 }
01053
01054 void Mem2d::DataWrite(std::ofstream &f){
01055 for(int x=0, y=0; y<data->GetNy(); y++)
01056 f.write((const char*)data->GetPtr(x,y), data->Zsize()*data->GetNx());
01057 }
01058
01059 void Mem2d::AutoHistogrammEvalMode (Point2D *p1, Point2D *p2, int delta, double epsilon){
01060 double dz_norm;
01061
01062 int med_bin, lo_bin, hi_bin;
01063 int nx0, nx, ny0, ny;
01064
01065
01066 HiLoMod (p1, p2, delta);
01067
01068
01069 if(p1 && p2){
01070 ny0 = MIN (MAX (MIN (p1->y, p2->y), 0), (data->GetNy()-1));
01071 ny = MIN (MAX (MAX ((p1->y-delta), (p2->y-delta)), 0), (data->GetNy()-1));
01072 nx0 = MIN (MAX (MIN (p1->x, p2->x), 0), (data->GetNx()-1));
01073 nx = MIN (MAX (MAX ((p1->x-delta), (p2->x-delta)), 0), (data->GetNx()-1));
01074 }
01075 else{
01076 nx0=ny0=0;
01077 nx=data->GetNy()-1;
01078 ny=data->GetNx()-1;
01079 }
01080
01081 if (nx0 > (nx-delta) || ny0 > (ny-delta))
01082 return;
01083
01084
01085 Zbin_num = MAX((int)(Zrange/3), 10);
01086 Zbin_width = Zrange / Zbin_num;
01087 dz_norm = 1./Zbin_width;
01088
01089 if (Zbins) delete [] Zbins;
01090 Zbins = new double[Zbin_num];
01091 if (Zibins) delete [] Zibins;
01092 Zibins = new double[Zbin_num];
01093
01094 for(int i = 0; i < Zbin_num; i++)
01095 Zbins[i] = 0.;
01096
01097 for(int col = nx0; col < nx; col+=delta)
01098 for(int line = ny0; line < ny; line+=delta){
01099 double f = (GetDataMode (col,line) - Zmin) / Zbin_width;
01100 int bin = (int)f;
01101
01102 if (bin < Zbin_num){
01103 double f1 = (bin+1) * Zbin_width;
01104 if ((f+dz_norm) > f1){
01105 Zbins [bin] += f-bin;
01106 ++bin;
01107 if (bin < Zbin_num)
01108 Zbins [bin] += bin-f;
01109 }
01110 else
01111 Zbins [bin] += 1.;
01112 }
01113 }
01114
01115
01116 Zibins [0] = Zbins [0];
01117 for(int i = 1; i < Zbin_num; i++)
01118 Zibins[i] = Zibins[i-1] + Zbins[i];
01119
01120 Zilow = Zibins [2];
01121 Zihigh = Zibins [Zbin_num-3];
01122 Zirange = Zihigh - Zilow;
01123 Zimed = Zilow + Zirange/2.;
01124
01125
01126 for (med_bin = 1; Zibins [med_bin] < Zimed; ++med_bin);
01127
01128 #define DIFF(I) (Zibins [I+1] - Zibins [I-1])
01129
01130
01131 for (lo_bin = med_bin; lo_bin > 2 && Zibins [lo_bin] > (Zilow+Zirange*epsilon); --lo_bin);
01132 for (hi_bin = med_bin; hi_bin < (Zbin_num-3) && Zibins [hi_bin] < (Zihigh-Zirange*epsilon); ++hi_bin);
01133
01134
01135
01136
01137
01138 Zmax = hi_bin*Zbin_width + Zmin;
01139 Zmin = lo_bin*Zbin_width + Zmin;
01140 Zrange=Zmax-Zmin;
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153 }
01154
01155
01156 void Mem2d::HiLoMod(Point2D *p1, Point2D *p2, int Delta){
01157 int i,j,i1,i2, j1,j2;
01158 double p, hi, lo;
01159
01160 if(p1 && p2){
01161 i1 = MIN (MAX (MIN (p1->y, p2->y), 0), (data->GetNy()-1));
01162 i2 = MIN (MAX (MAX ((p1->y-Delta), (p2->y-Delta)), 0), (data->GetNy()-1));
01163 j1 = MIN (MAX (MIN (p1->x, p2->x), 0), (data->GetNx()-1));
01164 j2 = MIN (MAX (MAX ((p1->x-Delta), (p2->x-Delta)), 0), (data->GetNx()-1));
01165 }
01166 else{
01167 i1=j1=0;
01168 i2=data->GetNy()-1;
01169 j2=data->GetNx()-1;
01170 }
01171
01172
01173
01174
01175 hi=lo=GetDataMode(j1, i1);
01176
01177 if (i1 > (i2-Delta) || j1 > (j2-Delta)){
01178
01179
01180 hi *= 1.2;
01181 lo *= 0.8;
01182 }
01183 else{
01184 for(i=i1; i<i2; i+=Delta)
01185 for(j=j1; j<j2; j+=Delta){
01186 p=GetDataMode(j, i);
01187 if(hi < p) hi = p;
01188 else
01189 if(lo > p) lo = p;
01190 }
01191 }
01192 Zmin=lo;
01193 Zmax=hi;
01194 Zrange=Zmax-Zmin;
01195 XSM_DEBUG (DBG_L6, "Mem2d::HiLoMod Zmin=" << Zmin << " Zmax=" << Zmax);
01196 }
01197
01198 void Mem2d::HiLo(double *hi, double *lo, int LinReg, Point2D *p1, Point2D *p2, int Delta){
01199 int i,j,i1,i2, j1,j2;
01200 double Hi, Lo;
01201 i1 = p1 ? p1->y : 0;
01202 i2 = p2 ? p2->y : data->GetNy();
01203 j1 = p1 ? p1->x : 0;
01204 j2 = p2 ? p2->x : data->GetNx();
01205 if(LinReg){
01206 double p;
01207 Hi=Lo=GetDataPktLineReg(j1, i1);
01208 for(i=i1; i<i2; i+=Delta)
01209 for(j=j1; j<j2; j+=Delta){
01210 p=GetDataPktLineReg(j, i);
01211 if(Hi < p) Hi = p;
01212 else
01213 if(Lo > p) Lo = p;
01214 }
01215 }
01216 else{
01217 data->SetPtr(j1,i1);
01218 Hi = Lo = data->GetThis();
01219 for(i=i1; i<i2; i++)
01220 for(j=j1, data->SetPtr(j1,i); j<j2; j++){
01221 if(Hi < data->GetThis()) Hi = data->GetNext();
01222 else
01223 if(Lo > data->GetThis()) Lo = data->GetNext();
01224 else
01225 data->GetNext();
01226 }
01227 }
01228
01229 *hi = (double)Hi;
01230 *lo = (double)Lo;
01231 }
01232
01233 int Mem2d::DataValid(){ return data_valid; }
01234
01235 int Mem2d::SetDataValid(){ return data_valid=1; }
01236
01237
01238
01239 MemDigiFilter::MemDigiFilter(double Xms, double Xns, int M, int N):Mem2d(2*N+1, 2*M+1, ZD_DOUBLE){
01240 xms=Xms, xns=Xns;
01241 n=N, m=M;
01242 XSM_DEBUG (DBG_L4, "MemDigiFilter: " << n << " " << m << " " << xms << " " << xns);
01243 }
01244
01245 gboolean MemDigiFilter::Convolve(Mem2d *Src, Mem2d *Dest){
01246 int i0=0;
01247 int i, j;
01248 int ii, jj;
01249 double scalefac, scalefaca;
01250 int mm, nn;
01251 int ms=m, ns=n;
01252 gboolean again = FALSE;
01253
01254 XSM_DEBUG (DBG_L2, "MemDigiFilter::Convolve: create kernel");
01255
01256 do{
01257 int ring_m, ring_n;
01258
01259 XSM_DEBUG (DBG_L6, "R " << (2*ms+1) << " " << (2*ns+1));
01260 n=ns; m=ms;
01261 Resize(2*ns+1, 2*ms+1);
01262 XSM_DEBUG (DBG_L6, "MemDigiFilter::Resize done.");
01263 CalcKernel();
01264 XSM_DEBUG (DBG_L6, "MemDigiFilter::CalcKernel done.");
01265 for(ring_m=ring_n=-1, scalefaca=scalefac=0., i=0; i<2*ms+1; i++)
01266 for(j=0; j<2*ns+1; j++){
01267 int tmp_ring_m = abs(i-m);
01268 int tmp_ring_n = abs(j-n);
01269 scalefaca += fabs(data->Z(j,i));
01270 scalefac += data->Z(j,i);
01271 if(data->Z(j,i)!=0. && ring_m<tmp_ring_m)
01272 ring_m = tmp_ring_m;
01273 if(data->Z(j,i)!=0. && ring_n<tmp_ring_n)
01274 ring_n = tmp_ring_n;
01275 }
01276 scalefac = scalefac!=0. ? fabs(scalefac) : (scalefaca+1.)/2.;
01277
01278 again = FALSE;
01279 if(ring_m < m) {again = TRUE; ms = ring_m;}
01280 if(ring_n < n) {again = TRUE; ns = ring_n;}
01281 }while(again);
01282 XSM_DEBUG (DBG_L6, "MemDigiFilter::Scaleing done.");
01283
01284 if(Src->data->GetNx()<1 && Src->data->GetNy()<1)
01285 return FALSE;
01286
01287 XSM_DEBUG (DBG_L6, "Prep Larger Src: " << (Src->data->GetNx() + 2*ns) << " " << (Src->data->GetNy() + 2*ms));
01288
01289 Mem2d x((nn=Src->data->GetNx()) + 2*ns, (mm=Src->data->GetNy()) + 2*ms, Src->GetTyp());
01290
01291
01292
01293 XSM_DEBUG (DBG_L6, "Mem2dDigi: " << ns << " " << ms << " " << nn << " " << mm);
01294 x.data->CopyFrom(Src->data, 0,0, ns,ms ,nn,mm);
01295
01296
01297
01298 for(i=0;i<mm;i++)
01299 for(j=0;j<ns;j++){
01300 x.data->Z(Src->data->Z(i0,i), j, i+ms);
01301 x.data->Z(Src->data->Z(nn-1,i),j+ns+nn,i+ms);
01302 }
01303
01304
01305 for(i=0;i<ms;i++){
01306 x.data->CopyFrom(Src->data, 0,0, ns,i ,nn);
01307 x.data->CopyFrom(Src->data, 0,mm-1, ns,i+mm+ms ,nn);
01308 }
01309
01310
01311 for(i=0;i<ms;i++)
01312 for(j=0;j<ns;j++){
01313 x.data->Z(Src->data->Z(i0,i0), j,i);
01314 x.data->Z(Src->data->Z(nn-1,i0), nn+ns+j,i);
01315 x.data->Z(Src->data->Z(i0,mm-1), j,mm+ms+i);
01316 x.data->Z(Src->data->Z(nn-1,mm-1), nn+ns+j,mm+ms+i);
01317 }
01318
01319 #ifdef SAVECONVOLKERN
01320
01321 std::ofstream fk;
01322 fk.open("/tmp/convolkern.dbl", ios::out|ios::bin);
01323 struct { short nx,ny; } fkh;
01324 fkh.nx=data->GetNx();
01325 fkh.ny=data->GetNy();
01326 fk.write((void*)&fkh, sizeof(fkh));
01327 DataWrite(fk);
01328 fk.close();
01329 #endif
01330 #ifdef SAVECONVOLSRC
01331
01332 std::ofstream fk;
01333 switch(x.GetTyp()){
01334 case ZD_BYTE: fk.open("/tmp/convolsrc.byt", std::ios::out); break;
01335 case ZD_SHORT: fk.open("/tmp/convolsrc.sht", std::ios::out); break;
01336 case ZD_LONG: fk.open("/tmp/convolsrc.byt", std::ios::out); break;
01337 case ZD_DOUBLE: fk.open("/tmp/convolsrc.dbl", std::ios::out); break;
01338 default: break;
01339 }
01340 if(fk.good()){
01341 struct { short nx,ny; } fkh;
01342 fkh.nx=x.data->GetNx();
01343 fkh.ny=x.data->GetNy();
01344 fk.write((const char*)&fkh, sizeof(fkh));
01345 x.DataWrite(fk);
01346 fk.close();
01347 }
01348 #endif
01349
01350
01351 if (Dest->GetNv () == 1)
01352 Dest->Resize(nn,mm);
01353
01354
01355 gint pcent=0;
01356 for(ii=0; ii<mm; ++ii){
01357 if(pcent < 100*ii/mm ){
01358 SET_PROGRESS((gfloat)ii/(gfloat)mm);
01359 }
01360
01361 for(jj=0; jj<nn; ++jj){
01362 const int tms = 2*ms;
01363 const int tns = 2*ns;
01364 const int sf2 = (int)scalefac/2;
01365 double sum;
01366
01367 for(sum=0., i=0; i<=tms; ++i){
01368 data->SetPtr(0,i);
01369 x.data->SetPtr(jj,i+ii);
01370 for(j=0; j<tns; j++)
01371 sum += x.data->GetNext() * data->GetNext();
01372 }
01373 Dest->data->Z(sum>0. ? ((sum+sf2)/scalefac) :
01374 -((-sum+sf2)/scalefac),
01375 jj,ii);
01376 }
01377 }
01378 SET_PROGRESS(0);
01379 XSM_DEBUG (DBG_L6, "done.");
01380 return TRUE;
01381 }