Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
121 changes: 54 additions & 67 deletions src/swt.c
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,7 @@ int formatarrayI(ptrdiff_t another[]);
int getlineJH(char line[], FILE *fdata);
int Interpolate(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float depth[], float distance[], float WavCombScale[], ptrdiff_t WavCombLeft[], ptrdiff_t WavCombRight[]);
int PostProcess(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float depth[], float distance[], float WavCombScale[], ptrdiff_t WavCombLeft[], ptrdiff_t WavCombRight[]);
float quick_select(float arr[], int n);
int ReadData (char *argv[], float distance[], float depth[], float longitude[], float latitude[]);
int tempInterp(float distance[], float depth[], float tempFiltered[],int start,int end);
int wavCompute(float distance[], float depth[], float WavCoeff[], float scale, float wavWeighting[], ptrdiff_t numberlines);
int wavInterpii(float WavCoeff[], ptrdiff_t numberlines, float CoeffThreshold, float distance[], float scale);
int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float depth[], float distance[], float WavCombScale[], ptrdiff_t WavCombLeft[], ptrdiff_t WavCombRight[]);
Expand Down Expand Up @@ -93,7 +91,8 @@ int formatarrayI(ptrdiff_t another[])
int getlineJH(char s[], FILE *fdata)

{
int c, i;
int c;
ptrdiff_t i;

for (i=0; (c=getc(fdata))!=EOF && c!='\n'; ++i)
s[i] = c;
Expand All @@ -116,7 +115,7 @@ and creates a filtered value.*/
int Interpolate(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float depth[], float distance[], float WavCombScale[], ptrdiff_t WavCombLeft[], ptrdiff_t WavCombRight[])

{
int e,f,g; /*loop variable*/
ptrdiff_t e,f,g; /*loop variable*/
float gradient; /*used in the linear interpolation of the space under objects*/


Expand Down Expand Up @@ -159,20 +158,20 @@ int Interpolate(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float d
int PostProcess(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float depth[], float distance[], float WavCombScale[], ptrdiff_t WavCombLeft[], ptrdiff_t WavCombRight[])

{
int b,c,d,e,f,h,i,j,k,l,m; /*loop variables*/
ptrdiff_t b,c,d,e,f,h,i,j,k,l,m; /*loop variables*/
float gradient; /*used in the linear interpolation of the space under objects*/
float area, outline; /*properties of the feature*/
float extradist, extraarea;
float ratio; /*ratio: outline length / Area */
float depthInterpL, depthInterpR; /*interpolated depths used in calculations*/
float maxratio; /*'success parameters' driving the routine*/
int moveL, moveR; /*number of data points to the left to try to shift the left hand edge*/
ptrdiff_t moveL, moveR; /*number of data points to the left to try to shift the left hand edge*/
int loopleft, loopright; /*loop logical variable*/
int mainloop, changes; /*logic variable for the main loop and number of changes*/
float objToDo[MAXCOLUMNS]; /*Array 1 = to do, 0 = done*/
int loopobjects; /*logic in the loop controlling the object that is being done*/
float MaxCoeff;
int MaxCoeffLocn; /*location of max remaining coefficient*/
ptrdiff_t MaxCoeffLocn; /*location of max remaining coefficient*/
int ObjectsCount; /*How many object are there still to do*/
int criteriaFulfilled;
float widthFract, MoveMax, Move;
Expand All @@ -182,7 +181,7 @@ int PostProcess(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float d
float height, heightmax, depthInterp;
float mingradient, RHSgradient, LHSgradient; /*gradients of the extremities and their minimum allowed value*/
float summitToEdge;
int InsidePosnLHS, InsidePosnRHS;
ptrdiff_t InsidePosnLHS, InsidePosnRHS;
float outsideFract, TooFlatFract; /*user-set geometric parameters for the slope based restriction*/
float HeightDiffFract, ScaleMulLimit; /*more user set parameters for real data*/
int TooBig, TooDeep; /*variable associated with ScaleMulLimit*/
Expand Down Expand Up @@ -425,7 +424,7 @@ int PostProcess(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float d
if (WavCombLeft[f] - (moveL +1) < 0)
{
criteriaFulfilled = 0;
printf("FAIL 1: Tries to move outside data\n");
/*("FAIL 1: Tries to move outside data\n");*/
}
/*Criterion 2*/
/*Version from Book16p26*/
Expand Down Expand Up @@ -454,7 +453,7 @@ int PostProcess(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float d
/*is between the limits of a larger object - in terms of its original coefficient*/
if ((WavCombLeft[f]-(moveL+1)) < WavCombRight[c] && (WavCombLeft[f]-(moveL+1)) > WavCombLeft[c] && WavComb[c] > WavComb[f])
{
/*printf("FAIL 3: Tries to move inside other larger object*/
/*("FAIL 3: Tries to move inside other larger object*/
criteriaFulfilled = 0;
}
}
Expand Down Expand Up @@ -538,8 +537,8 @@ int PostProcess(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float d
InsidePosnRHS--;
}
RHSgradient = (depth[InsidePosnRHS]-depth[(WavCombRight[f] + moveR)])/(distance[(WavCombRight[f] + moveR)]-distance[InsidePosnRHS]);
printf("Summit location: %td\n", summitlocation[f]);
printf("WavCombRight[%d] %td InsidePosnRHS: %d RHSgradient %f \n",f,WavCombRight[f]+moveR,InsidePosnRHS, RHSgradient);
/*("Summit location: %td\n", summitlocation[f]);*/
/*("WavCombRight[%d] %td InsidePosnRHS: %d RHSgradient %f \n",f,WavCombRight[f]+moveR,InsidePosnRHS, RHSgradient);*/

/*Make sure that the object hasn't expanded too much*/
TooBig = 0;
Expand Down Expand Up @@ -854,7 +853,7 @@ int ReadData (char *argv[], float distance[], float depth[], float longitude[],

{
int len;
int countlines;
ptrdiff_t countlines;
char line[MAXLINE];
float column1, column2, column3, column4;

Expand All @@ -874,19 +873,18 @@ int ReadData (char *argv[], float distance[], float depth[], float longitude[],
/*Reads 2 columns distance along track and depth*/
sscanf(line, "%f %f %f %f", &column1, &column2, &column3, &column4);
/* x,z in meters*/
/*printf("In main while,%f %f\n", column1, column2);*/
distance[countlines] = column1;
depth[countlines] = column2;
/*printf("%f %f %d\n", distance[countlines], depth[countlines],countlines);*/
/*("%f %f %d\n", distance[countlines], depth[countlines],countlines);*/
longitude[countlines] = column3;
latitude[countlines] = column4;
/*Read these in, although they're not used anywhere*/
/*Saves removing them from a lot of places*/
countlines++;
}

/*printf("distance[0] %f\n", distance[0]);
printf("distance[MAXCOLUMNS -1] %f\n", distance[MAXCOLUMNS -1]);*/
/*"(distance[0] %f\n", distance[0]);
"(distance[MAXCOLUMNS -1] %f\n", distance[MAXCOLUMNS -1]);*/

fclose (fdata);

Expand Down Expand Up @@ -1004,8 +1002,7 @@ int wavelet(float distance[], float depth[], char *argv[], ptrdiff_t numberlines



/*Combine the locations of objects across the scale*/
printf("Combining information from all the scales........\n");
/*Combine the locations of objects across the scales*/
wavInterpCombii(argv, WavComb, numberlines, distance, WavCombScale, ScaleInteract);

/*The final product of wavInterpCombii - send to file*/
Expand All @@ -1020,7 +1017,7 @@ int wavelet(float distance[], float depth[], char *argv[], ptrdiff_t numberlines

/*Turn the coefficients into a regional by finding the objects and doing linear
interpolation between their limits*/
printf("\nCreating a regional filter from the objects at all scales\n");
/*("\nCreating a regional filter from the objects at all scales\n");*/
wavInterpComb(WavComb, numberlines, WavFilt, depth, distance, WavCombScale, WavCombRight, WavCombLeft);

/*Do the second stage of processing - finding the exact limits of the objects*/
Expand All @@ -1029,15 +1026,6 @@ int wavelet(float distance[], float depth[], char *argv[], ptrdiff_t numberlines
/*Re-interpolate underneath the objects*/
Interpolate(WavComb, numberlines, WavFilt, depth, distance, WavCombScale, WavCombRight, WavCombLeft);



/*Got the blocked objects, now do some post processing to make their limits astetically acceptable*/
/*Firstly sort out any times when the cf pp 178 Book 2*/
/*WavPostPro(distance, depth, WavFilt, WavComb, numberlines, CoeffThreshold, WavCombScale);*/
/*WavPostProii(distance, depth, WavFilt, WavComb, numberlines, CoeffThreshold, WavCombScale, ScaleInteract);*/
/*WavPostPro(distance, depth, WavFilt, WavComb, numberlines, CoeffThreshold, WavCombScale);*/

/*Do I need to use WavPostPro again to check that WavPostProii hasn't messed it up????????????*/
/*And send Filtered to file*/
for(l = 0; l < numberlines; l++)
fprintf(fwavFilt, "%f %f %f\n", distance[l], depth[l], WavFilt[l]);
Expand Down Expand Up @@ -1287,7 +1275,7 @@ WavComb*/
int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float depth[], float distance[], float WavCombScale[], ptrdiff_t WavCombLeft[], ptrdiff_t WavCombRight[])

{
int a,b,c,d,e,f,g,k,h,i,j; /*loop variable*/
ptrdiff_t a,b,c,d,e,f,g,k,h,i,j; /*loop variable*/
float gradient; /*used in the linear interpolation of the space under objects*/
float area, outline; /*properties of the feature*/
float extradist, extraarea;
Expand All @@ -1299,22 +1287,21 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
gradient = 0;
countobjects = 0;

printf("Finding the limits of the objects\n");
/*Find the limits of the objects - store the location at the limit of the feature*/
/*i.e. where height of the feature is zero*/
/*NOTE: there can be no overlap at this stage*/
for (k = 0; k < numberlines; k++)
{
/*printf("WavComb[%d] %f\n",k, WavComb[k]);*/
/*("WavComb[%d] %f\n",k, WavComb[k])*/
if (WavComb[k] > 0)
{
countobjects++;
/*Find the left hand side limit*/
for (a = k; a >=0 && distance[a] > distance[k] - (WavCombScale[k]/4.0); a--)
{
/*printf("Scale/4: %f\n",(WavCombScale[k]/4.0));*/
/*("Scale/4: %f\n",(WavCombScale[k]/4.0))*/
WavCombLeft[k] = a;
/*printf("Distance of left limit [%d]: %f\n",a,distance[a]);*/
/*("Distance of left limit [%d]: %f\n",a,distance[a])*/
}
/*If the lhs hasn't moved, move it by one as this must be OK if a coefficient can be calcuated*/
if (WavCombLeft[k] == k)
Expand All @@ -1324,23 +1311,23 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
/*Find the right hand side limit*/
for (b = k; b < numberlines && distance[b] < distance[k] + (WavCombScale[k]/4.0); b++)
{
/*printf("Scale/4: %f\n",(WavCombScale[k]/4.0));*/
/*("Scale/4: %f\n",(WavCombScale[k]/4.0))*/
WavCombRight[k] = b;
/*printf("Distance of right limit[%d]: %f\n",b, distance[b]);*/
/*("Distance of right limit[%d]: %f\n",b, distance[b])*/
}
/*If the rhs hasn't moved, move it by one as this must be OK if a coefficient can be calcuated*/
if (WavCombRight[k] == k)
{
WavCombRight[k]++;
}
/*printf("Left limit: %d, %f km\n", WavCombLeft[k], distance[WavCombLeft[k]]);*/
/*printf("Right limit: %d, %f km\n",WavCombRight[k], distance[WavCombRight[k]]);*/
/*("Left limit: %d, %f km\n", WavCombLeft[k], distance[WavCombLeft[k]])*/
/*("Right limit: %d, %f km\n",WavCombRight[k], distance[WavCombRight[k]])*/
}
}
printf("No. OBJECTS is %d \n\n", countobjects);
/*("No. OBJECTS is %d \n\n", countobjects)*/


printf("Double-check that there is no overlap\n");
/*("Double-check that there is no overlap\n")*/
/*For each of the non-zero coefficients indicating an object ........*/
for (c = 0; c < numberlines; c++)
{
Expand All @@ -1352,16 +1339,16 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
/*And flag if any other left hand limit is between the limits for this object*/
if (WavCombLeft[d] >= WavCombLeft[c] && WavCombLeft[d] <= WavCombRight[c] && WavComb[d] > 0 && d != c)
{
printf("WARNING - whilst centre points not in each other's ranges, ranges overlap\n");
printf("Left limit located at [%d] is inside limits of object at [%d]\n",d,c);
printf("%f (lhs of right object) %f (rhs of left object)\n", distance[WavCombLeft[d]], distance[WavCombRight[c]]);
/*("WARNING - whilst centre points not in each other's ranges, ranges overlap\n")*/
/*("Left limit located at [%d] is inside limits of object at [%d]\n",d,c)*/
/*("%f (lhs of right object) %f (rhs of left object)\n", distance[WavCombLeft[d]], distance[WavCombRight[c]])*/
}
/*And flag if any other right hand limit is between the limits for this object*/
if (WavCombRight[d] >= WavCombLeft[c] && WavCombRight[d] <= WavCombRight[c] && WavComb[d] > 0 && d != c)
{
printf("WARNING - whilst centre points not in each other's ranges, ranges overlap\n");
printf("Right limit located at [%d] is inside limits of object at [%d]\n",d,c);
printf("%f (rhs of left object) %f (lhs of right object)\n", distance[WavCombRight[d]], distance[WavCombLeft[c]]);
/*("WARNING - whilst centre points not in each other's ranges, ranges overlap\n")*/
/*("Right limit located at [%d] is inside limits of object at [%d]\n",d,c)*/
/*("%f (rhs of left object) %f (lhs of right object)\n", distance[WavCombRight[d]], distance[WavCombLeft[c]])*/
}
}

Expand All @@ -1370,13 +1357,13 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float



printf("\nDoing initial interpolation.......\n");
/*("\nDoing initial interpolation.......\n")*/

/*Start of by setting WavFilt[] to depth[]*/
for (e = 0; e < numberlines; e++)
{
WavFilt[e] = depth[e];
/*printf("WavFilt[%d] %f\n",e,WavFilt[e]);*/
/*("WavFilt[%d] %f\n",e,WavFilt[e]);*/
}

/*In the objects, interpolate between the limits*/
Expand All @@ -1385,10 +1372,10 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
if (WavComb[f] > 0)
{
/*Gradient for the line under the object*/
/*printf("Extrapolate under object centred on [%d]\n",f);*/
/*printf("Limits are (%f,%f) (%f,%f)\n",depth[WavCombLeft[f]],distance[WavCombLeft[f]],depth[WavCombRight[f]],distance[WavCombRight[f]]);*/
/*("Extrapolate under object centred on [%d]\n",f);*/
/*("Limits are (%f,%f) (%f,%f)\n",depth[WavCombLeft[f]],distance[WavCombLeft[f]],depth[WavCombRight[f]],distance[WavCombRight[f]]);*/
gradient = (depth[WavCombLeft[f]]-depth[WavCombRight[f]])/(distance[WavCombLeft[f]]-distance[WavCombRight[f]]);
/*printf("gradient is %f\n",gradient);*/
/*("gradient is %f\n",gradient);*/
/*And do the interpolation*/
for (g = WavCombLeft[f]; g <= WavCombRight[f]; g++)
{
Expand All @@ -1397,7 +1384,7 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
}
}

printf("Calculating properties of the object\n");
/*("Calculating properties of the object\n");*/
/*In the objects, interpolate between the limits*/
for (h = 0; h < numberlines; h++)
{
Expand All @@ -1412,11 +1399,11 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
/*over the top - use pythag*/
extradist = sqrt(((depth[i]-depth[i+1])*(depth[i]-depth[i+1]))+((distance[i]-distance[i+1])*(distance[i]-distance[i+1])));
outline = outline + extradist;
/*printf("[%d] outline %f extradist %f\n",i, outline, extradist);*/
/*("[%d] outline %f extradist %f\n",i, outline, extradist);*/
}
/*and the line underneath - use pythag*/
outline = outline + sqrt(((depth[WavCombLeft[h]]-depth[WavCombRight[h]])*(depth[WavCombLeft[h]]-depth[WavCombRight[h]]))+((distance[WavCombLeft[h]]-distance[WavCombRight[h]])*(distance[WavCombLeft[h]]-distance[WavCombRight[h]])));
printf("[%d] total outline (inc. base) %f \n",i, outline);
/*("[%d] total outline (inc. base) %f \n",i, outline);*/
/*The area of the feature*/
/*There is no exception where the measured line crosses the extrapolated*/
for (j = WavCombLeft[h];j < WavCombRight[h]; j++ )
Expand All @@ -1425,11 +1412,11 @@ int wavInterpComb(float WavComb[], ptrdiff_t numberlines, float WavFilt[], float
depthInterpR = depth[WavCombLeft[h]] + ((distance[j+1]-distance[WavCombLeft[h]])*gradient);
extraarea = 0.5*((depth[j]-depthInterpL) + (depth [j+1] - depthInterpR))*(distance[j+1]-distance[j]);
area = area + extraarea;
/*printf("[%d] area %f extraarea %f\n",j, area, extraarea);*/
/*("[%d] area %f extraarea %f\n",j, area, extraarea);*/
}
printf("[%d] total area %f \n",j, area);
/*printf("[%d] total area %f \n",j, area);*/
ratio = area/outline;
printf("area/outline %f\n",ratio);
/*printf("area/outline %f\n",ratio);*/
}
}

Expand All @@ -1443,7 +1430,7 @@ int wavInterpCombii(char *argv[], float WavComb[], ptrdiff_t numberlines, float

{
int len;
int linenum;
ptrdiff_t linenum;
char line[MAXLINE];
float column1, column2, column3;
FILE *fdata; /*Interpreted Coefficients file - objects as selected max coeff, all scales sequentially*/
Expand Down Expand Up @@ -1504,7 +1491,7 @@ int wavInterpCombii_inner(int linenum, float WavComb[], float distance[], float
/*Go in reverse order through the array so that i'm taking the largest scale first*/
for (a = linenum -1; a >= 0; a--)
{
/*printf("Coeff to check: %f %f %f\n", distance[Arr1[a]], Arr2[a], Arr3[a]);*/
/*("Coeff to check: %f %f %f\n", distance[Arr1[a]], Arr2[a], Arr3[a]);*/
/*For each coefficient [a] remove it if it has a smaller coefficient than anything within its range*/
for (b = 0; b < linenum && Arr3[a] != 0; b++)
{
Expand All @@ -1513,14 +1500,14 @@ int wavInterpCombii_inner(int linenum, float WavComb[], float distance[], float
if (distdiff < Arr2[a]/4.0)
{
Rlogic = 1;
/*printf("Within range\n");*/
/*("Within range\n");*/
}
/*Now see if [a] should be eliminated by .....*/
/*1) Something thinner having a larger coefficient*/
if (Arr3[a] < Arr3[b] && Rlogic == 1 && a != b)
{
/*printf("%f %f %f vs %f %f %f\n", distance[Arr1[a]], Arr2[a], Arr3[a], distance[Arr1[b]], Arr2[b], Arr3[b]);*/
/*printf("coeff[a] %f is a smaller coeff than [b] %f \n\n", Arr3[a], Arr3[b]);*/
/*("%f %f %f vs %f %f %f\n", distance[Arr1[a]], Arr2[a], Arr3[a], distance[Arr1[b]], Arr2[b], Arr3[b]);*/
/*("coeff[a] %f is a smaller coeff than [b] %f \n\n", Arr3[a], Arr3[b]);*/
Arr3[a] = 0;
}
/*2) Something thinner being a decent size (e.g. 1/4 of coeff of larger scale) and substantially pointier*/
Expand All @@ -1532,8 +1519,8 @@ int wavInterpCombii_inner(int linenum, float WavComb[], float distance[], float
if (Arr3[a] != 0 && Arr3[a] < ((Arr3[b])*FlatnessCoeffFactor) && Rlogic == 1 && a != b && PairFlatnessFactor > FlatnessFactor)
{
Arr3[a] = 0;
/*printf("[a] eliminated by small and pointy\n");
printf("PairFlatnessFactor %f\n", PairFlatnessFactor);*/
/*("[a] eliminated by small and pointy\n");
("PairFlatnessFactor %f\n", PairFlatnessFactor);*/
}
}

Expand All @@ -1542,15 +1529,15 @@ int wavInterpCombii_inner(int linenum, float WavComb[], float distance[], float
/*If [a] is not to be eliminated, do some elimination*/
if (Arr3[a] > 0)
{
/*printf("NOT ELIMINATED %f %f %f\n", distance[Arr1[a]], Arr2[a], Arr3[a]);*/
/*("NOT ELIMINATED %f %f %f\n", distance[Arr1[a]], Arr2[a], Arr3[a]);*/
for (c = 0; c < linenum; c++)
{
Rlogic = 0;
distdiff = sqrt((distance[Arr1[a]] - distance[Arr1[c]])*(distance[Arr1[a]] - distance[Arr1[c]]));
if (distdiff < Arr2[a]/4.0)
{
Rlogic = 1;
/*printf("In Range: %f %f %f\n",distance[Arr1[c]], Arr2[c], Arr3[c]);*/
/*("In Range: %f %f %f\n",distance[Arr1[c]], Arr2[c], Arr3[c]);*/
}
/*Now, given the previous logic [c] should be eliminated if*/
/*in range with a smaller coefficient*/
Expand All @@ -1559,7 +1546,7 @@ int wavInterpCombii_inner(int linenum, float WavComb[], float distance[], float
if (Arr3[a] > Arr3[c] && a != c && Rlogic == 1 && Arr3[c] > 0)
{
Arr3[c] = 0;
/*printf("And eliminate\n");*/
/*("And eliminate\n");*/
}
}
}
Expand Down
Loading