VISH  0.2
ComputeMultiStreamLines.cpp
#include <ocean/plankton/VModules.hpp>
#include <ocean/plankton/VCreator.hpp>
#include <bone/FishField.hpp>
#include <bone/GridObject.hpp>
#include <ocean/GLvish/VRenderObject.hpp>
#include <ocean/GLvish/BoundingBox.hpp>
#include <ocean/eagle/PhysicalSpace.hpp>
#include <ocean/shrimp/VEnum.hpp>
#include <field/DirectProductArray.hpp>

#include <ocean/plankton/VOperatorCache.hpp>
#include <ocean/plankton/VTimer.hpp>

#include <fish/fiber/vector/Interpolate.hpp>
#include <fish/fiber/vector/LinearIpol.hpp>
#include <fish/fiber/baseop/UniGridMapper.hpp>
#include <fish/fiber/baseop/ExpandBBox.hpp>
#include <ocean/eagle/KDTree.hpp>
#include <fish/fiber/field/Fragment.hpp>
#include <fish/fiber/grid/Chart.hpp>

#include <bone/GridObject.hpp>

#include "EmitterPoints.hpp"
#include "LinesPoints.hpp"

//#define TIMEIT

using namespace Wizt;
using namespace Fiber;
using namespace Eagle;

void ________________________________________________getMTime(double&sec, VTimer&timer)
{
#ifdef TIMEIT
        sec = timer.secs();
#endif
}

void   ________________________________________________printMTime(const char*text, const double time)
{
#ifdef TIMEIT
        printf("----------------------------------------------------------------------------------- %s:\t%.4f\n", text, time);
#endif
}

