mem2d.C

Go to the documentation of this file.
00001 /* Gxsm - Gnome X Scanning Microscopy
00002  * universal STM/AFM/SARLS/SPALEED/... controlling and
00003  * data analysis software
00004  * 
00005  * Copyright (C) 1999,2000,2001,2002,2003 Percy Zahl
00006  *
00007  * Authors: Percy Zahl <zahl@users.sf.net>
00008  * additional features: Andreas Klust <klust@users.sf.net>
00009  * WWW Home: http://gxsm.sf.net
00010  *
00011  * This program is free software; you can redistribute it and/or modify
00012  * it under the terms of the GNU General Public License as published by
00013  * the Free Software Foundation; either version 2 of the License, or
00014  * (at your option) any later version.
00015  *
00016  * This program is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU General Public License
00022  * along with this program; if not, write to the Free Software
00023  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
00024  */
00025 
00026 /* -*- Mode: C++; indent-tabs-mode: nil; c-basic-offset: 8 c-style: "K&R" -*- */
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 // used to set progressbar of main window (P: 0..1)
00039 #define SET_PROGRESS(P) { gapp->SetProgress((gfloat)(P)); while (gtk_events_pending()) gtk_main_iteration(); }
00040 
00041 //#define SAVECONVOLKERN
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  * Low Level 2d Z Data Verwaltung:
00051  * beliebiger skalarer ZData-Typ "ZTYP"
00052  *
00053  * Objekt: ZData abstraktes Basisobject
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; // keep old 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  * Template Objekt: TZData
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 // accesses via real-lookup-table for future use
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 // set data
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 // set Z values to const in opt. area
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) // restrict to layer v?
00264                 vs=ve=v;
00265         if (!xs) xs = nx;
00266         if (!ys) ys = ny;
00267         // set one line
00268         for (int i=x0; i < (x0+xs); Zdat[y0*nv+vs][i++] = zval);
00269         // copy rest
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 // NetCDF Put / Get Methods "overloaded"
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  * 2d Memory Verwaltung:
00310  * class Mem2d
00311  *
00312  */
00313 
00314 // make new with size and type
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 // make copy or empty
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 // make with size from same type
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) // all
00384                 g_slist_foreach(scan_event_list, (GFunc)report_obj_func, gp);
00385         else{ // only within distance to xy[2]
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; // List of ScanEvents
00408         int p_index; // Event_Probe_####_... NcVar name counter, grouping of similar (same dimensions) data sets
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         // parse scan events
00437         sel = scan_event_list; // List of ScanEvents
00438         while (sel){
00439                 ScanEvent *se = (ScanEvent*)sel->data; // Event (contains different EventEntry types)
00440                 GSList *eel = se->event_list; // List of all EventEnties @ ScanEvent (one Location)
00441                 while (eel){
00442                         EventEntry *ee = (EventEntry*)eel->data; // EventEntry (contains Data)
00443                         switch (ee->description_id ()){
00444                         case 'P': // "Probe" Event
00445                           {
00446                                 ProbeEntry *pe = (ProbeEntry*)eel->data; // it's a ProbeEntry
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': // "User" Event, just store away
00455                           {     
00456                                 UserEntry *ue = (UserEntry*)eel->data; // it's a UserEntry
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                 // changed dimensions? Then write new nc-var set.
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                 // write data
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         // free ProbeEventsList and Elements!
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         // scan NC-file for ScanEvents
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         // scan NC-file for UserEvents
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 /* Speicher für 2D Daten anfordern */
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                 //  case ZD_ULONG:  data = new TZData<unsigned long>(Nx, Ny, Nv); break;
00560                 // ... since not supportted by NetCDF it has to be excluded !
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 /* Speicher für 2D Daten freigeben */
00574 void Mem2d::Mdelete(){
00575         XSM_DEBUG (DBG_L6, "Mem2d::Mdelete: delete data");
00576         delete data;
00577         data=NULL;
00578         // free statistics stuff
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 /* Speichermodus setzen */
00603 void Mem2d::Modus(int mod){
00604         Mod = mod;
00605 }
00606 
00607 /* Neue Speichergröße für 2D Daten */
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 // Get Data with LineReg !
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 // Get Data with Horizontmethode !
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 // Get Data with "koehler-type" differential filter method !
00645 // KLEN=16 (default) may be adjusted for better/faster results
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         // erst hier die Sicherheit, damit nicht f1 = f2 = f3 = f4 = 0 
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         // maybe call HiLoMod() before !
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 // OutputFeld[i] = InputFeld[i] - (a * i + b)
00808 //
00809 {
00810         int y;
00811         double sumi, sumy, sumiy, sumyq, sumiq;
00812         int i, istep;
00813         // short *lpi, *lpo ;
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; // etwas Abstand (10%) vom Rand !
00820                 Xmin = data->GetNx()/10;
00821                 istep = (long) MAX(1,(Xmax-Xmin)/30); // 30 Samples / Line !!!
00822                 for (i=Xmin; i < Xmax; i+=istep) { /* Rev 2 */
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   // Clipping
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   // Delta mit Vorzeichen !!!
00859   dx = (double)(end->x - start->x);
00860   dy = (double)(end->y - start->y);
00861   
00862   // Laenge der Linie in Pixeln
00863   LineLen = (int)sqrt(dx*dx+dy*dy);
00864   
00865   Resize(LineLen,1);
00866   
00867   // double Schrittweite im Feld 
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         // Radius
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                                 // correct endianess
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 // I'm able to import D2D !!!!
01024 // Read: 0..65536 =(^2)=> 0 .. 2^32/2500 = 1.7e6
01025 // CNT =  tau * XCPS*XCPS/2500
01026 void Mem2d::DataD2DRead(std::ifstream &f, double gate){
01027         int i,j;
01028         double fac, val;
01029         //  unsigned short *buf;
01030         //  buf = new unsigned short[(nx+1)*(ny+1)];
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         //  fac = (spa->data.GateTime*1e-3)/2500.; // spa4
01035         fac = (gate*1e-3)/2500.; // spa4
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                         //      d2cnt[i][j] = (CNT)*buf*(CNT)*buf;
01040                         if (WORDS_BIGENDIAN)
01041                                 val = (double) (GINT16_FROM_LE (*buf)) + 32768.;
01042                         else
01043                                 val = (double) (*buf) + 32768.;
01044                         //      d2cnt[i][j] = (CNT)(val*val*fac);
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 //      double m_imed;
01062         int    med_bin, lo_bin, hi_bin;
01063         int nx0, nx, ny0, ny;
01064 
01065         // find full data range in area
01066         HiLoMod (p1, p2, delta);
01067 
01068         // check area
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; // ERROR!!
01083 
01084         // calculate histogramm, auto bin size
01085         Zbin_num   = MAX((int)(Zrange/3), 10);  // +/-1 dz (3dz) in ein bin per default
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){ // partiell in bin, immer rechter Rand, da bin>0 && bin < f
01105                                         Zbins [bin] += f-bin;
01106                                         ++bin;
01107                                         if (bin < Zbin_num)
01108                                                 Zbins [bin] += bin-f; // 1.-(f-(bin-1))
01109                                 }
01110                                 else // full inside of bin
01111                                         Zbins [bin] += 1.;
01112                         }
01113                 }
01114 
01115         // integrate histogramm data and auto evaluate relevant area around medium
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 //      m_imed = DIFF(med_bin);
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 // not reliable
01135 //      for (lo_bin = med_bin; lo_bin > 2 && DIFF(lo_bin) > 0.; --lo_bin);
01136 //      for (hi_bin = med_bin; hi_bin < (Zbin_num-3) && DIFF(hi_bin) > 0.; ++hi_bin);
01137 
01138         Zmax = hi_bin*Zbin_width + Zmin;
01139         Zmin = lo_bin*Zbin_width + Zmin;
01140         Zrange=Zmax-Zmin;
01141 
01142 /*
01143         std::ofstream h; h.open ("/tmp/hh", std::ios::out);
01144         h << "#ZAutoStatistics  Zmin=" << Zmin << " Zmax=" << Zmax << "\n";
01145         h << "#ZAutoStatistics  Zilow=" << Zilow << " Zimed=" << Zimed << " Zihigh=" << Zihigh << "\n";
01146 
01147         { int i=0; h << (i*Zbin_width + Zmin) << " " << Zbins[i] << " " << Zibins[i] << " " << "0 -2" << " " << i << "\n"; }
01148         for (int i=1; i<Zbin_num-1; ++i)
01149                 h << (i*Zbin_width + Zmin) << " " << Zbins[i] << " " << Zibins[i] << " " << DIFF(i) << " " << (i>lo_bin && i<hi_bin ? -1:-2) << " " << i << "\n";
01150         { int i=Zbin_num-1; h << (i*Zbin_width + Zmin) << " " << Zbins[i] << " " << Zibins[i] << " " << "0 -2" << " " << i << "\n"; }
01151         h.close ();
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 //      std::cerr << "Mem2d::HiLoMod " << i1 << " " << i2 << " " << j1 << " " << j2 << " :: " << data->GetNx() << ":" << data->GetNy() << std::endl;
01173         
01174 
01175         hi=lo=GetDataMode(j1, i1);
01176 
01177         if (i1 > (i2-Delta) || j1 > (j2-Delta)){
01178 //              std::cerr << "Mem2d::HiLoMod ERROR" << std::endl;
01179                 // bad area, make some defaults
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   //  if(Hi < 1) Hi=1L; // Div Zero prevent....
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 // MemDigiFilter - public Mem2d
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   // create the convolution kernel
01256   do{
01257     int ring_m, ring_n;
01258     // (Re)Calc Kernel ... xms, xns, ms, ns
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   // mm * nn is now size of Src:
01289   Mem2d x((nn=Src->data->GetNx()) + 2*ns, (mm=Src->data->GetNy()) + 2*ms, Src->GetTyp());
01290  
01291   // *** WORKINGMARKER *** 11/2/1999 PZ ***
01292   // fill the central part of the x matrix with the data
01293   XSM_DEBUG (DBG_L6, "Mem2dDigi: " << ns << " " << ms << " " << nn << " " << mm);
01294   x.data->CopyFrom(Src->data, 0,0, ns,ms ,nn,mm); // !!!!!!!!!!!! tot
01295   
01296   // now fill edges and corners with copies of edge data
01297   // edge left / right
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); // data[i+ms][j]       = Src->data[i][0];
01301       x.data->Z(Src->data->Z(nn-1,i),j+ns+nn,i+ms); // data[i+ms][j+ns+nn] = Src->data[i][nn-1];
01302     }
01303   
01304   // edge top / bottom
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   // corners
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   // save Kernel to /tmp/convolkern.dbl
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   // save datasrc with added bounds to /tmp/convolsrc.xxx
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   // do convolution !
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 }

Generated on Sat Apr 1 09:04:01 2006 for GXSM by  doxygen 1.4.6