59
65
POS_ERR_FINAL // Not actually an error, but tells us where the end of the num is.
62
#ifdef __LITTLE_ENDIAN__
64
void floatSwapBytes(float *inFloat);
68
//Dump point in form (x,y,z) to stderr
69
void DumpPoint(const Point3D &pt);
71
//!Make a pos file from a set of points, assigning them a mass of "mass"
72
void makePos(const vector<Point3D> &points, float mass, const char *name);
75
69
//!make a pos file from a set of a set of IonHits
76
void makePos(const vector<IonHit> &points, const char *name);
70
unsigned int IonVectorToPos(const vector<IonHit> &points, const std::string &name);
78
72
//!make/append to a pos file from a set of a set of IonHits
79
73
void appendPos(const vector<IonHit> &points, const char *name);
81
//!Print short summmary of the ions in a list
82
void PrintSummary(std::vector<IonHit> &posIons, RangeFile &rangeFile);
85
//!Returns the limits of a dataset of Ions
86
void dataLimits(const vector<IonHit> &posIons, Point3D &low, Point3D &upper);
91
//! Returns the shortest distance between a facet and a Point
92
/* The inputs are the facet points (ABC) and the point P.
93
* distance is shortest using standard plane version
94
* \f$ D = \vec{AB} \cdot \vec{n} \f$
95
* iff dot products to each combination of \f$ \left( AP,BP,CP \right) \leq 0 \f$
96
* otherwise closest point is on the boundary of the simplex.
97
* tested by shortest distance to each line segment (E is shortest pt. AB is line segement)
98
* \f$ \vec{E} = \frac{\vec{AB}}{|\vec{AB}|}
99
* ( \vec{PB} \cdot \vec{AB})\f$
101
float distanceToFacet(const Point3D &fA, const Point3D &fB,
102
const Point3D &fC, const Point3D &p, const Point3D &normal);
105
//! Returns the shortest distance between a line segment and a given point
106
/* The inpus are the ends of the line segment and the point. Uses the formula that
108
* D = \abs{\vec{PE}}\f$
110
* \mathrm{~if~} \vec{PA} \cdot \vec{AB} > 0
111
* \rightarrow \vec{PE} = \vec{A}
114
* \mathrm{~if~} \vec{AB} \cdot \vec{PB} > 0 ~\&~ \vec{PA} \cdot \vec{AB} < 0
115
* \rightarrow \vec{PB} \cdot \frac{\vec{AB}}{\abs{\vec{AB}}}
118
* \mathrm{~if~} \vec{PB} \cdot \vec{AB} < 0
119
* \rightarrow \vec{B}
122
float distanceToSegment(const Point3D &fA, const Point3D &fB, const Point3D &p);
125
//!Check which way vectors attached to two 3D points "point", - "together", "apart" or "in common"
127
unsigned int vectorPointDir(const Point3D &fA, const Point3D &fB,
128
const Point3D &vA, const Point3D &vB);
131
void DumpIonInfo(const std::vector<IonHit> &vec,
134
//!Return the minimum of 2 objects(inline)
135
template<class T> inline T min2(T T1, T T2)
145
//!Return the minimum of 3 objects (inline)
146
template<class T> inline T min3(T T1,T T2,T T3)
148
//check Tor T1 smallest
149
if(T1 < T2 && T1 < T3)
152
//either T2 or T3 is smallest
75
//!Set the bounds from an array of ion hits
76
BoundCube getIonDataLimits(const vector<IonHit> &p);//
159
78
//Range file strucutres
161
//!Data holder for colour as float
169
81
enum{ RANGE_ERR_OPEN =1,
385
278
* x,y,z,mass/charge.
387
280
//!Load a pos file into a T of IonHits
389
unsigned int GenericLoadFloatFile(int inputnumcols, int outputnumcols, int index[], T &posIons,const char *posFile, unsigned int &progress, bool (*callback)())
391
//buffersize must be a power of two and at least sizeof(IONHIT)
392
const unsigned int NUMROWS=512;
393
const unsigned int BUFFERSIZE=inputnumcols * sizeof(float) * NUMROWS;
394
const unsigned int BUFFERSIZE2=outputnumcols * sizeof(float) * NUMROWS;
395
char *buffer=new char[BUFFERSIZE];
396
char *buffer2=new char[BUFFERSIZE2];
399
return POS_ALLOC_FAIL;
402
std::ifstream CFile(posFile,std::ios::binary);
408
return POS_OPEN_FAIL;
411
CFile.seekg(0,std::ios::end);
412
size_t fileSize=CFile.tellg();
418
return POS_EMPTY_FAIL;
421
CFile.seekg(0,std::ios::beg);
423
//calculate the number of points stored in the POS file
428
size_t curBufferSize=BUFFERSIZE;
429
size_t curBufferSize2=BUFFERSIZE2;
431
if(fileSize % (inputnumcols * sizeof(float)))
435
return POS_SIZE_MODULUS_ERR;
440
posIons.resize(fileSize/sizeof(IONHIT));
442
catch(std::bad_alloc)
446
return POS_ALLOC_FAIL;
450
while(fileSize < curBufferSize) {
451
curBufferSize = curBufferSize >> 1;
452
curBufferSize2 = curBufferSize2 >> 1;
455
//Technically this is dependant upon the buffer size.
456
unsigned int curProg = 10000;
458
int maxCols = inputnumcols * sizeof(float);
459
int maxPosCols = outputnumcols * sizeof(float);
462
while((size_t)CFile.tellg() <= fileSize-curBufferSize)
464
CFile.read(buffer,curBufferSize);
469
return POS_READ_FAIL;
472
for (int j = 0; j < NUMROWS; j++) // iterate through rows
474
for (int i = 0; i < outputnumcols; i++) // iterate through floats
476
memcpy(&(buffer2[j * maxPosCols + i * sizeof(float)]),
477
&(buffer[j * maxCols + index[i] * sizeof(float)]), sizeof(float));
481
hitStruct = (IONHIT *)buffer2;
483
for(ui=0; ui<curBufferSize2; ui+=(sizeof(IONHIT)))
485
hit.setHit(hitStruct);
486
//Data bytes stored in pos files are big
487
//endian. flip as required
488
#ifdef __LITTLE_ENDIAN__
496
return POS_NAN_LOAD_ERROR;
507
progress= (unsigned int)((float)(CFile.tellg())/((float)fileSize)*100.0f);
508
curProg=PROGRESS_REDUCE;
514
return POS_ABORT_FAIL;
521
curBufferSize = curBufferSize >> 1 ;
522
curBufferSize2 = curBufferSize2 >> 1 ;
523
}while(curBufferSize2 >= sizeof(IONHIT));
525
ASSERT((unsigned int)CFile.tellg() == fileSize);
534
unsigned int LimitLoadPosFile(int inputnumcols, int outputnumcols, int index[], T &posIons,const char *posFile, size_t limitCount,
535
unsigned int &progress, bool (*callback)())
537
//buffersize must be a power of two and at least sizeof(IONHIT)
538
const unsigned int NUMROWS=1;
539
const unsigned int BUFFERSIZE=inputnumcols * sizeof(float) * NUMROWS;
540
const unsigned int BUFFERSIZE2=outputnumcols * sizeof(float) * NUMROWS;
541
char *buffer=new char[BUFFERSIZE];
542
char *buffer2=new char[BUFFERSIZE2];
545
return POS_ALLOC_FAIL;
548
std::ifstream CFile(posFile,std::ios::binary);
554
return POS_OPEN_FAIL;
557
CFile.seekg(0,std::ios::end);
558
size_t fileSize=CFile.tellg();
564
return POS_EMPTY_FAIL;
567
CFile.seekg(0,std::ios::beg);
569
//calculate the number of points stored in the POS file
572
size_t maxCols = inputnumcols * sizeof(float);
575
if(fileSize % maxCols)
579
return POS_SIZE_MODULUS_ERR;
582
maxIons =fileSize/maxCols;
583
limitCount=std::min(limitCount,maxIons);
585
//If we are going to load the whole file, don't use a sampling method to do it.
586
if(limitCount == maxIons)
592
//Try opening it using the normal functions
593
return GenericLoadFloatFile(inputnumcols, outputnumcols, index, posIons,posFile,progress, callback);
596
//Use a sampling method to load the pos file
597
std::vector<size_t> ionsToLoad;
600
posIons.resize(limitCount);
605
randomDigitSelection(ionsToLoad,maxIons,rng,
606
limitCount,dummy,callback);
608
catch(std::bad_alloc)
612
return POS_ALLOC_FAIL;
617
GreaterWithCallback<size_t> g(callback,PROGRESS_REDUCE);
618
std::sort(ionsToLoad.begin(),ionsToLoad.end(),g);
620
unsigned int curProg = PROGRESS_REDUCE;
622
//TODO: probably not very nice to the disk drive. would be better to
623
//scan ahead for contigous data regions, and load that where possible.
624
//Or switch between different algorithms based upon ionsToLoad.size()/
626
std::ios::pos_type nextIonPos;
627
for(size_t ui=0;ui<ionsToLoad.size(); ui++)
629
nextIonPos = ionsToLoad[ui]*maxCols;
631
if(CFile.tellg() !=nextIonPos )
632
CFile.seekg(nextIonPos);
634
CFile.read(buffer,BUFFERSIZE);
636
for (int i = 0; i < outputnumcols; i++) // iterate through floats
637
memcpy(&(buffer2[i * sizeof(float)]), &(buffer[index[i] * sizeof(float)]), sizeof(float));
639
memcpy((char *)(&hit), buffer2, sizeof(IONHIT));
642
return POS_READ_FAIL;
643
posIons[ui].setHit(&hit);
644
//Data bytes stored in pos files are big
645
//endian. flip as required
646
#ifdef __LITTLE_ENDIAN__
647
posIons[ui].switchEndian();
650
if(posIons[ui].hasNaN())
654
return POS_NAN_LOAD_ERROR;
661
progress= (unsigned int)((float)(CFile.tellg())/((float)fileSize)*100.0f);
662
curProg=PROGRESS_REDUCE;
667
return POS_ABORT_FAIL;
680
unsigned int LimitRestrictLoadPosFile(T &posIons,const char *posFile, size_t limitCount,
681
const BoundCube &bound, unsigned int &progress, bool (*callback)())
683
//buffersize must be a power of two and at least sizeof(IONHIT)
684
const unsigned int BUFFERSIZE=8192;
685
char *buffer=new char[BUFFERSIZE];
688
return POS_ALLOC_FAIL;
691
std::ifstream CFile(posFile,std::ios::binary);
696
return POS_OPEN_FAIL;
699
CFile.seekg(0,std::ios::end);
700
size_t fileSize=CFile.tellg();
705
return POS_EMPTY_FAIL;
708
CFile.seekg(0,std::ios::beg);
710
//calculate the number of points stored in the POS file
715
if(fileSize % sizeof(IONHIT))
718
return POS_SIZE_MODULUS_ERR;
721
maxIons =fileSize/sizeof(IONHIT);
722
limitCount=std::min(limitCount,maxIons);
724
//If we are going to load the whole file, don't use a sampling method to do it.
725
if(limitCount == maxIons)
730
//Try opening it using the normal functions
734
return GenericLoadFloatFile(4, 4, index, posIons,posFile,progress, callback);
737
//Use a sampling method to load the pos file
740
posIons.resize(limitCount);
742
catch(std::bad_alloc)
744
return POS_ALLOC_FAIL;
751
const size_t RECORD_WIDTH=16;
753
//FIXME: If limitCount is near maxIons, this will be
754
//very inefficient better to remember only the gaps in that case.
756
//OK, this algorithm is a bit tricky. We need to try to load
757
//the limit size, but discard ions that lie outside the
758
//desired bounding volume. We don't know in advance how many we will get
759
//until we inspect the data, so we have to guess that all of them fit, then
760
//keep trying more until we run out of ions.
762
//So keep adding new ions, limitCount at a time, checking that the new
763
//indicies are unique in the load.
765
//Random sequence to load for this pass and already loaded sequence respectively
766
std::vector<size_t> ionsToLoad,ionsLoaded;
769
//Step 1: Generate random sequence
772
//Create an array of random numbers to load
773
for(size_t ui=0; ui<limitCount; ui++)
774
ionsToLoad[ui] = (rng.genUniformDev()*(maxIons-1));
776
//Sort & Remove repeated indices
777
std::sort(ionsToLoad.begin(),ionsToLoad.end());
778
std::unique(ionsToLoad.begin(),ionsToLoad.end());
780
//Step 2: Remove any entries we have already loaded from the list
785
//Sort already loaded ions before doing binary search below
786
std::sort(ionsLoaded.begin(),ionsLoaded.end());
788
for(unsigned int ui=0;ui<ionsToLoad.size()-seenCount;ui++)
790
//Keep swapping out bad entries until we find a good one
791
while(std::binary_search(ionsLoaded.begin(),ionsLoaded.end(),
792
ionsToLoad[ui]) && ui!= ionsToLoad.size()-seenCount)
794
//Swap this "bad" ion out for a fresh one
795
std::swap(ionsToLoad[ui],ionsToLoad[ionsToLoad.size()-seenCount]);
801
//Sort the vector front again, as we have disrupted the order during discard
802
//(should be near sorted, so not too costly)
803
//Ignore the back indicies,a as we wont use them anymore.
804
//Being sorted means that we don't have to spin the disc backwards, and can jump forwards
806
std::sort(ionsToLoad.begin(),ionsToLoad.end()-seenCount);
810
unsigned int curProg = PROGRESS_REDUCE;
812
std::ios::pos_type nextIonPos;
813
unsigned int ionStart=0;
814
for(size_t ui=0;ui<ionsToLoad.size()-seenCount; ui++)
816
nextIonPos = ionsToLoad[ui]*sizeof(IONHIT);
817
ionsLoaded.push_back(ionsToLoad[ui]);
819
if(CFile.tellg() !=nextIonPos )
820
CFile.seekg(nextIonPos);
823
CFile.read((char *)&hit,sizeof(IONHIT));
825
return POS_READ_FAIL;
828
//Reject the ion if it is not in the
830
if(bound.containsPt(Point3D(hit.pos[0],hit.pos[1],hit.pos[2])))
832
posIons[pointCount].setHit(&hit);
833
//Data bytes stored in pos files are big
834
//endian. flip as required
835
#ifdef __LITTLE_ENDIAN__
836
posIons[pointCount].switchEndian();
838
//Check for nan. If we find one, abort
839
if(posIons[pointCount].hasNaN())
842
return POS_NAN_LOAD_ERROR;
848
//Update progress bar
852
progress= (unsigned int)((float)(pointCount)/((float)limitCount)*100.0f);
853
curProg=PROGRESS_REDUCE;
859
return POS_ABORT_FAIL;
865
//Keep loading until we can do no more, or we have enough.
866
}while(pointCount < limitCount && ionsLoaded.size() != fileSize/RECORD_WIDTH);
874
//!Load a pos file into a T of IonHits, excluding ions not within the bounds
876
unsigned int RestrictLoadPosFile(T &posIons,const char *posFile, const BoundCube &b,
877
unsigned int &progress, bool (*callback)())
881
//buffersize must be a power of two and at least sizeof(IONHIT)
882
const unsigned int BUFFERSIZE=8192;
883
char *buffer=new char[BUFFERSIZE];
886
return POS_ALLOC_FAIL;
889
std::ifstream CFile(posFile,std::ios::binary);
894
return POS_OPEN_FAIL;
897
CFile.seekg(0,std::ios::end);
898
size_t fileSize=CFile.tellg();
903
return POS_EMPTY_FAIL;
906
CFile.seekg(0,std::ios::beg);
908
//calculate the number of points stored in the POS file
913
size_t curBufferSize=BUFFERSIZE;
915
if(fileSize % sizeof(IONHIT))
918
return POS_SIZE_MODULUS_ERR;
923
posIons.resize(fileSize/sizeof(IONHIT));
925
catch(std::bad_alloc)
927
return POS_ALLOC_FAIL;
931
while(fileSize < curBufferSize)
932
curBufferSize = curBufferSize >> 1;
934
//Technically this is dependant upon the buffer size.
935
unsigned int curProg = PROGRESS_REDUCE;
939
while((size_t)CFile.tellg() <= fileSize-curBufferSize)
941
CFile.read(buffer,curBufferSize);
945
return POS_READ_FAIL;
948
hitStruct = (IONHIT *)buffer;
950
for(ui=0; ui<curBufferSize; ui+=(sizeof(IONHIT)))
952
hit.setHit(hitStruct);
953
//Data bytes stored in pos files are big
954
//endian. flip as required
955
#ifdef __LITTLE_ENDIAN__
962
return POS_NAN_LOAD_ERROR;
965
if(b.containsPt(hit.getPos()))
978
progress= (unsigned int)((float)(CFile.tellg())/((float)fileSize)*100.0f);
979
curProg=PROGRESS_REDUCE;
984
return POS_ABORT_FAIL;
991
curBufferSize = curBufferSize >> 1 ;
992
}while(curBufferSize >= sizeof(IONHIT));
994
ASSERT((unsigned int)CFile.tellg() == fileSize);
1001
//!Set the bounds from an array of ion hits
1002
BoundCube getIonDataLimits(const vector<IonHit> &p);//
281
unsigned int GenericLoadFloatFile(unsigned int inputnumcols, unsigned int outputnumcols,
282
unsigned int index[], vector<IonHit> &posIons,const char *posFile,
283
unsigned int &progress, bool (*callback)());
286
unsigned int LimitLoadPosFile(unsigned int inputnumcols, unsigned int outputnumcols, unsigned int index[],
287
vector<IonHit> &posIons,const char *posFile, size_t limitCount,
288
unsigned int &progress, bool (*callback)(),bool strongRandom);
added
removed
src/APTClasses.h