class   ComputeMultiStreamLines : public virtual VRenderObject,
                                  public virtual Fish<Grid>,
                                  public virtual Fish<Skeleton>,
                                  public virtual Fish<Field>
{ 

public:
        typedef MemArray<3, Eagle::tvector3>           VectorArray_t; 
        typedef MemArray<3, Eagle::point3>             CoordsArray_t; 
        typedef DirectProductMemArray<Eagle::point3>   ProcArray_t; 

        typedef CreativeIterator<Eagle::point3>         CoordsIterator_t;
        typedef CreativeIterator<double         >       ScalarIterator_t; 
        typedef CreativeIterator<Eagle::tvector3>       VectorIterator_t;
        typedef std::vector<point> Line_t;

        typedef pair<point, point> line_t;
        typedef pair<line_t, rgba_float_t> colorLine_t;
        typedef pair<point, rgba_float_t> colorPoint_t;         
        

        struct  PerBlockInfo : public ReferenceBase<PerBlockInfo>
        {
                point                  MidPoint;
                double                 radius;
                RefPtr<UniGridMapper>  UniMap;
                RefPtr<FragmentID>     FragID;
                MultiIndex<3>          Offset;

                PerBlockInfo()
                : ReferenceBase<PerBlockInfo>(this)
                        {}
        };      

        

        struct  MyState : State
        {
                int RefinementLevel;

                bool check_all;

                std::vector<int> candidates;

                //std::vector<MultiIndex<3> > cand_coord;
                //std::vector<int > cand_coord_cand;

                // used for debug drawing only
                std::vector<colorLine_t> lines;
                std::vector<colorPoint_t> points;

                RefPtr<Field> Coords;
                
                // Iteraror saves, wheter we have curvilinear or procedural coordinates in the fragments here.
                bool have_curvi_coords;
                
                bool have_fragmented_coords;
                bool have_fragmented_data;

                // The field that is used to store the GL cache
                RefPtr<Field> F;

                
                //  store midpoints. later this should be exchanged by the kdtree
                typedef std::vector<RefPtr<PerBlockInfo> > Blocks_t;
                Blocks_t                                Blocks;

                RefPtr<KDTree<3, int> >                 Tree;
                double                                  max_radius, last_radius, tree_fac;
                

                MyState():have_curvi_coords(0), have_fragmented_coords(0), have_fragmented_data(0), max_radius(0.0), last_radius(0.0){}
        };

        // iterator only needed to prepare info for a kd tree and insert info into
        // the kdtree
        struct  It : FieldFragmentIterator
        {
                int c;
                RefPtr<MyState> S;
                
                It(const RefPtr<MyState>_S):c(0), S(_S){}
                
                bool apply(const RefPtr<FragmentID>&f,
                           const MemCore::StrongPtr<Fiber::CreativeArrayBase>&DC)
                        
                        {
                                printf("entering iterator\n");fflush(stdout);
                                assert(S);
                                if(!S->Tree)
                                        S->Tree = new KDTree<3, int>(); 
                                assert(S->Tree); 

                                RefPtr<PerBlockInfo> VS = new PerBlockInfo(); 

                                // store Fragment id in State for random accessability of fragments in the update functino via the state object 
                                VS->FragID = f; 
                                
                                // get dimensions of fragment
                        MultiIndex<3> fi, li, lix, liy, liz; 
                        bool test; 
                                test = makeDimensions( li, DC); 
                                li = li[0]-1, li[1]-1, li[2]-1;  // <---------------------------------------------- ging nit in Optimize 
                                if(!test)
                                {
                                        printf("It::apply() Error getting fragment dimsensions"); 
                                        return false;
                                } 

                                // get Coordinate offset of Fragment (in case there's the same Positions field used for multiple data Blocks)
                         MultiIndex<3> Offset; 
                                if ( RefPtr<FragmentLocation<3> > F3 = interface_cast<FragmentLocation<3> > (f) )
                                {
                                        Offset = F3->Offset; 
                                } 
                                else
                                {
                                        Offset = 0, 0, 0;
                                } 

                                VS->Offset = Offset;
                                
                                // prepare Multiindizes for coordinate data access
                                fi = 0,0,0;          fi  += Offset; 
                                lix = li[0], 0, 0;   lix += Offset; 
                                liy = 0, li[1], 0;   liy += Offset; 
                                liz = 0, 0, li[2];   liz += Offset;

                                // store MidPoint and radius of bounding sphere of each fragment, this info is going to be used for storing 
                                // fragment ids in a kd tree structure. Alternatively, one could use the less tight but existing bounding info of the fragment. 
                                // This should definitely be tried out (performance). 
                        tvector p000, p010, p110, p100, p001, p011, p111, p101; 


                                // if f == NullPtr data is gotten from the unfragmented field 
                                // f == NullPtr if there are no fragments 
                        RefPtr<ProcArray_t>   PCrds; 
                        RefPtr<CoordsArray_t> Crds; 

                        // try getting fragmented coords (using the frag id of the data field, which in case should always be consistent)
                                PCrds = S->Coords->getData(f); 
                                Crds  = S->Coords->getData(f); 

                                S->have_fragmented_coords = true; 

                                if(!f)
                                        S->have_fragmented_coords = false;      

                        // if no fragmented found try the unfragmented ones
                                if(!PCrds)
                                {
                                        PCrds = S->Coords->getData(); 
                                        if(PCrds)
                                                S->have_fragmented_coords = false;
                                } 
                                if(!Crds)
                                {
                                        Crds  = S->Coords->getData(); 
                                        if(Crds)
                                                S->have_fragmented_coords = false;
                                }
                                //cudaws01.cct.lsu.edu
                                if( Crds )
                                {
                                        p000 = tvector((*Crds)[fi]); 
                                        p010 = tvector((*Crds)[fi + liy]); 
                                        p110 = tvector((*Crds)[fi + lix + liy]); 
                                        p100 = tvector((*Crds)[fi + lix]); 
                                        p001 = tvector((*Crds)[fi + liz]); 
                                        p011 = tvector((*Crds)[fi + liy + liz]); 
                                        p111 = tvector((*Crds)[fi + lix + liy + liz]); 
                                        p101 = tvector((*Crds)[fi + lix + liz]); 

                                        S->have_curvi_coords = true;
                                } 
                                else
                                {
                                        if(!PCrds)
                                        {
                                                puts("no proc and no curvi coords found"); 
                                                return false;
                                        } 

                                        p000 = tvector((*PCrds)[fi]); 
                                        p010 = tvector((*PCrds)[fi + liy]); 
                                        p110 = tvector((*PCrds)[fi + lix + liy]); 
                                        p100 = tvector((*PCrds)[fi + lix]); 
                                        p001 = tvector((*PCrds)[fi + liz]); 
                                        p011 = tvector((*PCrds)[fi + liy + liz]); 
                                        p111 = tvector((*PCrds)[fi + lix + liy + liz]); 
                                        p101 = tvector((*PCrds)[fi + lix + liz]); 
                                        
                                        S->have_curvi_coords == false;
                                }

                                //cout << p000 << p010 << p110 << p100 << p001 << p011 << p111 << p101 << endl;

                                VS->MidPoint = point(0.125 * (p000 + p011 + p110 + p100 + p001 + p011 + p111 + p101)); 
                                VS->radius   = norm( point(p000) - VS->MidPoint ); 

                                S->max_radius = max( S->last_radius, VS->radius ); 
                                S->last_radius = VS->radius;

                                S->Blocks.push_back( VS ); 

                                cout << "ITERATE inserting " << VS->MidPoint << endl;

                                S->Tree->insert(VS->MidPoint, S->Blocks.size()-1 ); 

                                //S->Tree->speak(); 

                                c++;
                                
                                return true;
                        }
        } ; 
        

        static  inline double max(const double a, const double b)
        {
                return (a>b)?a:b;
        }
        
        static inline bool pointInBB(point& P, RefPtr<BoundingBox> BBox )
        {
                if( BBox->maxcoord()[0] > P[0] && BBox->mincoord()[0] <= P[0] &&
                    BBox->maxcoord()[1] > P[1] && BBox->mincoord()[1] <= P[1] &&
                    BBox->maxcoord()[2] > P[2] && BBox->mincoord()[2] <= P[2] )
                        return true;

                return false;
        }

        //norm 
        
        
        class BBSelectorBase
        {
        protected:
                RefPtr<MyState>S;
                
        public:
                BBSelectorBase(const RefPtr<MyState>&_S):S(_S) {assert(S);};
                virtual ~BBSelectorBase(){}
        
                virtual unsigned int getNext() = 0;
                virtual bool hasNext() = 0;
                virtual void reset(const int _last, const point& pos) = 0;
        };

        class BBSelectorMemoryAll : public BBSelectorBase
        {
                int size, current, last, count_up, count_down, count;
        public:

                BBSelectorMemoryAll(const RefPtr<MyState>_S) : BBSelectorBase(_S), current(0), last(0), count_up(1), count_down(0), count(0) 
                {
                        size = S->Blocks.size();
                }
                
                unsigned int getNext();
                bool hasNext();
                
                void reset(const int _last, const point& pos) { current = 0; last = _last; count = 0; count_up = 1; count_down = 0; }
        };
        

        class BBSelectorMemoryTree : public BBSelectorBase
        {
                int last, current, size;
                point reset_pos;
                map<double, int>m;
                map<double, int>::iterator it;
        public:
                BBSelectorMemoryTree(const RefPtr<MyState>_S) :  BBSelectorBase(_S), reset_pos(point(0,0,0))
                        {
                        m.clear(); 
                        last=0; 
                        current=0; 
                        size=0;
                }
                
                unsigned int getNext();
                bool hasNext();
                void reset( const int _last, const point& pos );
        };
        
        class CellTracing
        {
        protected:
                RefPtr<CoordsArray_t>          CoordsArr;
                RefPtr<VectorArray_t>          VecsArr;
                MultiIndex<3> iAx, iAy, iAz; 
        
        public:
                CellTracing()
                        {
                                iAx = 1, 0, 0; 
                                iAy = 0, 1, 0; 
                                iAz = 0, 0, 1;
                        }
                
                virtual ~CellTracing(){}
                
                void setCoords(const RefPtr<CoordsArray_t>  CA)  { CoordsArr = CA; }
                void setVecs  (const RefPtr<VectorArray_t>  VA)  { VecsArr   = VA; }

                inline tvector cross(const tvector&t1, const tvector&t2)
                        {
                        tvector res; 
                                res =  t1[1]*t2[2] - t1[2]*t2[1],
                                       -( t1[0]*t2[2] - t1[2]*t2[0] ),
                                       t1[0]*t2[1] - t1[1]*t2[0]; 

                                return res;
                        }

                virtual int    pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines) = 0;
                virtual point  localUVW   (const point&p, const int debug, std::vector<colorLine_t>& debuglines) = 0;
        };

        class HexaHedralCellNewton : public CellTracing
        {
                RefPtr<UniGridMapper> UniMap;

                point uvw;
                int point_type;
        public:
                HexaHedralCellNewton():CellTracing()
                {}

 
                virtual ~HexaHedralCellNewton(){}

                void  setUniMap  ( RefPtr<UniGridMapper> unimap ) { UniMap = unimap; }
                int   pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines);
                point localUVW   (const point&p, const int debug, std::vector<colorLine_t>& debuglines);
        };



        class HexaHedralCellXYTracing : public CellTracing
        {

        // objects calced in pointInCell() for use in localUVW()  
        std::vector<point>          interval_xy; 
        std::vector<MultiIndex<3> > interval_index_xy, 
                                    interval_index_xy_perm;                             
        MultiIndex<3> Z, uZ; 
        public:
                HexaHedralCellXYTracing()
                {
                        interval_xy.resize(4); 
                        interval_index_xy.resize(4); 
                        interval_index_xy_perm.resize(4);
                }
                virtual ~HexaHedralCellXYTracing(){}

                static  inline double dist_xy(const point& p1, const point& p2)
                {
                        return sqrt( (p2[0]-p1[0])*(p2[0]-p1[0]) + (p2[1]-p1[1])*(p2[1]-p1[1]) );
                }
        
                // intersect line segment and half-line
                int testIntersect(const point& p1, const point& p2, const point& p, const tvector& d);
                int intersect(const point& p1, const point& p2, const point& p, const tvector& d, double*lambda,  double*mue); //costs // + 7 * 10 
                int pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines);
                
                // method must only be called after successful pointInCell test!
                point localUVW   (const point&CurrentPoint, const int debug, std::vector<colorLine_t>& debuglines); 





        };      
        

        override RefPtr<State> newState() const
        {
                return new MyState();
        }

        // prepare caching of computed lines
        typedef VOperatorCache< ::IntegralLines::LinesPoints > OpCache_t;
        
        TypedSlot<int> LineLength;
        TypedSlot<double> StepSize;
        TypedSlot< ::IntegralLines::EmitterPoints> PointsInput;

        TypedSlot<double>       SliceLocation;

        TypedSlot<Enum>         ScalingOptions;
        TypedSlot<Enum>         Debug;
        TypedSlot<Enum>         UVWType;
        TypedSlot<Enum>         Shifter;
        TypedSlot<Enum>         TreeAcc;

        TypedSlot<double>       Scale,
                                Thickness,
                                ArrowAngle,
                                RelativeSize,
                                AtanLength,
                                TreeQueryFac;

        TypedSlot<int>          FirstFragment,
                                LastFragment;

        TypedSlot<int>         nX, nY, nZ;

        VOutput< ::IntegralLines::LinesPoints> LinesOutput;

//      TypeList_t      AcceptedVectorfieldTypes;

        VOutput<BundlePtr>      myBundle;
        VOutput<Grid>           myIntegralLines;

        ComputeMultiStreamLines(const string&name, int p, const RefPtr<VCreationPreferences>&VP)
        : VRenderObject(name, p, VP) 
        , Fish<VObject>(this)
        , Fish<Field>("datafield")
        , LineLength(this, "linelength", 40)
        , StepSize(this, "stepsize", 0.02)
        , PointsInput(this, "points_input", ::IntegralLines::EmitterPoints())
        , SliceLocation(this, "slicecoords", 0.0, 3)
        , ScalingOptions(this,"normalize", Enum("linear", "norm", "atan"), 3)
        , Debug(this, "debug", Enum("off", "on"))
        , UVWType(this, "uvw_type", Enum("XY", "Newton"))
        , Shifter(this, "shifter", Enum("off", "on"))
        , TreeAcc(this, "kdtree", Enum("off", "on"))
        , Scale        (this, "scale"      , 0.0, 3)
        , Thickness    (this, "thickness"  , 1.0  ,NullPtr() , 2)
        , ArrowAngle   (this, "arrowangle" , 45.0 ,NullPtr() , 3)
        , RelativeSize (this, "headsize"   , 0.3  ,NullPtr() , 3)
        , AtanLength   (this, "ascale"     , 1.0  ,NullPtr(), 3)
        , TreeQueryFac(this, "treefac", 1.0)
        , FirstFragment(this, "firstfragment",   0)
        , LastFragment (this, "lastfragment" ,  10, 4)
        , nX(this, "iX",   0, 2)
        , nY(this, "iY",   0, 2)
        , nZ(this, "iZ",   0, 2)
        , LinesOutput(self(), "line_output", ::IntegralLines::LinesPoints() ) 

        , myBundle(self(), "outspacetime", BundlePtr())
        , myIntegralLines(self(), "integrallines", GridSelector() )
        {
                printf("ComputeMultiStreamLines: I'm Here!!\n");fflush(stdout);
                FirstFragment.setProperty("max", 2087);
                LastFragment.setProperty("max", 4000); 

                LineLength.setProperty("max", 500);

                Scale.setProperty( "min",-10.0);
                Scale.setProperty( "max", 10.0);

                Thickness .setProperty( "max", 10.0);
                ArrowAngle.setProperty( "max", 90.0); 

                AtanLength.setProperty( "min",-10.0);
                AtanLength.setProperty( "max", 10.0); 

//              AcceptedVectorfieldTypes.insert( FiberType<Eagle::tvector3>::getFiberType().self() ); 
                acceptType<Eagle::tvector3> ();
        }

        override bool update(VRequest&R, double precision);
        override void render(VRenderContext&Context) const;     
};

//#define VERBOSE


bool ComputeMultiStreamLines::update(VRequest&Context, double precision)
{
        puts("update");fflush(stdout);
VTimer Tim_update; 
double secs; 

        printf("---- ComputeMultiStreamLines::update\n");fflush(stdout); 
        if (!updateSelection(Context) )
        {
                puts("ComputeMultiStreamLines::update(): No field!"); 
                removeState(Context); 
                return true;
        } 

RefPtr<MyState> S = myState(Context); 

// BundlePtr BP; 
//      SpacetimeSlot() << Context >> BP; 
//      pair<double, RefPtr<Fiber::Grid> > grid = getGrid( Context ); 
//      if(grid.second)
//              printf("------------------ found grid: \n" ); 
        

int     ML = getMaxLevel( Context );

BBoxExpander BBE( getBoundingBall( Context ) ); 

        for(int rl = 0; rl <= ML; rl++)
        {
        RefPtr<Skeleton> Level = getRefinementLevel(SpacetimeSlot(), Context, rl).second;
                if (!Level) 
                {
                        printf("ComputeMultiStreamLines: update(): No level found!\n"); 
                        break;
                }
        RefPtr<Representation> LevelRep = Level->DefaultRepresentation(); 
        RefPtr<Field> Coords  = LevelRep->Positions();
                if (!Coords)
                {
                        puts("No coord field!");
                        return true;
                } 
                S->Coords = Coords;
                BBE.expand( Coords );

        RefPtr<Field> F = getField( LevelRep, Context );
                if (!F)
                {
                        printf("No field for [%s]!\n",
                               Fieldname( Context ).c_str() );
                        return true; 
                }
                
                 // We use the field in the finest refinement level for opengl
                 // cacheing, assuming that this is the one changing fastest
                 // in time.
                 
                S->F = F; 

                if(S->F->nFragments() > 1)
                        S->have_fragmented_data = true; 
                else
                        S->have_fragmented_data = false;

                S->RefinementLevel = rl; 

                printf("We have %d Fragments here!\n",  S->F->nFragments() );

                // skip line if already done before, iterate over all Blocks and build KDTree 
                if(S->Blocks.size() == 0)
                {
                bool test;
                It MyIt(S); 

                        puts("before iterate");fflush(stdout);
                        ________________________________________________getMTime(secs, Tim_update);
                        test = F->iterate(MyIt);
                        ________________________________________________printMTime("Coords iteration", Tim_update.secs() - secs); 

                        if(!test)
                        {
                                puts("iteration failure, abording update"); 
                                return false;
                        }
                        printf("xxxxxxxxxxxxxxxxxxxxxxxxx heAr the tree speAK xxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n");fflush(stdout); 
                        S->Tree->speak();
                } 
                
        } 

        //S->adjustSliceRange(SliceLocation);
        setBoundingBall( Context, BBE.BBox); 

        // use caching field for caching the lines 
        RefPtr<OpCache_t>&OCLines = OpCache_t::retrieve( *S->F ) [ self() ]; 
        if (!OCLines) OCLines = new OpCache_t(); 

        ::IntegralLines::EmitterPoints StartPoints; 
        PointsInput << Context >> StartPoints; 

RefPtr<ValueSet> Changeables = new ValueSet(); 
Changeables & Context += LineLength, StepSize, PointsInput, Debug, UVWType, TreeAcc, TreeQueryFac; 

        if (OCLines->unchanged( Changeables ) )
        {
                printf("ComputeStreamLines::update(): Line Calc Skipped, Lines in Cache...\n"); 
        ::IntegralLines::LinesPoints tmp = (*OCLines); 
                LinesOutput << Context << tmp; 
                return true;
        }       

        OCLines->Lines.resize( StartPoints.Vertices.size() ); // resize std::vector to number of lines 
        GLCache::clear( *OCLines ); 

        // if no emitter points available -> quit
        if( StartPoints.Vertices.size() <= 0 ) return true; 

        //printf("ComputeMultiStreamLines:update(): StartPoint[0]: %f %f %f\n", StartPoints.Vertices[0][0], StartPoints.Vertices[0][1], StartPoints.Vertices[0][2]); 

        S->lines.clear(); //these are just lines drawn for debugging 



string grid_name = getGridname(Context) + "_ComputeMultiStreamLines"; 

pair< double, RefPtr< Fiber::Slice > >
        ActualSlice = getSlice (SpacetimeSlot(), Context); 

// create a new Grid in the actual slice of the bundle 
//Grid&LineGrid = (*ActualSlice.second)[grid_name]; 
Grid&LineGrid = ActualSlice.second->newGrid(grid_name); 

// Create skeleton for the vertices 
Skeleton&LineVertices = LineGrid.makeVertices(3); 

// Create a skeleton ID of dimension 1, and index depth 1
Skeleton&LineEdges = LineGrid[ SkeletonID(1,1) ]; 

RefPtr<Chart> myCartesian = LineGrid.makeChart( typeid(CartesianChart3D) ); 
Representation&CartesianVertices = LineVertices[ myCartesian ]; 
Representation&EdgesAsVertices   = LineEdges[ LineVertices ]; 

RefPtr<Field> LineCoords    = CartesianVertices[ FIBER_POSITIONS ]; 
RefPtr<Field> LineSetEdges  = EdgesAsVertices  [ FIBER_POSITIONS ]; 

MultiIndex<1> NumberOfLineVertices(0);


        {
        BundlePtr BPtr;
                findSpacetimeSlot() << Context >> BPtr; 
                myBundle << Context << BPtr; 

        GridSelector GS( new StringSelection( grid_name ) ); 
                myIntegralLines << Context << GS;
        }




double step_size; 
Enum Deb, UV, Shift, TreeA; 
int line_length = 0; 

        LineLength << Context >> line_length;
        StepSize   << Context >> step_size; 
        Debug      << Context >> Deb; 
        UVWType    << Context >> UV; 
        Shifter    << Context >> Shift; 
        TreeAcc    << Context >> TreeA; 
        TreeQueryFac<< Context >> S->tree_fac;

HexaHedralCellXYTracing HexaXY; 
HexaHedralCellNewton    HexaNewton; 

BBSelectorMemoryAll  BBSelectAll( S ); 
BBSelectorMemoryTree BBTreeSel( S ); 

BBSelectorBase* BBSelect; 

        if(TreeA() == 0)
                BBSelect = &BBSelectAll; 
        else
                BBSelect = &BBTreeSel;

                

int last_cand = 0; 

RefPtr<VectorArray_t> VecsArr; 

double bound_secs = 0.0; 
double tree_secs = 0.0; 
double alllines_secs; 
        ________________________________________________getMTime(alllines_secs, Tim_update); 

unsigned int nth =  StartPoints.Vertices.size()/10; 
        nth = (nth == 0)?5:nth;

for( unsigned int vs = 0; vs < StartPoints.Vertices.size(); vs++ ) 
{
unsigned percent; 
char buf[16];
        if( vs % nth == 0)
        {
                percent = (unsigned)ceil( (double)vs / StartPoints.Vertices.size() * 100); 
                sprintf(buf, "%d %% done", percent); 
                VActionNotifier::ProgressInfo (string(buf)); 
        }
 
Line_t& myLine = (*OCLines).Lines[vs]; 
        myLine.clear(); 

        printf("ComputeMultiStreamLines:update(): Entering StartPointvertex: %d ---------------------------------------------- \n", vs); fflush(stdout); 

        S->candidates.clear(); 


point CurrentPoint = StartPoints.Vertices[vs]; 

        myLine.push_back(CurrentPoint); 

        bool abreis_shifted = false;

        for(int len = 0; len < line_length; len++) 
        {
                ________________________________________________getMTime(secs, Tim_update);
        int real_candidate = 0; 
        int real_cand_index = 0; 
        int cell_found = 0; 
        
        //reset block selector 
double tree_tmp_secs, tmp_secs; 
        ________________________________________________getMTime(tree_tmp_secs, Tim_update);
        BBSelect->reset( last_cand, CurrentPoint ); 
        ________________________________________________getMTime(tmp_secs, Tim_update);
        tree_secs += (tmp_secs - tree_tmp_secs); 

//      printf("Has Next %d\n", BBSelect->hasNext()  );

        while(  BBSelect->hasNext() )
        {
        int j = BBSelect->getNext(); 

                real_cand_index = j; 

                //printf("getNext() %d\n", j);fflush(stdout); 

                //assert(S->Blocks[j]->FragID);

        // get the bounding box of the fragment
        RefPtr< CreativeArrayBase > MetaData = S->F->getCreator( S->Blocks[j]->FragID ); 
        assert( MetaData ); 

double bound_tmp_secs; 
        ________________________________________________getMTime(bound_tmp_secs, Tim_update); 
        RefPtr<BoundingBox> BBox; 
                BBox= getFragmentBBox(S->Blocks[j]->FragID, MetaData, *(S->Coords) );
                assert( BBox );
                ________________________________________________getMTime(tmp_secs, Tim_update); 
                bound_secs += (tmp_secs - bound_tmp_secs);


                VecsArr = S->F->getCreator( S->Blocks[j]->FragID )->create(); 
                assert( VecsArr ); 


                //printf("status of state: have_curvi=%d, have_fragmented_coords=%d, have_fragmented_data=%d\n, ", S->have_curvi_coords, S->have_fragmented_coords, S->have_fragmented_data);

                
                // check if point is inside the selected BBox 
                if (BBox->inside(CurrentPoint))
                {                                 
                        S->candidates.push_back(j); 
                        //printf("ComputeMultiStreamLines:update(): index: %d\n", j);fflush(stdout); 

                        if( S->have_curvi_coords )
                        {
                                RefPtr<CoordsArray_t> CoordsArr = S->Coords->getData(S->Blocks[j]->FragID); 
                                assert( CoordsArr ); 

                                // prepare and invoke search for the cell in the selected block 
                                if( UV() == 0 )
                                {
                                        HexaXY.setCoords( ( CoordsArr ) ); 
                                        HexaXY.setVecs  ( ( VecsArr) ); 
                                        cell_found = HexaXY.pointInCell(CurrentPoint, Deb(), S->lines); 
                                } 
                                else if( UV() ==1 )
                                {
                                        HexaNewton.setCoords( (CoordsArr) ); 
                                        HexaNewton.setVecs  ( (VecsArr) ); 

                                double uni_secs; 
                                        ________________________________________________getMTime(uni_secs, Tim_update); 
                                        if( !(S->Blocks[ real_cand_index ]->UniMap) )
                                        {
                                                S->Blocks[ real_cand_index ]->UniMap = new UniGridMapper( (CoordsArr) ); 
                                        } 
                                        ________________________________________________printMTime("Unigridmapper:", Tim_update.secs() - uni_secs); 
                                        HexaNewton.setUniMap( (S->Blocks[ real_cand_index ]->UniMap)  ); 
                                        cell_found = HexaNewton.pointInCell(CurrentPoint, Deb(), S->lines); 
                                } 
                        } 
                        else // do the search for a procedural coords
                        {
                                // direct product arrays are alsways aligned to world xyz coordinates, thus
                                // a positive check of the bounding box guaratnies to be in the coordinate range 
                                // of the block.
                                cell_found = true;
                        }

                        last_cand = j;

                        if( cell_found )
                        {
                                break;
                        } 
                } 
        } 


        ________________________________________________printMTime("Find real cand", Tim_update.secs() - secs); 



        ________________________________________________getMTime(secs, Tim_update);

//      printf("cell_found: %d in %d\n", cell_found, real_cand_index); 

        // compute UVW coords and new line point, if a cell was found 
        if( cell_found )
        {
        point float_index; 

                abreis_shifted = false; 

                if(S->have_curvi_coords)
                {
                        if(UV() == 0)
                                float_index = HexaXY.localUVW(CurrentPoint, Deb(), S->lines); 
                        else if(UV() == 1)
                                float_index = HexaNewton.localUVW(CurrentPoint, Deb(), S->lines);
                } 
                else
                {
                        RefPtr<ProcArray_t> PCrd;
                        if( S->have_fragmented_coords )
                                PCrd = S->Coords->getData( S->Blocks[last_cand]->FragID ); 
                        else
                                PCrd = S->Coords->getData();
                        assert(PCrd); 

                        float_index = CurrentPoint; 

                bool test = PCrd->locate(float_index); 
                        assert(test);

                        // if unfragged coords correct float index by Offset of fragment
                if( !S->have_fragmented_coords && S->have_fragmented_data)
                for(unsigned i = 0; i < 3; i++)
                        float_index[i] -=  S->Blocks[last_cand]->Offset[i];

                }

                Interpolate<3, tvector, LinearIpol  <tvector> > VecField( *(VecsArr), float_index); 
                tvector DataDir = VecField.eval(); 

                        //std::cout << "current point: " << CurrentPoint << std::endl; 
                        //std::cout << "direction:"      << DataDir      << std::endl;

                point OldPoint = CurrentPoint;

                CurrentPoint = CurrentPoint + step_size * DataDir; 

                if(Deb() == 1) { S->lines.push_back( colorLine_t( line_t( OldPoint, CurrentPoint ), RED) ); } 

                real_candidate = 0; 

                myLine.push_back(CurrentPoint); 

                printf("ComputeMultiStreamLines:update(): Proceeding to next line point %d in line %d\n\n", len, vs); 

        }       
        else
        {
                // Step a half step back if no point was found to ensure, stepping over tight block borders.
                // This does not seem neccessary at all after correction of an index bug in the UniGridMapper...
                if( UV() == 1 && abreis_shifted == false && Shift() == 1)
                {
                        abreis_shifted = true; 

                        printf("Trying little shift forward instead!\n"); 

                point a, b; 
                unsigned max = myLine.size();
                        b = myLine[max-1]; 
                        a = myLine[max-2]; 

                tvector dir = b-a; 

                        b = a + dir * 0.5; 
                        std::cout<< a <<" + 1.2 * " << dir << std::endl; 

                        myLine.pop_back();
                        myLine.push_back(b);
                        CurrentPoint = b; 
                        std::cout<< "CerrentPoint " << CurrentPoint << std::endl; 

                        real_candidate = 0;
                        
                }
                else
                {
                        last_cand = 0; 
                        printf("ComputeMultiStreamLines:update(): No real candidate found, Breaking Line\n"); 
                        break; 
                }
        } 



        ________________________________________________printMTime("Compute cand", Tim_update.secs() - secs); 

        } //for len 
} //for v 
________________________________________________printMTime("BoundingBox access", bound_secs);
________________________________________________printMTime("TreeQuery access", tree_secs); 
________________________________________________printMTime("All lines comp", Tim_update.secs() - alllines_secs); 

::IntegralLines::LinesPoints&computed_lines = (*OCLines);  

        printf("ComputeStreamLines::update(); computed_lines size %d\n", computed_lines.Lines.size()); 

        LinesOutput << Context << computed_lines; 

        //printf("All Lines Done\nDeleting empyt lines(horrible!)\n"); 
        for( std::vector< std::vector<point> >::iterator it = computed_lines.Lines.begin(); it != computed_lines.Lines.end(); it++ )
        {
                if( it->size() == 1)
                {
                        it->push_back((*it)[0]);
                } 
                if( it->size() == 0)
                {
                        it->push_back(point(0,0,0)); 
                        it->push_back(point(0,0,0));
                }
        } 



        ________________________________________________getMTime(secs, Tim_update);

        for( index_t ii = 0; ii < computed_lines.Lines.size(); ii++ )
        {
                MultiIndex<1> computed_linesInd( computed_lines.Lines[ii].size() ); 
                NumberOfLineVertices += computed_linesInd; 
                printf("line %d, has %d verteces\n", ii,  computed_lines.Lines[ii].size());
        } 

        printf( "Number of Vertices: %d\n", NumberOfLineVertices[0] ); fflush(stdout);

        typedef MemArray<1, point> LinesCoordsArr_t; 
        RefPtr<LinesCoordsArr_t> LinesCoordsData = new LinesCoordsArr_t(NumberOfLineVertices); 
        MultiArray<1,point > LineCrds = *LinesCoordsData; 
        {
        index_t vert_count = 0;
                for(index_t lin = 0; lin <  computed_lines.Lines.size(); lin++ )
                {
                        printf("\nEntering line %d of %d\n", lin, computed_lines.Lines.size());fflush(stdout); 
                        for(index_t ver = 0; ver < computed_lines.Lines[lin].size() ; ver++ )
                        {
                                printf(" V(%d %d %d) ", ver, computed_lines.Lines[lin].size(), lin);fflush(stdout); 
                                printf(" C(%f %f %f)\n", computed_lines.Lines[lin][ver][0], computed_lines.Lines[lin][ver][1], computed_lines.Lines[lin][ver][2]);
                        MultiIndex<1> n( vert_count );  
                                
                                LineCrds[ n ] =  computed_lines.Lines[lin][ver];
                                vert_count++; 
                        } 
                } 
                //printf("\n");
        } 
        LineCoords->setPersistentData( LinesCoordsData ); 

//        printf("Coords done!\n");

        typedef MemArray<1, std::vector<index_t> > EdgesArray_t; 
        RefPtr<EdgesArray_t> EdgesArray = new EdgesArray_t( computed_lines.Lines.size()  ); 
        MultiArray<1, std::vector<index_t> >&LinesEdges = *EdgesArray; 
        {
        index_t vert_count = 0;
                for(index_t lin = 0; lin <  computed_lines.Lines.size(); lin++ )
                {
                        //printf("\nEntering line %d of %d\n", lin, computed_lines.Lines.size());fflush(stdout); 
                        std::vector<index_t>&E = LinesEdges[ lin ]; 
                        E.resize( computed_lines.Lines[lin].size() ); 
                        //printf("2\n");fflush(stdout);
                        for(index_t ver = 0; ver < computed_lines.Lines[lin].size() ; ver++ )
                        {
                                //printf("(E(%d %d)\n", ver, vert_count);fflush(stdout); 
                                E[ver] = vert_count; 
                                vert_count++; 
                        } 
                } 
        }
  
        LineSetEdges->setPersistentData( EdgesArray ); 

        
        ________________________________________________printMTime("Copy in Bundle", Tim_update.secs() - secs); 

        // the bundle data object has to be touched, to be updated correctly
        if(RefPtr<VObject> vobj = SpacetimeSlot()->Source() )
                vobj->touch(); 


        ________________________________________________printMTime("update Time", Tim_update.secs()); 


        printf("ComputMultipleStreamLines::update() Written into Bundle, setPersitentData and touched Spacetimeslot()\n");
        VActionNotifier::ProgressInfo (string("Streamline computation complete"));
        return true; 
}


static  inline void glLine(const ComputeMultiStreamLines::point&A, const ComputeMultiStreamLines::tvector&dir)
        {
//              glNormal( *dir  );
                glVertex( A     );
                glVertex( A+dir );
        }


void ComputeMultiStreamLines::render(VRenderContext&Context) const
{
/*
Enum Deb; 

        Debug << Context >> Deb;

        if(Deb() == 1)
        {
                
                
                RefPtr<FieldState> state = getState(Context); 
                if (!state)
                        return; 

                if (state->Blocks.size() < 1)
                {
                        puts("override void ComputeMultiStreamLines::render(VRenderContext&Context) const: NOTHING to render..."); 
                        return;
                } 

                
                assert(state); 


        int count = 0; 
                glBegin(GL_POINTS); 
                for( std::vector<colorPoint_t>::iterator it = state->points.begin(); it != state->points.end(); it++)
                {
                        glColor3f((it->second)[0], (it->second)[1], (it->second)[2]); 
                        glVertex(it->first); 
                        count++;
                } 
                glEnd(); 

                //printf("%d points drawed\n", count); 

                glBegin(GL_LINES); 
                for( std::vector<colorLine_t>::iterator it = state->lines.begin(); it != state->lines.end(); it++)
                {
                        glColor3f((it->second)[0], (it->second)[1], (it->second)[2]); 
                        glVertex((it->first).first); 
                        glVertex((it->first).second); 
                } 
                glEnd();                
        }
        // render candidate boxes, only candidates of one line are rendered, of the last one. 
        // So, if the last line is not in data range (only emit point visible) no candidate boxes 
        // are visible. 
        // for(unsigned int i = 0; i < myState->candidates.size(); i++ )
//      {
//              RefPtr<BoundingBox> BBox = myState->Blocks[myState->candidates[i]]->BBox; 

//              glColor3f(0.4,0.4,0.4); 
//              glDrawBox( point(BBox->maxcoord()), point(BBox->mincoord()) ); 
//      } 

        //printf("Render: draw box ? sizes: %d %d\n", myState->cand_coord_cand.size() , myState->cand_coord.size()); 


        // render cand_intervals 
        //printf("ComputeMultiStreamLines:render(): sizes: cand_coord %d, cand_coord_cand %d\n", myState->cand_coord.size(), myState->cand_coord_cand.size()); 

//              for(unsigned int j = 0; j < myState->cand_coord_cand.size(); j = j+4 ) 
//              {
//                      MultiIndex<3> m1, m2, m3, m4; 
//                              m1 = myState->cand_coord[j];  
//                              m2 = myState->cand_coord[j+1]; 
//                              m3 = myState->cand_coord[j+2]; 
//                              m4 = myState->cand_coord[j+3]; 
                                
//                              //printf("Render: draw box ix %d iy  %d in %d\n", m1[0], m1[1], myState->cand_coord_cand[j]); 
//                              glBegin( GL_LINES ); 
//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m1 ]  ); 
//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m2 ]  ); 

//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m2 ]  ); 
//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m3 ]  ); 

//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m3 ]  ); 
//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m4 ]  ); 

//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m4 ]  ); 
//                              glVertex( (*(myState->Blocks[ myState->cand_coord_cand[j] ]->CoordsArr))[ m1 ]  ); 
//                              glEnd(); 
//                      } 
                
// 
                
//                      glLineWidth(2);
                        
//                      glBegin( GL_LINES );
//                      for(unsigned int j = 0; j < myState->lines.size(); j+=2 )
//                      {
//                              glColor3f(0, 1, 1); 
//                              glVertex ( myState->lines[j]); 
//                              glVertex ( myState->lines[j+1]);
//                      }

//                      glEnd(); 

//                      glPointSize(4); 
//                      glColor3f(1.0, 1.0, 0.0); 

//                      glBegin( GL_POINTS );
//                      for(unsigned int j = 0; j < myState->lines.size(); j+=2 )
//                      {
//                              glColor3f(0, 1, 1); 
//                              glVertex ( myState->lines[j]); 
//                      } 
//                      glEnd(); 


//               int nx, ny, nz; 


//                       nX << Context >> nx; 
//                       nY << Context >> ny; 
//                       nZ << Context >> nz;

//                      MultiIndex<3> bla; 

//                      bla = nx, ny, nz; 

//                      glColor3f(0, 1.0, 1.0); 
//                      glBegin(GL_POINTS);
//                      glVertex( (*(myState->Blocks[ First ]->CoordsArr))[ bla ]  ); 
//                      glEnd();


//                      // glBegin( GL_POINTS );
// //                   for(unsigned int j = 0; j < myState->points.size(); j++ )
// //                   {
// //                           glColor3f(1, 1, 0); 
// //                           glVertex(myState->points[j]);
// //                   } 
// //                   glEnd();
                        
// //   puts(" render ok");
//      } 

//      glDisable(GL_LINE_SMOOTH); 
//      }
*/
}


int ComputeMultiStreamLines::HexaHedralCellXYTracing::testIntersect(const point& p1, const point& p2, const point& p, const tvector& d)
{
        tvector d_seg = p2 - p1; 
        double lambda = 0; 
        double mue = 0; 

                lambda = ( p[1]*d[0] + p1[0] * d[1] - p[0] * d[1] - p1[1]*d[0]) / ( d_seg[1]*d[0] - d_seg[0]*d[1] ); 
                if( lambda < 0.0 || lambda > 1.0 )
                        return 0; 

                 mue = (p1[0] + lambda * d_seg[0] - p[0]) / d[0]; 
                 if( mue < 0.0 )
                         return 0; 

                 return 1;
}

int ComputeMultiStreamLines::HexaHedralCellXYTracing::intersect(const point& p1, const point& p2, const point& p, const tvector& d, double*lambda,  double*mue)
{
tvector d_seg = p2 - p1; 
        *lambda = 0; 
        *mue = 0; 

        *lambda = ( p[1]*d[0] + p1[0] * d[1] - p[0] * d[1] - p1[1]*d[0]) / ( d_seg[1]*d[0] - d_seg[0]*d[1] ); 
        if( *lambda < 0.0 || *lambda > 1.0 )
                return 0; 

        *mue = (p1[0] + *lambda * d_seg[0] - p[0]) / d[0]; 
        if( *mue < 0.0 )
                return 0; 

        return 1;       
}

int ComputeMultiStreamLines::HexaHedralCellXYTracing::pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines)
{
        MultiIndex<3> FieldSize = CoordsArr->Size(); 


        //printf("ComputeMultiStreamLines:update(): Candidate size: %d %d %d\n", FieldSize[0], FieldSize[1], FieldSize[2] ); 
        MultiIndex<3> ifirst, ilast, ifx, ify, ifz; 

        point float_index, max_float_index, p1, p2; 

        unsigned int ix, iy; 

                max_float_index = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1; 

                ilast = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1; 
                ifirst = 0,0,0; 
                ifx = 1, 0, 0; 
                ify = 0, 1, 0;
                ifz = 0, 0, 1; 

                p1 = (*CoordsArr)[ifirst]; 
                p2 = (*CoordsArr)[ilast]; 

        tvector vx, vy, vz; 

                vx = (*CoordsArr)[ifx] - p1; 
                vy = (*CoordsArr)[ify] - p1; 
                vz = (*CoordsArr)[ifz] - p1; 
                // if ilast... < 1 --> break 

                // printf("vx: %f %f %f \n", vx[0], vx[1], vx[2]); 
//                      printf("vy: %f %f %f \n", vy[0], vy[1], vy[2]); 
//                      printf("vz: %f %f %f \n", vz[0], vz[1], vz[2]); 

                // here we are looking for the index that represents the z-Axis, since we are only looking at 
                // the xy coords for the cell finding later. 

                // re-index-index Z and un-re-index-index uZ for reindexing und unreindexing z-index to the 3rd index,of MultiIndizes, 
                // which happens not to be the index for the z-direction generally.
                if( vx[2] != 0.0)
                {
                        Z  = 2, 1, 0; 
                        uZ = 2, 1, 0; 
                        //printf("x is z\n");   
                } 
                else if( vy[2] != 0.0)
                {
                        Z  = 2, 0, 1; 
                        uZ = 1, 2, 0; 
                        //printf("y is z\n");
                }
                else if( vz[2] != 0.0)
                {
                        Z = 0, 1, 2; 
                        uZ = 0, 1, 2; 
                        //printf("z is z\n");   
                } 
                else
                {
                        Z = 0, 1, 2; 
                        uZ = 0, 1, 2; 
                        //printf("nothing is z\n");     
                } 

                //printf("p1: %f %f %f, p2: %f %f %f\n", p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]); 
 
                
        tvector testdir(0.77753, 0.2274, 0.111); 

                //tvector testdir(0.111, 0.45, 0); 
                // only x and y coords are crucial to decide whether the point is really contained (cylinrical coords) 
                // find by intersection of an line in z-plane intersection the bouning rectangle of a coordinate interval in the z-plane 

                //printf("Testing Intervals-------: %d %d %d\n", ilast[0], ilast[1], ilast[2]); 

                //printf("Flipped Intervals-------: %d %d %d\n", ilast[Z[0]], ilast[Z[1]], ilast[Z[2]]); 

         int real_candidate;

                for(ix = 0; ix < ilast[ Z[0]  ]; ix++)
                {
                        for( iy = 0; iy < ilast[ Z[1] ]; iy++)
                        {
                                //printf("Testing Interval: %d of %d, %d of %d\n", ix, ilast[ Z[0]],  iy, ilast[ Z[1] ] );
                                //fflush(stdout); 
                                real_candidate = 0; 

                                interval_index_xy[0] = ix, iy, 0; 
                                interval_index_xy[1] = ix+1, iy, 0; 
                                interval_index_xy[2] = ix+1, iy+1, 0; 
                                interval_index_xy[3] = ix, iy+1, 0; 

                                //printf("1");fflush(stdout); 
                                
                                //for(int i = 0; i < 4; i++) 
                                interval_index_xy_perm[0] = interval_index_xy[0][ uZ[0] ], interval_index_xy[0][ uZ[1] ], interval_index_xy[0][ uZ[2] ]; 
                                interval_index_xy_perm[1] = interval_index_xy[1][ uZ[0] ], interval_index_xy[1][ uZ[1] ], interval_index_xy[1][ uZ[2] ]; 
                                interval_index_xy_perm[2] = interval_index_xy[2][ uZ[0] ], interval_index_xy[2][ uZ[1] ], interval_index_xy[2][ uZ[2] ]; 
                                interval_index_xy_perm[3] = interval_index_xy[3][ uZ[0] ], interval_index_xy[3][ uZ[1] ], interval_index_xy[3][ uZ[2] ]; 

                                //printf("2");fflush(stdout); 

                                //for(int i = 0 ; i < 4 ; i ++) 
                                //      printf("permindex[%d]: %d %d %d\n", i, interval_index_xy_perm[i][0], interval_index_xy_perm[i][1], interval_index_xy_perm[i][2]); 

                                interval_xy[0] = (*CoordsArr)[interval_index_xy_perm[0]];  //t1 
                                interval_xy[1] = (*CoordsArr)[interval_index_xy_perm[1]];  //t2 
                                interval_xy[2] = (*CoordsArr)[interval_index_xy_perm[2]];  //t3 
                                interval_xy[3] = (*CoordsArr)[interval_index_xy_perm[3]];  //t4 

                                //printf("3");fflush(stdout); 

                                real_candidate += testIntersect(interval_xy[0], interval_xy[1], P, testdir); 
                                real_candidate += testIntersect(interval_xy[1], interval_xy[2], P, testdir); 
                                real_candidate += testIntersect(interval_xy[2], interval_xy[3], P, testdir); 
                                real_candidate += testIntersect(interval_xy[3], interval_xy[0], P, testdir); 

                                //printf("4");fflush(stdout);   

                                // S->lines.push_back( interval_xy[0]  ); S->lines.push_back( interval_xy[1] ); 
//                              S->lines.push_back( interval_xy[1]  ); S->lines.push_back( interval_xy[2] ); 
//                              S->lines.push_back( interval_xy[2]  ); S->lines.push_back( interval_xy[3] );
//                              S->lines.push_back( interval_xy[3]  ); S->lines.push_back( interval_xy[0] ); 

                                
                                if(real_candidate == 1)
                                { 
                                        //printf("REAk CAND found at index %d %d %d\n", interval_index_xy_perm[0][0], interval_index_xy_perm[0][1], interval_index_xy_perm[0][2]); 

                                        if(debug == 1)
                                        {
                                                debuglines.push_back( colorLine_t( line_t( interval_xy[0], interval_xy[1]), GREEN) ); 
                                                debuglines.push_back( colorLine_t( line_t( interval_xy[1], interval_xy[2]), GREEN) ); 
                                                debuglines.push_back( colorLine_t( line_t( interval_xy[2], interval_xy[3]), GREEN) ); 
                                                debuglines.push_back( colorLine_t( line_t( interval_xy[3], interval_xy[0]), GREEN) ); 
                                        } 

                                        return 1;

                                } 
                        } 

                        
                } 

        return 0;       
        }
                
point ComputeMultiStreamLines::HexaHedralCellXYTracing::localUVW   (const point&CurrentPoint, const int debug, std::vector<colorLine_t>& debuglines)
{
        if( !VecsArr)
                return point(0,0,0); 

        if( !CoordsArr)
                return point(0,0,0); 

        if( interval_xy.size() != 4)
                return point(0,0,0); 

        if( interval_index_xy.size() != 4)
                return point(0,0,0); 

        if( interval_index_xy_perm.size() != 4)
                return point(0,0,0);

        //printf("ComputeMultiStreamLines:update(): YEES real cand found at %d ,%d in %d\n", ix, iy, real_cand_index); 
        //printf("ComputeMultiStreamLines:update(): sizes: %d %d \n", S->cand_coord.size(), S->cand_coord_cand.size()  ); 
                


        // decide in wich triangle the point is: 

tvector r = interval_xy[2] - interval_xy[0]; 
double k = r[1]/r[0]; 
double d = -k * interval_xy[0][0] + interval_xy[0][1]; 

std::vector<int> in_eras; // delete via reindexing, to make it a bit confuseing 
        in_eras.resize(4); 
        in_eras[0] = 0; 
        in_eras[1] = 1; 
        in_eras[2] = 2; 
        in_eras[3] = 3; 

        //S->points.push_back(point(CurrentPoint[0], CurrentPoint[1], interval_xy[0][2] )); 

        // is point over line? 
        if( CurrentPoint[1] > k * CurrentPoint[0] + d  )
        {
                if( interval_xy[1][1] > k *  interval_xy[1][0] + d )
                { // interval point 1 is also over -> erase point 3
                        
                        in_eras.erase(in_eras.begin() + 3); 
//                                      S->points.push_back( interval_xy[3] );
                } 
                else
                {
                        in_eras.erase(in_eras.begin() + 1); 
//                                      S->points.push_back( interval_xy[1] );
                }
        } 
        else // is under line
        {
                if( interval_xy[1][1] <= k *  interval_xy[1][0] + d )
                { // interval point 1 is also under -> erase point 3
                        in_eras.erase(in_eras.begin() + 3); 
//                           S->points.push_back( interval_xy[3] );
                } 
                else
                {
                        in_eras.erase(in_eras.begin() + 1); 
//                        S->points.push_back( interval_xy[1] );
                }        
        } 


        //S->lines.push_back( interval_xy[ in_eras[0] ] ); 
        //S->lines.push_back( interval_xy[ in_eras[1] ] ); 

        //S->lines.push_back( interval_xy[ in_eras[1] ] ); 
        //S->lines.push_back( interval_xy[ in_eras[2] ] ); 

        //S->lines.push_back( interval_xy[ in_eras[2] ] ); 
        //S->lines.push_back( interval_xy[ in_eras[0] ] ); 


        //in_eras.erase(in_eras.begin() + dist_max_index); 

        //printf("ineras: %d %d %d, size %d\n", in_eras[0], in_eras[1], in_eras[2], in_eras.size()); 

        // interpol via cutting lines in z plane 

point pv; 
        pv = CurrentPoint[0], CurrentPoint[1], interval_xy[ in_eras[0] ][2]; 

tvector dir = pv - interval_xy[ in_eras[0] ]; 

double lambda, mue; 
        intersect( interval_xy[ in_eras[1] ], interval_xy[ in_eras[2] ], CurrentPoint, dir, &lambda, &mue  ); 

        //printf("lambda %f, mue %f\n", lambda, mue); 

        if(debug == 1)
        {
                //S->points.push_back( interval_xy[ in_eras[0] ] +  mue * dir  ); 
                //S->points.push_back( interval_xy[ in_eras[0] ]); 
                //S->points.push_back( interval_xy[ in_eras[1] ] +  lambda * (interval_xy[ in_eras[2] ] - interval_xy[ in_eras[1] ])  ); 
        } 

tvector tmp_pos=   lambda*      tvector( interval_index_xy_perm[ in_eras[1] ][0], interval_index_xy_perm[ in_eras[1] ][1], interval_index_xy_perm[ in_eras[1] ][2] )
                   + (1 - lambda)* tvector( interval_index_xy_perm[ in_eras[2] ][0], interval_index_xy_perm[ in_eras[2] ][1], interval_index_xy_perm[ in_eras[2] ][2] ); 

point s = interval_xy[ in_eras[1] ] +  lambda * (interval_xy[ in_eras[2] ] - interval_xy[ in_eras[1] ]); 

double n1 = dist_xy( s, pv ); 
double n2 = dist_xy( pv, interval_xy[ in_eras[0] ]  ); 
double ng = n1 + n2; 

        n1 = n1 / ng; 

point float_index; 
        float_index =   n1 * tmp_pos + 
                        (1 - n1)* tvector( interval_index_xy_perm[ in_eras[0] ][0], interval_index_xy_perm[ in_eras[0] ][1], interval_index_xy_perm[ in_eras[0] ][2] ); 



// printf("lambda %f, mue %f, tmp_pos %f %f\n", lambda, mue, float_index[0], float_index[1]); 

MultiIndex<3> FieldSize = CoordsArr->Size(); 
MultiIndex<3> ifirst, ilast; 
ilast = FieldSize[0] - 1, FieldSize[1] - 1, FieldSize[2] - 1; 
ifirst = 0,0,0; 

point p1, p2, max_float_index; 
        p1 = (*CoordsArr)[ifirst]; 
        p2 = (*CoordsArr)[ilast]; 
        max_float_index = ilast[0], ilast[1], ilast[2]; 

// calc accurate z-float index (cheating by assuming kind of cylindrical coordinates with z as axis) 
double z_height = p2[2] - p1[2]; 

        // interpolated float index 
float_index[ Z[2] ] = (CurrentPoint[2] - p1[2]) / z_height * max_float_index[Z[2]] ; 


//printf("z_height: %f, Current_z - p1z: %f, max_ind: %f\n", z_height, (CurrentPoint[2] - p1[2]), max_float_index[Z[2]]); 
//printf("ComputeMultiStreamLines:update(): index: %f, %f, %f\n", float_index[0], float_index[1], float_index[2] ); 



MultiIndex<3> c0, c1, c2, c3, c4, c5, c6, c7; 
         c0 = (int)round(floor(float_index[0])), (int)round(floor(float_index[1])), (int)round(floor(float_index[2])); 


         for(int i = 0; i < (int)round(max_float_index[Z[2]]) - 1; i++ )
         {
         MultiIndex<3> da, db; 

                 da = c0[0], c0[1], c0[2]; 
                 db = c0[0], c0[1], c0[2]; 

                 da[Z[2]] = i; 
                 db[Z[2]] = i+1; 

         point pa, pb; 
                 pa = (*CoordsArr)[ da ]; 
                 pb = (*CoordsArr)[ db ]; 

                 //printf("da %d %d %d, db %d %d %d\n", da[0], da[1], da[2], db[0], db[1], db[2]); 
                 //printf("dists [%d]: norm %f, z-dist: %f\n", i, norm((pb-pa)), pb[2]-pa[2]); 

                 //printf("Its %f not here! betweeb %f and %f\n", CurrentPoint[2], pa[2], pb[2]); 

                 if( pa[2] <= CurrentPoint[2] && CurrentPoint[2] <= pb[2] )
                 {
                         //printf("Its heeeeeere! betweeb %f and %f\n", pa[2], pb[2]); 
                         break;
                 } 

                 if (pa[2] >= CurrentPoint[2] && CurrentPoint[2] >= pb[2]  )
                 {
                         //printf("Its heeeeeere! betweeb %f and %f\n", pa[2], pb[2]); 
                         
                         float_index[Z[2]] = 1 - (CurrentPoint[2] - pb[2]) / (pa[2] - pb[2])  +  i; 

                         //printf("zind: %f\n", float_index[Z[2]]); 

                         break;
                 } 
                 
         } 

         //printf("ComputeMultiStreamLines:update(): index: %f, %f, %f\n", float_index[0], float_index[1], float_index[2] ); 

         if(debug == 1)
         {
                 c0[Z[2]] = (int)round(floor(float_index[Z[2]])); 
                 //printf("c0: %d %d %d\n", c0[0], c0[1], c0[2]); 

                 c1 = c0[0],   c0[1]+1, c0[2]; 
                 c2 = c0[0]+1, c0[1]+1, c0[2]; 
                 c3 = c0[0]+1, c0[1],   c0[2]; 
                 c4 = c0[0],   c0[1],   c0[2]+1; 
                 c5 = c0[0],   c0[1]+1, c0[2]+1; 
                 c6 = c0[0]+1, c0[1]+1, c0[2]+1; 
                 c7 = c0[0]+1, c0[1],   c0[2]+1; 

                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c0], (*CoordsArr)[c1]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c1], (*CoordsArr)[c2]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c2], (*CoordsArr)[c3]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c3], (*CoordsArr)[c0]), YELLOW) ); 

                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c4], (*CoordsArr)[c5]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c5], (*CoordsArr)[c6]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c6], (*CoordsArr)[c7]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c7], (*CoordsArr)[c4]), YELLOW) ); 

                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c0], (*CoordsArr)[c4]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c1], (*CoordsArr)[c5]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c2], (*CoordsArr)[c6]), YELLOW) ); 
                 debuglines.push_back( colorLine_t( line_t( (*CoordsArr)[c3], (*CoordsArr)[c7]), YELLOW) ); 


         } 

         return float_index;
}


int ComputeMultiStreamLines::HexaHedralCellNewton::pointInCell(const point&P, const int debug, std::vector<colorLine_t>& debuglines)
{
        if(!UniMap)
        {
                printf("No unimap available");fflush(stdout);
                return -1;
        } 

        point_type = UniMap->localCellCoordinatesFromCurviGrid( P, uvw );

        return point_type;
}

point ComputeMultiStreamLines::HexaHedralCellNewton::localUVW   (const point&p, const int debug, std::vector<colorLine_t>& debuglines)
{
        return uvw;
}


unsigned int ComputeMultiStreamLines::BBSelectorMemoryTree::getNext()
{
        if(current == 0)
        {
//              printf("returned last, without tree\n");
                current++;
                return last; 
        } 
        if(current == 1) // query the tree here
        {
                //printf("querying creating tree\n");
                m.clear();
                S->Tree->nearest_range(reset_pos, S->max_radius * S->tree_fac, m); 
                //printf("Map Size %d\n", m.size()); 

                // for(map<double,int>::iterator i = m.begin(); i != m.end(); i++)
//                      printf("map: %f %d\n", i->first ,i->second);
                
                size = m.size()+1; 
                current++; 

                if(m.size() == 0) { return 0; } 

                it = m.begin(); 
         int tmp =  it->second; 
                 it++; 
                 return tmp;
        }
        else
        {
        int tmp = it->second; 
                it++; 
                current++;
                return tmp;
        }
}

void ComputeMultiStreamLines::BBSelectorMemoryTree::reset(const int _last, const point& pos)
{
        last = _last; 
        reset_pos = pos;

        current = 0;

        it = m.begin(); 
}

bool ComputeMultiStreamLines::BBSelectorMemoryTree::hasNext()
{
        if(current == 0) // theres always a last to test
                return true; 

        if(current == 1) // tree query has to be done
                return true;

        if( current < size)
                return true; 
        else
                return false;
}

unsigned int ComputeMultiStreamLines::BBSelectorMemoryAll::getNext()
{
        if( count % 2 == 0 )
        {
                if( (last + count_up) < size )
                {
                        current = last + count_up; 
                        count_up++;
                } 
                else
                {
                        current = last - count_down; 
                        count_down++;
                }
        } 
        else
        {
                if( (last - count_down) >= 0 )
                {
                        current = last - count_down; 
                        count_down++;
                } 
                else
                {
                        current = last + count_up; 
                        count_up++;
                }       
        } 

        count++; 
        return current;
}


bool ComputeMultiStreamLines::BBSelectorMemoryAll::hasNext()
{
        if( count < size )
                return true; 
        else
                return false;
}



static Ref<VCreator<ComputeMultiStreamLines, AcceptList<Fiber::Field> > > MyComputeMultiStreamLinesCreator("IntegralLines/ComputeMultiStreamLines");