973 lines
41 KiB
Python
973 lines
41 KiB
Python
#!/usr/bin/env python
|
|
|
|
# Surface Tension Calculation via Pendant Drop Method (Bashforth-Adams)
|
|
# Created by : Ren Tristan A. de la Cruz
|
|
# Created on : 2016 July 26
|
|
# Last Modified : 2025 January 27
|
|
# Contact : rentristandelacruz@gmail.com
|
|
# Contact : radelacruz@up.edu.ph
|
|
|
|
# ========== IMPORTS ==========
|
|
|
|
import math
|
|
import copy
|
|
import time
|
|
import sys
|
|
import cv2
|
|
import multiprocessing
|
|
import numpy as np
|
|
import scipy.optimize #as spo
|
|
|
|
# ========== READ THE IMAGE & EXTRACT THE EDGE COORDINATES ==========
|
|
|
|
def ExtractEdgeCoordinate(_ImageFilePath, _OutputFileName, _SyringeWidth):
|
|
|
|
# Detect the edges of the image.
|
|
# All edge pixels will white(255) and rest of the pixels will be black(0).
|
|
Image = cv2.imread(_ImageFilePath, cv2.IMREAD_COLOR)
|
|
EdgeImage = cv2.Canny(Image, 100, 200)
|
|
|
|
# Find the first edge point (white pixel).
|
|
# It is assumed that first edge point part of the syringe.
|
|
NumOfRow, NumOfCol = EdgeImage.shape
|
|
EdgePoint = (0, 0)
|
|
EdgePointFlag = 0
|
|
for RowIndex in range(NumOfRow):
|
|
for ColIndex in range(NumOfCol):
|
|
if EdgeImage[RowIndex, ColIndex] == 255:
|
|
EdgePoint = (RowIndex, ColIndex)
|
|
EdgePointFlag = 1
|
|
break
|
|
if EdgePointFlag == 1:
|
|
break
|
|
|
|
|
|
# Color the all the syringe and droplet pixels differently (70).
|
|
# Determine the left-most, right-most, and bottom droplet pixels.
|
|
EdgeImage[EdgePoint] = 70
|
|
DropletEdgePoints = []
|
|
DropletEdgePoints.append(EdgePoint)
|
|
LeftMostCol = 9999999999
|
|
RightMostCol = 0
|
|
BottomMostRow = 0
|
|
BottomMostRow2 = 0
|
|
CenterCol = 0
|
|
|
|
while (len(DropletEdgePoints) != 0):
|
|
CurrentPoint = DropletEdgePoints[0]
|
|
DropletEdgePoints.remove(CurrentPoint)
|
|
NeighborPoints = []
|
|
for RowOffset in range(-3,4, 1):
|
|
for ColOffset in range(-3,4,1):
|
|
if(RowOffset != 0 or ColOffset != 0):
|
|
NeighborPoints.append((CurrentPoint[0]+RowOffset, CurrentPoint[1]+ColOffset))
|
|
for PointRow, PointCol in NeighborPoints:
|
|
if (0 <= PointRow and PointRow < NumOfRow) and (0 <= PointCol and PointCol < NumOfCol):
|
|
if EdgeImage[PointRow][PointCol] == 255:
|
|
DropletEdgePoints.append((PointRow, PointCol))
|
|
EdgeImage[PointRow, PointCol] = 70
|
|
if PointRow > BottomMostRow:
|
|
BottomMostRow = PointRow
|
|
if PointCol > RightMostCol:
|
|
RightMostCol = PointCol
|
|
if PointCol < LeftMostCol:
|
|
LeftMostCol = PointCol
|
|
|
|
CenterCol = (LeftMostCol + RightMostCol)//2
|
|
for Row in range(BottomMostRow + 1):
|
|
if EdgeImage[Row][CenterCol] == 70:
|
|
BottomMostRow2 = Row
|
|
BottomMostRow = BottomMostRow2
|
|
|
|
|
|
# Determine the transition row between the syringe pixels and the droplet pixels.
|
|
PreviousWidth = 0
|
|
WidthIncCount = 5
|
|
WidthHistory = [0]*(WidthIncCount+1)
|
|
TransitionRow = 0
|
|
|
|
for Row in range(BottomMostRow):
|
|
LeftEdgeCol = 0
|
|
RightEdgeCol = 9999999
|
|
for Col in range(LeftMostCol-1, RightMostCol+1, 1):
|
|
if EdgeImage[Row,Col] == 70:
|
|
LeftEdgeCol = Col
|
|
break
|
|
for Col in range(RightMostCol+1, LeftMostCol-1, -1):
|
|
if EdgeImage[Row,Col] == 70:
|
|
RightEdgeCol = Col
|
|
break
|
|
CurrentWidth = RightEdgeCol - LeftEdgeCol
|
|
if (CurrentWidth != PreviousWidth) and (0 < CurrentWidth) and (CurrentWidth < (9999999-NumOfCol)):
|
|
for i in range(WidthIncCount):
|
|
WidthHistory[i] = WidthHistory[i+1]
|
|
WidthHistory[WidthIncCount] = CurrentWidth
|
|
WidthIncreaseCount = 0
|
|
for i in range(WidthIncCount):
|
|
if WidthHistory[i] < WidthHistory[i+1]:
|
|
WidthIncreaseCount += 1
|
|
if WidthIncreaseCount == WidthIncCount:
|
|
TransitionRow = Row
|
|
break
|
|
if (0 < CurrentWidth) and (CurrentWidth < (9999999-NumOfCol)):
|
|
PreviousWidth = CurrentWidth
|
|
|
|
|
|
# Determine the width of the syringe in terms of pixel count
|
|
WidthSum = 0
|
|
WidthCount = 0
|
|
SyringeLevel = int(0.75 * TransitionRow)
|
|
for Row in range(TransitionRow):
|
|
RightEdgeCol = 9999999
|
|
LeftEdgeCol = 0
|
|
for Col in range(LeftMostCol, CenterCol + 1, 1):
|
|
if EdgeImage[Row, Col] == 70:
|
|
EdgeImage[Row, Col] = 255
|
|
LeftEdgeCol = Col
|
|
break
|
|
for Col in range(RightMostCol, CenterCol - 1, -1):
|
|
if EdgeImage[Row, Col] == 70:
|
|
EdgeImage[Row, Col] = 255
|
|
RightEdgeCol = Col
|
|
break
|
|
CurrentWidth = RightEdgeCol - LeftEdgeCol
|
|
if WidthCount <= SyringeLevel:
|
|
if (0 < CurrentWidth) and (CurrentWidth < (9999999-NumOfCol)):
|
|
WidthSum += CurrentWidth
|
|
WidthCount += 1
|
|
AverageWidth = float(WidthSum)/float(WidthCount)
|
|
|
|
|
|
# Extract the pixel coordinates of all droplet pixels (color: 70)
|
|
PositiveEdgePoints = []
|
|
NegativeEdgePoints = []
|
|
for Row in range(BottomMostRow, TransitionRow-1 ,-1):
|
|
for Col in range(LeftMostCol, CenterCol + 1, 1):
|
|
if EdgeImage[Row, Col] == 70:
|
|
NegativeEdgePoints.append((Row, Col))
|
|
break
|
|
for Col in range(RightMostCol, CenterCol, -1):
|
|
if EdgeImage[Row, Col] == 70:
|
|
PositiveEdgePoints.append((Row, Col))
|
|
break
|
|
|
|
|
|
# Scale the point coordinates by using the syringe length
|
|
SyringeWidth = _SyringeWidth
|
|
ScalingFactor = SyringeWidth/AverageWidth
|
|
ScaledNegativeEdgePoints = []
|
|
NegativeCoordinatesFile = open(_OutputFileName+"_Negative.dat", 'w')
|
|
for Row, Col in NegativeEdgePoints:
|
|
Row2 = (ScalingFactor*float(BottomMostRow-Row)) + 0.0001
|
|
Col2 = ScalingFactor*float(Col - CenterCol)
|
|
NegativeCoordinatesFile.write(str(Col2) + " " + str(Row2) + "\n")
|
|
ScaledNegativeEdgePoints.append((Col2,Row2))
|
|
NegativeCoordinatesFile.close()
|
|
|
|
ScaledPositiveEdgePoints = []
|
|
PositiveCoordinatesFile = open(_OutputFileName+"_Positive.dat", 'w')
|
|
for Row, Col in PositiveEdgePoints:
|
|
Row2 = (ScalingFactor*float(BottomMostRow-Row)) + 0.0001
|
|
Col2 = ScalingFactor*float(Col - CenterCol)
|
|
PositiveCoordinatesFile.write(str(Col2) + " " + str(Row2) + "\n")
|
|
ScaledPositiveEdgePoints.append((Col2, Row2))
|
|
PositiveCoordinatesFile.close()
|
|
|
|
|
|
# Draw the detected edge points of the droplet.
|
|
Image2 = copy.deepcopy(Image)
|
|
for Row in range(0, NumOfRow):
|
|
for Col in range(0, NumOfCol):
|
|
Image[(Row,Col)] = 0
|
|
EdgeImage[(Row,Col)] = 0
|
|
for Point in PositiveEdgePoints:
|
|
Image[Point] = (0, 0, 255)
|
|
EdgeImage[Point] = 255
|
|
for Point in NegativeEdgePoints:
|
|
Image[Point] = (0, 0, 255)
|
|
EdgeImage[Point] = 255
|
|
|
|
# Draw the center column of the drop
|
|
cv2.line(Image, (CenterCol, TransitionRow), (CenterCol, BottomMostRow), (0, 255, 0), 1) # image, start, end, color(BGR), thickness
|
|
|
|
MaxWidthRatio = 0.0
|
|
for Row in range(BottomMostRow, TransitionRow-1 ,-1):
|
|
RightCol = 99999
|
|
LeftCol = -1
|
|
for Col in range(LeftMostCol, CenterCol + 1, 1):
|
|
if EdgeImage[Row, Col] == 255:
|
|
LeftCol = Col
|
|
break
|
|
for Col in range(RightMostCol, CenterCol, -1):
|
|
if EdgeImage[Row, Col] == 255:
|
|
RightCol = Col
|
|
break
|
|
WidthRatio = float(RightCol - LeftCol)/float(RightMostCol-LeftMostCol)
|
|
if WidthRatio < 1.0 and WidthRatio > MaxWidthRatio:
|
|
MaxWidthRatio = WidthRatio
|
|
|
|
RowWithMaxWidth = []
|
|
for Row in range(BottomMostRow, TransitionRow-1 ,-1):
|
|
RightCol = 99999
|
|
LeftCol = -1
|
|
for Col in range(LeftMostCol, CenterCol + 1, 1):
|
|
if EdgeImage[Row, Col] == 255:
|
|
LeftCol = Col
|
|
break
|
|
for Col in range(RightMostCol, CenterCol, -1):
|
|
if EdgeImage[Row, Col] == 255:
|
|
RightCol = Col
|
|
break
|
|
WidthRatio = float(RightCol - LeftCol)/float(RightMostCol-LeftMostCol)
|
|
if WidthRatio == MaxWidthRatio:
|
|
RowWithMaxWidth.append((Row, LeftCol, RightCol))
|
|
|
|
# Draw the main diameter
|
|
n = len(RowWithMaxWidth)
|
|
n = int(n/2)
|
|
Row = RowWithMaxWidth[n][0]
|
|
LeftCol = RowWithMaxWidth[n][1]
|
|
RightCol = RowWithMaxWidth[n][2]
|
|
cv2.line(Image, (LeftCol, Row), (RightCol, Row), (0, 255, 0), 1)
|
|
|
|
# Measure and draw the ten diameters
|
|
Diameters = []
|
|
MaxWidth = RightMostCol - LeftMostCol
|
|
for i in range(1, 13, 1):
|
|
for Col in range(LeftMostCol, CenterCol + 1, 1):
|
|
if EdgeImage[BottomMostRow-int(i*MaxWidth/10.0), Col] == 255:
|
|
LeftCol = Col
|
|
break
|
|
for Col in range(RightMostCol, CenterCol, -1):
|
|
if EdgeImage[BottomMostRow-int(i*MaxWidth/10.0), Col] == 255:
|
|
RightCol = Col
|
|
break
|
|
Diameters.append(RightCol-LeftCol)
|
|
cv2.line(Image, (LeftCol, BottomMostRow-int(i*MaxWidth/10.0)), (RightCol, BottomMostRow-int(i*MaxWidth/10.0)), (255, 0, 0), 1)
|
|
|
|
#DELETE THIS
|
|
#cv2.imwrite("1.jpg", Image)
|
|
#sys.exit()
|
|
|
|
|
|
# Scale the 12 drop diameters
|
|
for i in range(12):
|
|
Diameters[i] *= ScalingFactor
|
|
|
|
ScaledMaxWidth = ScalingFactor*float(MaxWidth)
|
|
ScaledUpWidth = ScalingFactor*float(RightCol - LeftCol)
|
|
return ScaledPositiveEdgePoints, ScaledNegativeEdgePoints, ScaledMaxWidth, Diameters
|
|
|
|
|
|
|
|
# ========== RK4 for BASHFORTH-ADAMS FORMULATION ==========
|
|
|
|
def BashforthAdamsRK4(_afRK4Table, _fApexCurvature, _fShapeFactor, _fSign):
|
|
A = _fApexCurvature
|
|
A2 = 2.0/A
|
|
B = _fShapeFactor
|
|
KX0 = 1.0
|
|
KT0 = 1.0
|
|
Inv6 = 1.0/6.0
|
|
fTotalDiffSquared = 0.0
|
|
n = 1
|
|
while n < len(_afRK4Table):
|
|
x = _afRK4Table[n-1][3] # predicted x
|
|
y = _afRK4Table[n-1][1] # actual y
|
|
t = _afRK4Table[n-1][4] # predicted angle
|
|
dy = _afRK4Table[n][2] # actual y diff
|
|
AG = 1.0 #AG: Added
|
|
ydy = y + 0.5 * dy
|
|
KX1 = dy * (1.0 / math.tan(t))
|
|
KT1 = dy * (1.0 / math.sin(t)) * (A2 - (B * y) - (math.sin(t) / (AG * x))) # A -> AG
|
|
KX2 = dy * (1.0 / math.tan(t + 0.5 * KT1))
|
|
STKT1 = math.sin(t + 0.5 * KT1)
|
|
KT2 = dy * (1.0 / STKT1) * (A2 - (B * ydy) - (STKT1 / (AG * (x + 0.5 * KX1)))) # A -> AG
|
|
KX3 = dy * (1.0 / math.tan(t + 0.5 * KT2))
|
|
STKT2 = math.sin(t + 0.5 * KT2)
|
|
KT3 = dy * (1.0 / STKT2) * (A2 - (B * ydy) - (STKT2 / (AG * (x + 0.5 * KX2)))) # A -> AG
|
|
KX4 = dy * (1.0 / math.tan(t + KT3))
|
|
STKT3 = math.sin(t + KT3)
|
|
KT4 = dy * (1.0 / STKT3) * (A2 - (B * (y + dy)) - (STKT3 / (AG * (x + KX3)))) # A -> AG
|
|
KX = (KX1 + 2.0 * (KX2 + KX3) + KX4) * Inv6
|
|
KT = (KT1 + 2.0 * (KT2 + KT3) + KT4) * Inv6
|
|
_afRK4Table[n][3] = x + KX
|
|
_afRK4Table[n][4] = t + KT
|
|
fTotalDiffSquared += (_afRK4Table[n][0] - (_fSign * _afRK4Table[n][3])) ** 2
|
|
n += 1
|
|
return fTotalDiffSquared
|
|
|
|
def MakeRK4table(_edgePoints):
|
|
InitialValues = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0]
|
|
|
|
rk4table = []
|
|
rk4table.append(InitialValues)
|
|
yPrev = 0.0
|
|
for x, y in _edgePoints:
|
|
EntryValues = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
|
|
EntryValues[0] = x
|
|
EntryValues[1] = y
|
|
EntryValues[2] = y - yPrev
|
|
yPrev = y
|
|
rk4table.append(EntryValues)
|
|
|
|
return rk4table
|
|
"""
|
|
RK4NegTable = []
|
|
RK4NegTable.append(InitialValues)
|
|
yPrev = 0.0
|
|
for x, y in _NegEdgePoints:
|
|
EntryValues = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
|
|
EntryValues[0] = x
|
|
EntryValues[1] = y
|
|
EntryValues[2] = y - yPrev
|
|
yPrev = y
|
|
RK4NegTable.append(EntryValues)
|
|
|
|
return RK4NegTable
|
|
"""
|
|
|
|
# ========== FINDING OPTIMAL PARAMETERS for BASHFORTH-ADAMS SOLUTION ==========
|
|
def FindOptimalParameters(_PosEdgePoints, _NegEdgePoints, _ARange_radCurv, _BRange_shapeFact, _XRange, _TRange_angle, _Params, _Ticks):
|
|
|
|
RK4PosTable = MakeRK4table(_PosEdgePoints)
|
|
RK4NegTable = MakeRK4table(_NegEdgePoints)
|
|
|
|
MinimumRMSD = 999999.0
|
|
BestParameters = [0, 0, 0, 0]
|
|
for A_radCurv in _ARange_radCurv:
|
|
for B_shapeFact in _BRange_shapeFact:
|
|
for X in _XRange:
|
|
for T_angle in _TRange_angle:
|
|
# x, y, dy, x*, t*, KX1, KT1, KX2, KT2, KX3, KT3, KX4, KT4
|
|
fTotalDiffSquared = 0.0
|
|
RK4PosTable[0][3] = X
|
|
RK4PosTable[0][4] = T_angle
|
|
RK4NegTable[0][3] = X
|
|
RK4NegTable[0][4] = T_angle
|
|
fTotalDiffSquared += BashforthAdamsRK4(RK4PosTable, A_radCurv, B_shapeFact, 1.0)
|
|
fTotalDiffSquared += BashforthAdamsRK4(RK4NegTable, A_radCurv, B_shapeFact, -1.0)
|
|
|
|
CurrentRMSD = (fTotalDiffSquared/(len(_PosEdgePoints) + len(_NegEdgePoints)))**0.5
|
|
if CurrentRMSD < MinimumRMSD:
|
|
MinimumRMSD = CurrentRMSD
|
|
BestParameters = [A_radCurv, B_shapeFact, X, T_angle]
|
|
_Ticks.put("#")
|
|
|
|
A_radCurv = BestParameters[0]
|
|
B_shapeFact = BestParameters[1]
|
|
X = BestParameters[2]
|
|
T_angle = BestParameters[3]
|
|
ST = (9.8 * A_radCurv) / B_shapeFact
|
|
_Params.put((MinimumRMSD, A_radCurv, B_shapeFact, X, T_angle))
|
|
|
|
|
|
# ========== FINDING OPTIMAL PARAMETERS USING MINIMZATION ==========
|
|
def OptimizeParameters(_PosEdgePoints, _NegEdgePoints, _AGuess_radCurv, _BGuess_shapeFact, _XGuess, _TGuess_angle):
|
|
|
|
RK4PosTable = MakeRK4table(_PosEdgePoints)
|
|
RK4NegTable = MakeRK4table(_NegEdgePoints)
|
|
"""
|
|
ARG_RAD = 0
|
|
ARG_SHAPE = 1
|
|
ARG_X = 2
|
|
ARG_ANGLE = 3
|
|
"""
|
|
|
|
def msd(args):
|
|
A_radCurv = args[0]
|
|
B_shapeFact = args[1]
|
|
X = args[2]
|
|
T_angle = args[3]
|
|
|
|
RK4PosTable[0][3] = X
|
|
RK4PosTable[0][4] = T_angle
|
|
RK4NegTable[0][3] = X
|
|
RK4NegTable[0][4] = T_angle
|
|
fTotalDiffSquared = BashforthAdamsRK4(RK4PosTable, A_radCurv, B_shapeFact, 1.0) + BashforthAdamsRK4(RK4NegTable, A_radCurv, B_shapeFact, -1.0)
|
|
|
|
return (fTotalDiffSquared/(len(_PosEdgePoints) + len(_NegEdgePoints)))**0.5
|
|
|
|
|
|
optResult = scipy.optimize.minimize(msd, [_AGuess_radCurv, _BGuess_shapeFact, _XGuess, _TGuess_angle])
|
|
if optResult.success:
|
|
print("minimization success")
|
|
#print(optResult.x)
|
|
#print(optResult.fun)
|
|
return [float(optResult.fun)] + optResult.x.tolist()
|
|
else:
|
|
return [9999999,0,0,0,0]
|
|
|
|
|
|
# ========== Juza's equation for 1/H ==========
|
|
|
|
def HJuza(_S_diamRatio, _K):
|
|
# LowerS, UpperS, Ke, Kex, Kep, K3, K2, K1, K0
|
|
|
|
# 0.82 - 0.999
|
|
Coefficient8 = [[ 0.82, 0.83, 1791452.8, 533.85850, 170.55837, -81676.15812, 202187.89690, -166836.27680, 45888.22973],
|
|
[ 0.83, 0.84, 189.07578, 272.52455, 72.730022, -19963.60152, 50019.79936, -41775.38186, 11629.85610],
|
|
[ 0.84, 0.85, 7.1549709, 165.45945, 35.285687, -6944.66848, 17609.43832, -14883.84379, 4193.34232],
|
|
[ 0.85, 0.87, 1.1496269, 95.066407, 12.600013, -2158.91585, 5571.60176, -4792.82331, 1374.26272],
|
|
[ 0.87, 0.895, 0.47873040, 50.357240, 0.089787883, -567.76534, 1503.51828, -1327.11835, 390.45562],
|
|
[ 0.895, 0.93, 0.35255000, 25.498893, -5.4176608, -165.99710, 454.58851, -414.93772, 126.23908],
|
|
[ 0.93, 0.97, 0.32002037, 6.8474560, -8.0901766, -102.84970, 293.25377, -278.69176, 88.27639],
|
|
[ 0.97, 0.99, 0.30845061, -32.343947, -10.455428, -475.69091, 1398.86173, -1371.17931, 448.00538],
|
|
[ 0.99, 0.999, 0.30110729, -333.50440, -15.711260, -11334.69334, 33822.93507, -33642.61426, 11154.37157]]
|
|
|
|
# 0.635 - 0.997
|
|
Coefficient9 = [[ 0.635, 0.645, 58249804, 119.64583, 89.483167, -28072.11478, 53918.28040, -34519.94034, 7366.77050],
|
|
[ 0.645, 0.66, 5100.2910, 70.920100, 46.811007, -9051.10879, 17726.62108, -15572.21470, 2518.10214],
|
|
[ 0.66, 0.685, 19.518159, 38.509198, 19.951285, -2338.00251, 4719.64936, -3175.58038, 712.17229],
|
|
[ 0.685, 0.72, 1.3823760, 20.055606, 5.9746092, -522.67397, 1102.26575, -774.75596, 181.49582],
|
|
[ 0.72, 0.77, 0.49074291, 10.484929, -0.31885445, -111.99730, 250.54286, -186.78287, 46.40581],
|
|
[ 0.77, 0.84, 0.34607161, 5.4063548, -2.9788620, -24.21100, 58.53312, -47.15244, 12.65668],
|
|
[ 0.84, 0.93, 0.31342828, 2.1140055, -4.1151543, -9.16969, 24.37544, -21.58758, 6.36954],
|
|
[ 0.93, 0.98, 0.30397966, -3.6334200, -4.9395699, -29.80626, 85.43510, -81.61766, 25.98669],
|
|
[ 0.98, 0.997, 0.30007321, -34.095653, -6.1531194, -434.11439, 1287.65238, -1273.10762, 419.56923]]
|
|
|
|
# 0.17 - 0.983
|
|
Coefficient10 = [[ 0.17, 0.30, 0.31081678, -0.086278355, -2.7023254, -13.95071, 10.30398, -2.49619, 0.19805],
|
|
[ 0.30, 0.45, 0.30636442, -0.094871613, -2.7246428, 0.38766, -0.44128, 0.16613, -0.02068],
|
|
[ 0.45, 0.68, 0.31156188, -0.063734909, -2.6789763, 0.41537, -0.71168, 0.40321, -0.07551],
|
|
[ 0.68, 0.90, 0.31195754, -0.092720991, -2.6863859, -0.44813, 1.06637, -0.84260, 0.22106],
|
|
[ 0.90, 0.983, 0.30712046, -1.5619311, -2.9809169, -5.74538, 16.23781, -15.29128, 4.79808]]
|
|
|
|
# 0.06 - 0.953
|
|
Coefficient11 = [[ 0.06, 0.08, 1.6030433, 0.043380701, -0.15208454, -1341.32950, 282.27311, -19.71010, 0.45664],
|
|
[ 0.08, 0.13, 0.85491410, -0.056761738, -0.65372414, -160.97221, 50.59168, -5.23625, 0.17842],
|
|
[ 0.13, 0.20, 0.57401485, -0.15265815, -1.0443069, -42.79943, 21.22139, -3.47357, 0.18765],
|
|
[ 0.20, 0.30, 0.41622542, -0.27871823, -1.4464670, -14.08741, 10.58863, -2.63159, 0.21621],
|
|
[ 0.30, 0.44, 0.33587066, -0.42674255, -1.8022482, -3.22540, 3.58964, -1.32194, 0.16106],
|
|
[ 0.44, 0.75, 0.31504207, -0.51197565, -1.9499037, 0.07609, -0.13681, 0.08090, -0.01572],
|
|
[ 0.75, 0.953, 0.31658938, -0.49808609, -1.9290368, -0.24018, 0.61450, -0.52258, 0.14771]]
|
|
|
|
# 0.10 - 0.902
|
|
Coefficient12 = [[ 0.10, 0.12, 103.56308, 1.2019976, 4.4869062, -6625.95443, 2292.46751, -264.28560, 10.15218],
|
|
[ 0.12, 0.14, 2.1585254, 0.34129627, 0.83653215, -742.15388, 289.52312, -37.59979, 1.62555],
|
|
[ 0.14, 0.19, 0.65654048, 0.031877859, -0.37693028, -66.92134, 33.18430, -5.45788, 0.29773],
|
|
[ 0.19, 0.27, 0.45076908, -0.10378254, -0.82836350, -10.66154, 7.36263, -1.68453, 0.12768],
|
|
[ 0.27, 0.40, 0.37085833, -0.21965481, -1.1287016, -2.87493, 2.88810, -0.95966, 0.10546],
|
|
[ 0.40, 0.55, 0.33605594, -0.33317467, -1.3397591, -0.70394, 1.00402, -0.47492, 0.07450],
|
|
[ 0.55, 0.902, 0.32790337, -0.39597877, -1.4180887, 0.00490, -0.01070, 0.00769, -0.00182]]
|
|
|
|
if ( _K == 0.8):
|
|
Coeff = Coefficient8
|
|
if ( _K == 0.9):
|
|
Coeff = Coefficient9
|
|
if ( _K == 1.0):
|
|
Coeff = Coefficient10
|
|
if ( _K == 1.1):
|
|
Coeff = Coefficient11
|
|
if ( _K == 1.2):
|
|
Coeff = Coefficient12
|
|
|
|
# Juza 1996/1997 Equation: 8
|
|
for Entry_k in Coeff:
|
|
if (Entry_k[0] <= _S_diamRatio) and (_S_diamRatio <= Entry_k[1]):
|
|
HInv = Entry_k[2] * _S_diamRatio**(Entry_k[3] * math.log(_S_diamRatio) + Entry_k[4]) + Entry_k[5] * (_S_diamRatio**3.0) + Entry_k[6] * (_S_diamRatio **2.0) + Entry_k[7] * _S_diamRatio + Entry_k[8]
|
|
return HInv
|
|
break
|
|
|
|
return 0
|
|
|
|
|
|
|
|
# ========== USAGE PRINTER ==========
|
|
|
|
def PrintCorrectUsage():
|
|
print("\nDESCRIPTION:\n")
|
|
print("This program will take an image or a set of files (coordinate points) representing")
|
|
print("a droplet and will process the data to calculate the surface tension of the droplet.")
|
|
print("When an image file is set, the program will extract the coordinate of the points")
|
|
print("representing the edge of the droplet. When the files are set, it is assume that they")
|
|
print("will already contain the points representing the edge of the droplet. RK4 will be")
|
|
print("use to solve the Bashforth-Adams set of equations to get the parameters for surface")
|
|
print("tension calculation. parameters: A (apex radius curvature) and B (shape factor)")
|
|
print("\nUSAGE:\n")
|
|
print("Syntax 1: ./img.py -option1 value1 -option2 value2")
|
|
print("Syntax 2: python img.py -option1 value1 -option2 value2\n")
|
|
print("Options:")
|
|
print("-help : Prints this help message. ")
|
|
print("-out : Reference name for all the output files (coordinate files, RK4 dump, etc.).")
|
|
print(" : This should always be set. ")
|
|
print("-img : Image file to be processed. ")
|
|
print(" : When this is set, the file name inputs (-filep, -filen) should not be set.")
|
|
print("-filep : File containing the positive(x) point coordinates to be processed.")
|
|
print(" : If this is set, -filen should also be set and -img should not be set.")
|
|
print("-filen : File containing the negative(x) point coordinates to be processed.")
|
|
print(" : If this is set, -filep should also be set and -img should not be set.")
|
|
print("-proc : The number of processor to use.")
|
|
print(" : If not specified the default value is 2.")
|
|
print("-syrg : This is the syringe width in millimeters.")
|
|
print(" : If not specified the default value is 0.805 mm.")
|
|
print("-a : RK4 Parameter: minimum value for alpha (radius of curvature at apex).")
|
|
print(" : If not specified the default value is 0.0001")
|
|
print("-A : RK4 Parameter: maximum value for alpha (radius of curvature at apex).")
|
|
print(" : If not specified the default value is 2.0001")
|
|
print("-b : RK4 Parameter: minimum value for beta (shape factor).")
|
|
print(" : If not specified the default value is 0.0001")
|
|
print("-B : RK4 Parameter: minimum value for beta (shape factor).")
|
|
print(" : If not specified the default value is 0.5001.")
|
|
print("-x : RK4 Parameter: minimum value for x (initial guessed x coordinate).")
|
|
print(" : If not specified the default value is 0.0001.")
|
|
print("-X : RK4 Parameter: maximum value for x (initial guessed x coordinate).")
|
|
print(" : If not specified the default value is 0.5001.")
|
|
print("-t : RK4 Parameter: minimum value for t (initial guessed angle).")
|
|
print(" : If not specified the default value is 0.0001.")
|
|
print("-T : RK4 Parameter: maximum value for t (initial guessed angle).")
|
|
print(" : If not specified the default value is 0.5001.")
|
|
print("-ai : RK4 Parameter: alpha interval - the number of values to consider between min alpha and max alpha.")
|
|
print(" : If not specified the default value is 20.")
|
|
print("-bi : RK4 Parameter: beta interval - the number of values to consider between min beta and max beta.")
|
|
print(" : If not specified the default value is 20.")
|
|
print("-xi : RK4 Parameter: x interval - the number of values to consider between min x and max x.")
|
|
print(" : If not specified the default value is 5.")
|
|
print("-ti : RK4 Parameter: angle interval - the number of values to consider between min angle and max angle.")
|
|
print(" : If not specified the default value is 5.")
|
|
print("-solver : 1 for divide and conquer scanning.")
|
|
print(" : If not specified, optimization will be used to solve the parameters.")
|
|
print("\nSAMPLE USAGE:\n")
|
|
print("1. Calculate surface tension of a drop using default values.")
|
|
print("./img.py -img water-drop.jpg -out water-drop")
|
|
print("python img.py -img water-drop.jpg -out water-drop\n")
|
|
print("2. Uses 4 processors for calculation.")
|
|
print("./img.py -img water-drop.jpg -out water-drop -proc 4\n")
|
|
print("3. Setting interval sizes for A and B parameters (search for 30 possible values for A and B).")
|
|
print("./img.py -img water-drop.jpg -out water-drop -ai 30 -bi 30\n")
|
|
print("4. Setting min (0.00001) and max (1.5) values for parameter.")
|
|
print("./img.py -img water-drop.jpg -out water-drop -a 0.00001 -A 1.5\n")
|
|
|
|
def ComputeMain(ImageFilePath):
|
|
OutputFileName = ""
|
|
PosPointsFileName = ""
|
|
NegPointsFileName = ""
|
|
ProcessCount = 2
|
|
SyringeWidth = 0.70 * 1.15
|
|
MinA_radCurv = 0.0001
|
|
MaxA_radCurv = 2.0001
|
|
MinB_shapeFact = 0.0001
|
|
MaxB_shapeFact = 1.0001
|
|
MinX = 0.0001
|
|
MaxX = 0.5001
|
|
MinT_angle = 0.0001
|
|
MaxT_angle = 0.5001
|
|
AStep = 20
|
|
BStep = 20
|
|
XStep = 5
|
|
TStep = 5
|
|
|
|
ImageFlag = 0
|
|
PosFileFlag = 0
|
|
NegFileFlag = 0
|
|
OutFlag = 0
|
|
Counter = 1
|
|
SolveMethod = 0
|
|
|
|
ARange_radCurv = []
|
|
a = MinA_radCurv
|
|
ahop = float(MaxA_radCurv - MinA_radCurv)/float(AStep-1)
|
|
while (a <= MaxA_radCurv):
|
|
ARange_radCurv.append(a)
|
|
a += ahop
|
|
if ahop == 0.0:
|
|
break
|
|
if(np.abs(MaxA_radCurv - ARange_radCurv[len(ARange_radCurv)-1] - ahop) < 0.000001 and ahop != 0.0):
|
|
ARange_radCurv.append(MaxA_radCurv)
|
|
|
|
BRange_shapeFact = []
|
|
b = MinB_shapeFact
|
|
bhop = float(MaxB_shapeFact - MinB_shapeFact)/float(BStep-1)
|
|
while (b <= MaxB_shapeFact):
|
|
BRange_shapeFact.append(b)
|
|
b += bhop
|
|
if bhop == 0:
|
|
break
|
|
if(np.abs(MaxB_shapeFact - BRange_shapeFact[len(BRange_shapeFact)-1] - bhop) < 0.000001 and bhop != 0.0):
|
|
BRange_shapeFact.append(MaxB_shapeFact)
|
|
|
|
XRange = []
|
|
x = MinX
|
|
xhop = float(MaxX - MinX)/float(XStep-1)
|
|
while (x <= MaxX):
|
|
XRange.append(x)
|
|
x += xhop
|
|
if xhop == 0:
|
|
break
|
|
if(np.abs(MaxX - XRange[len(XRange)-1] - xhop) < 0.000001 and xhop != 0.0):
|
|
XRange.append(MaxX)
|
|
|
|
TRange_angle = []
|
|
t = MinT_angle
|
|
thop = float(MaxT_angle - MinT_angle)/float(TStep-1)
|
|
while (t <= MaxT_angle):
|
|
TRange_angle.append(t)
|
|
t += thop
|
|
if thop == 0.0:
|
|
break
|
|
if(np.abs(MaxT_angle - TRange_angle[len(TRange_angle)-1] - thop) < 0.000001 and thop != 0.0):
|
|
TRange_angle.append(MaxT_angle)
|
|
|
|
Pos_edgePts, Neg_edgePts, De_scaledMaxWidth, Dms_diamters = ExtractEdgeCoordinate(ImageFilePath, OutputFileName, SyringeWidth)
|
|
HInv1 = HJuza(Dms_diamters[7]/De_scaledMaxWidth, 0.8)
|
|
HInv2 = HJuza(Dms_diamters[8]/De_scaledMaxWidth, 0.9)
|
|
HInv3 = HJuza(Dms_diamters[9]/De_scaledMaxWidth, 1.0)
|
|
HInv4 = HJuza(Dms_diamters[10]/De_scaledMaxWidth, 1.1)
|
|
HInv5 = HJuza(Dms_diamters[11]/De_scaledMaxWidth, 1.2)
|
|
HInvM = (0.20)*(HInv1 + HInv2 + HInv3 + HInv4 + HInv5)
|
|
HInvSD = math.sqrt((0.20)*( (HInv1-HInvM)**2.0 + (HInv2-HInvM)**2.0 + (HInv3-HInvM)**2.0 + (HInv4-HInvM)**2.0 + (HInv5-HInvM)**2.0))
|
|
|
|
def BruteForce():
|
|
PosList = []
|
|
NegList = []
|
|
for ProcessNo in range(ProcessCount):
|
|
PosList.append(copy.deepcopy(Pos_edgePts))
|
|
NegList.append(copy.deepcopy(Neg_edgePts))
|
|
#print(ProcessNo)
|
|
|
|
ARanges = []
|
|
for ProcessNo in range(ProcessCount):
|
|
ARanges.append([])
|
|
|
|
i = 0
|
|
for A in ARange_radCurv:
|
|
ARanges[i%ProcessCount].append(A)
|
|
i += 1
|
|
|
|
BestParams = multiprocessing.Queue()
|
|
TicksQueue = multiprocessing.Queue()
|
|
Processes = []
|
|
|
|
for ProcessNo in range(ProcessCount):
|
|
Processes.append(multiprocessing.Process(name='Search ' + str(ProcessNo), target=FindOptimalParameters, args=(PosList[ProcessNo], NegList[ProcessNo], ARanges[ProcessNo], BRange_shapeFact, XRange, TRange_angle, BestParams, TicksQueue)))
|
|
|
|
for aProcess in Processes:
|
|
aProcess.start()
|
|
|
|
sys.stdout.write("\rParameter Search Progress : 0%")
|
|
sys.stdout.flush()
|
|
Count = 0
|
|
while(Count < len(ARange_radCurv)):
|
|
t = TicksQueue.get()
|
|
Count += 1
|
|
sys.stdout.write("\rParameter Search Progress : " + str(int(100*float(Count)/float(len(ARange_radCurv)))) + "%")
|
|
sys.stdout.flush()
|
|
|
|
for aProcess in Processes:
|
|
aProcess.join()
|
|
|
|
BestOutput = (999999999, 1, 1, 1 ,1)
|
|
while not BestParams.empty():
|
|
Output = BestParams.get()
|
|
if Output[0] < BestOutput[0]:
|
|
BestOutput = Output
|
|
return BestOutput
|
|
|
|
def Optimize():
|
|
return OptimizeParameters(Pos_edgePts, Neg_edgePts, 0.95, 0.3, 0.1, 0.2)
|
|
|
|
starttime = time.time()
|
|
|
|
BestOutput = Optimize() if 1 != SolveMethod else BruteForce()
|
|
|
|
# print("Negative Cooridnate File : " + OutputFileName + "_Negative.dat")
|
|
# print("Positive Cooridnate File : " + OutputFileName + "_Positive.dat")
|
|
# print("Syringe Width (millimeter): " + str(SyringeWidth))
|
|
# print("Radius of Curvature Range : " + str(MinA_radCurv) + "-" + str(MaxA_radCurv) + " (interval: "+str(len(ARange_radCurv))+")")
|
|
# print("Shape Factor Range : " + str(MinB_shapeFact) + "-" + str(MaxB_shapeFact) + " (interval: "+str(len(BRange_shapeFact))+")")
|
|
# print("Initial X Coordinate Range: " + str(MinX) + "-" + str(MaxX) + " (interval: "+str(len(XRange))+")")
|
|
# print("Initial Angle Range : " + str(MinT_angle) + "-" + str(MaxT_angle) + " (interval: "+str(len(TRange_angle))+")")
|
|
# print("Parameter Combinations : " + str(len(ARange_radCurv) * len(BRange_shapeFact) * len(XRange) * len(TRange_angle)))
|
|
# print("Number of Process to Use : " + str(ProcessCount))
|
|
# print("Juza 5-Plane (1/H) Std.Dev: " + str(HInvSD))
|
|
|
|
# print("best output: " + str(BestOutput))
|
|
|
|
# print("\nLowest RMSD : " + str(BestOutput[0]))
|
|
# print("Computed Surface Tension : " + str(9.8*BestOutput[1]/BestOutput[2]))
|
|
# print("Best Radius of Curvature : " + str(BestOutput[1]))
|
|
# print("Best Shape Factor : " + str(BestOutput[2]))
|
|
# print("Best Initial X Coordinate : " + str(BestOutput[3]))
|
|
# print("Best Initial Angle : " + str(BestOutput[4]))
|
|
# print("Lapsed Time in Seconds : " + str(time.time() -starttime))
|
|
|
|
return {
|
|
"negCoordFile": OutputFileName + "_Negative.dat",
|
|
"posCoordFile": OutputFileName + "_Positive.dat",
|
|
"syringeWidth": SyringeWidth,
|
|
"radiusCurvatureRange": str(MinA_radCurv) + "-" + str(MaxA_radCurv) + " (interval: "+str(len(ARange_radCurv))+")",
|
|
"shapeFactorRange": str(MinB_shapeFact) + "-" + str(MaxB_shapeFact) + " (interval: "+str(len(BRange_shapeFact))+")",
|
|
"initialXCoordinateRange": str(MinX) + "-" + str(MaxX) + " (interval: "+str(len(XRange))+")",
|
|
"initialAngleRange": str(MinT_angle) + "-" + str(MaxT_angle) + " (interval: "+str(len(TRange_angle))+")",
|
|
"paramCombinations": len(ARange_radCurv) * len(BRange_shapeFact) * len(XRange) * len(TRange_angle),
|
|
"numProcess": ProcessCount,
|
|
"juza": HInvSD,
|
|
"lowestRMSD": BestOutput[0],
|
|
"computedSurfaceTension": 9.8*BestOutput[1]/BestOutput[2],
|
|
"bestRadiusCurvature": BestOutput[1],
|
|
"bestShapeFactor": BestOutput[2],
|
|
"bestInitialXCoordinate": BestOutput[3],
|
|
"bestInitialAngle": BestOutput[4],
|
|
"lapsedTime": time.time()-starttime
|
|
}
|
|
|
|
|
|
# ========== MAIN ==========
|
|
|
|
if __name__ == '__main__':
|
|
ImageFilePath = '/Users/fsantelices/Desktop/Pendant Drop App/Surface Tension Data/data/ambient-6-26-15-mf-3-30s-MORPHED.jpg'
|
|
print(ComputeMain(ImageFilePath))
|
|
|
|
# ImageFilePath = ""
|
|
# OutputFileName = ""
|
|
# NegPointsFileName = ""
|
|
# ProcessCount = 2
|
|
# SyringeWidth = 0.70 * 1.15
|
|
# MinA_radCurv = 0.0001
|
|
# MaxA_radCurv = 2.0001
|
|
# MinB_shapeFact = 0.0001
|
|
# MaxB_shapeFact = 1.0001
|
|
# MinX = 0.0001
|
|
# MaxX = 0.5001
|
|
# MinT_angle = 0.0001
|
|
# MaxT_angle = 0.5001
|
|
# AStep = 20
|
|
# BStep = 20
|
|
# XStep = 5
|
|
# TStep = 5
|
|
|
|
# ImageFlag = 0
|
|
# PosFileFlag = 0
|
|
# NegFileFlag = 0
|
|
# OutFlag = 0
|
|
# Counter = 1
|
|
# SolveMethod = 0
|
|
|
|
# if (len(sys.argv) <= 1):
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
|
|
# while (Counter <= len(sys.argv)-1):
|
|
|
|
# CurrentOption = sys.argv[Counter]
|
|
# try:
|
|
# CurrentValue = sys.argv[Counter+1]
|
|
# except:
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
|
|
# if CurrentOption == "-help":
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
# elif CurrentOption == "-img":
|
|
# ImageFilePath = CurrentValue
|
|
# ImageFlag = 1
|
|
# elif CurrentOption == "-filep":
|
|
# PosPointsFileName = CurrentValue
|
|
# PosFileFlag = 1
|
|
# elif CurrentOption == "-filen":
|
|
# NegPointsFileName = CurrentValue
|
|
# NegFileFlag = 1
|
|
# elif CurrentOption == "-out":
|
|
# OutputFileName = CurrentValue
|
|
# OutFlag = 1
|
|
# elif CurrentOption == "-proc":
|
|
# ProcessCount = int(CurrentValue)
|
|
# elif CurrentOption == "-syrg":
|
|
# SyringeWidth = float(CurrentValue)
|
|
# elif CurrentOption == "-a":
|
|
# MinA_radCurv = float(CurrentValue)
|
|
# elif CurrentOption == "-A":
|
|
# MaxA_radCurv = float(CurrentValue)
|
|
# elif CurrentOption == "-b":
|
|
# MinB_shapeFact = float(CurrentValue)
|
|
# elif CurrentOption == "-B":
|
|
# MaxB_shapeFact = float(CurrentValue)
|
|
# elif CurrentOption == "-x":
|
|
# MinX = float(CurrentValue)
|
|
# elif CurrentOption == "-X":
|
|
# MaxX = float(CurrentValue)
|
|
# elif CurrentOption == "-t":
|
|
# MinT_angle = float(CurrentValue)
|
|
# elif CurrentOption == "-T":
|
|
# MaxT_angle = float(CurrentValue)
|
|
# elif CurrentOption == "-ai":
|
|
# AStep = int(CurrentValue)
|
|
# elif CurrentOption == "-bi":
|
|
# BStep = int(CurrentValue)
|
|
# elif CurrentOption == "-xi":
|
|
# XStep = int(CurrentValue)
|
|
# elif CurrentOption == "-ti":
|
|
# TStep = int(CurrentValue)
|
|
# elif CurrentOption == "-solver":
|
|
# SolveMethod = int(CurrentValue)
|
|
# else:
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
# Counter = Counter + 2
|
|
|
|
# if (OutFlag == 0 and ImageFlag == 1):
|
|
# print("Output file name not set. (-out)")
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
# if (ImageFlag == 1 and ((PosFileFlag == 1) or (NegFileFlag == 1))):
|
|
# print("Image(-img) can not be set with files option (-filep, -filen).")
|
|
# print("Either set image only or files only.")
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
# if (ImageFlag == 0 and ((PosFileFlag == 1 and NegFileFlag == 0) or (PosFileFlag == 0 and NegFileFlag == 1))):
|
|
# print("File inputs (-filep, -filen) should both be set if image input is not set.")
|
|
# PrintCorrectUsage()
|
|
# sys.exit()
|
|
|
|
# ARange_radCurv = []
|
|
# a = MinA_radCurv
|
|
# ahop = float(MaxA_radCurv - MinA_radCurv)/float(AStep-1)
|
|
# while (a <= MaxA_radCurv):
|
|
# ARange_radCurv.append(a)
|
|
# a += ahop
|
|
# if ahop == 0.0:
|
|
# break
|
|
# if(np.abs(MaxA_radCurv - ARange_radCurv[len(ARange_radCurv)-1] - ahop) < 0.000001 and ahop != 0.0):
|
|
# ARange_radCurv.append(MaxA_radCurv)
|
|
|
|
# BRange_shapeFact = []
|
|
# b = MinB_shapeFact
|
|
# bhop = float(MaxB_shapeFact - MinB_shapeFact)/float(BStep-1)
|
|
# while (b <= MaxB_shapeFact):
|
|
# BRange_shapeFact.append(b)
|
|
# b += bhop
|
|
# if bhop == 0:
|
|
# break
|
|
# if(np.abs(MaxB_shapeFact - BRange_shapeFact[len(BRange_shapeFact)-1] - bhop) < 0.000001 and bhop != 0.0):
|
|
# BRange_shapeFact.append(MaxB_shapeFact)
|
|
|
|
# XRange = []
|
|
# x = MinX
|
|
# xhop = float(MaxX - MinX)/float(XStep-1)
|
|
# while (x <= MaxX):
|
|
# XRange.append(x)
|
|
# x += xhop
|
|
# if xhop == 0:
|
|
# break
|
|
# if(np.abs(MaxX - XRange[len(XRange)-1] - xhop) < 0.000001 and xhop != 0.0):
|
|
# XRange.append(MaxX)
|
|
|
|
# TRange_angle = []
|
|
# t = MinT_angle
|
|
# thop = float(MaxT_angle - MinT_angle)/float(TStep-1)
|
|
# while (t <= MaxT_angle):
|
|
# TRange_angle.append(t)
|
|
# t += thop
|
|
# if thop == 0.0:
|
|
# break
|
|
# if(np.abs(MaxT_angle - TRange_angle[len(TRange_angle)-1] - thop) < 0.000001 and thop != 0.0):
|
|
# TRange_angle.append(MaxT_angle)
|
|
|
|
|
|
# print("\n")
|
|
# if (ImageFlag == 1):
|
|
# print("Image File : " + ImageFilePath)
|
|
# print("Output Name : " + OutputFileName)
|
|
# print("Negative Cooridnate File : " + OutputFileName + "_Negative.dat")
|
|
# print("Positive Cooridnate File : " + OutputFileName + "_Positive.dat")
|
|
# print("Syringe Width (millimeter): " + str(SyringeWidth))
|
|
# else:
|
|
# print("Negative Cooridnate File : " + NegPointsFileName)
|
|
# print("Positive Cooridnate File : " + PosPointsFileName)
|
|
# print("Radius of Curvature Range : " + str(MinA_radCurv) + "-" + str(MaxA_radCurv) + " (interval: "+str(len(ARange_radCurv))+")")
|
|
# print("Shape Factor Range : " + str(MinB_shapeFact) + "-" + str(MaxB_shapeFact) + " (interval: "+str(len(BRange_shapeFact))+")")
|
|
# print("Initial X Coordinate Range: " + str(MinX) + "-" + str(MaxX) + " (interval: "+str(len(XRange))+")")
|
|
# print("Initial Angle Range : " + str(MinT_angle) + "-" + str(MaxT_angle) + " (interval: "+str(len(TRange_angle))+")")
|
|
# print("Parameter Combinations : " + str(len(ARange_radCurv) * len(BRange_shapeFact) * len(XRange) * len(TRange_angle)))
|
|
# print("Number of Process to Use : " + str(ProcessCount))
|
|
|
|
# if (ImageFlag == 1):
|
|
# Pos_edgePts, Neg_edgePts, De_scaledMaxWidth, Dms_diamters = ExtractEdgeCoordinate(ImageFilePath, OutputFileName, SyringeWidth)
|
|
# HInv1 = HJuza(Dms_diamters[7]/De_scaledMaxWidth, 0.8)
|
|
# HInv2 = HJuza(Dms_diamters[8]/De_scaledMaxWidth, 0.9)
|
|
# HInv3 = HJuza(Dms_diamters[9]/De_scaledMaxWidth, 1.0)
|
|
# HInv4 = HJuza(Dms_diamters[10]/De_scaledMaxWidth, 1.1)
|
|
# HInv5 = HJuza(Dms_diamters[11]/De_scaledMaxWidth, 1.2)
|
|
# HInvM = (0.20)*(HInv1 + HInv2 + HInv3 + HInv4 + HInv5)
|
|
# HInvSD = math.sqrt((0.20)*( (HInv1-HInvM)**2.0 + (HInv2-HInvM)**2.0 + (HInv3-HInvM)**2.0 + (HInv4-HInvM)**2.0 + (HInv5-HInvM)**2.0))
|
|
# print("Juza 5-Plane (1/H) Std.Dev: " + str(HInvSD))
|
|
# else:
|
|
# Pos_edgePts = []
|
|
# PositivePointFile = open(PosPointsFileName,'r')
|
|
# for Line in PositivePointFile:
|
|
# Coordinate = Line.split()
|
|
# #print(Coordinate)
|
|
# Pos_edgePts.append((float(Coordinate[0]), float(Coordinate[1])))
|
|
# PositivePointFile.close()
|
|
# Neg_edgePts = []
|
|
# NegativePointFile = open(NegPointsFileName,'r')
|
|
# for Line in NegativePointFile:
|
|
# Coordinate = Line.split()
|
|
# #print(Coordinate)
|
|
# Neg_edgePts.append((float(Coordinate[0]), float(Coordinate[1])))
|
|
# NegativePointFile.close()
|
|
|
|
# def BruteForce():
|
|
# PosList = []
|
|
# NegList = []
|
|
# for ProcessNo in range(ProcessCount):
|
|
# PosList.append(copy.deepcopy(Pos_edgePts))
|
|
# NegList.append(copy.deepcopy(Neg_edgePts))
|
|
# #print(ProcessNo)
|
|
|
|
# ARanges = []
|
|
# for ProcessNo in range(ProcessCount):
|
|
# ARanges.append([])
|
|
|
|
# i = 0
|
|
# for A in ARange_radCurv:
|
|
# ARanges[i%ProcessCount].append(A)
|
|
# i += 1
|
|
|
|
# BestParams = multiprocessing.Queue()
|
|
# TicksQueue = multiprocessing.Queue()
|
|
# Processes = []
|
|
|
|
# for ProcessNo in range(ProcessCount):
|
|
# Processes.append(multiprocessing.Process(name='Search ' + str(ProcessNo), target=FindOptimalParameters, args=(PosList[ProcessNo], NegList[ProcessNo], ARanges[ProcessNo], BRange_shapeFact, XRange, TRange_angle, BestParams, TicksQueue)))
|
|
|
|
# for aProcess in Processes:
|
|
# aProcess.start()
|
|
|
|
# sys.stdout.write("\rParameter Search Progress : 0%")
|
|
# sys.stdout.flush()
|
|
# Count = 0
|
|
# while(Count < len(ARange_radCurv)):
|
|
# t = TicksQueue.get()
|
|
# Count += 1
|
|
# sys.stdout.write("\rParameter Search Progress : " + str(int(100*float(Count)/float(len(ARange_radCurv)))) + "%")
|
|
# sys.stdout.flush()
|
|
|
|
# for aProcess in Processes:
|
|
# aProcess.join()
|
|
|
|
# BestOutput = (999999999, 1, 1, 1 ,1)
|
|
# while not BestParams.empty():
|
|
# Output = BestParams.get()
|
|
# if Output[0] < BestOutput[0]:
|
|
# BestOutput = Output
|
|
# return BestOutput
|
|
|
|
# def Optimize():
|
|
# return OptimizeParameters(Pos_edgePts, Neg_edgePts, 0.95, 0.3, 0.1, 0.2)
|
|
|
|
# starttime = time.time()
|
|
|
|
# BestOutput = Optimize() if 1 != SolveMethod else BruteForce()
|
|
|
|
# print("best output: " + str(BestOutput))
|
|
|
|
# print("\nLowest RMSD : " + str(BestOutput[0]))
|
|
# print("Computed Surface Tension : " + str(9.8*BestOutput[1]/BestOutput[2]))
|
|
# print("Best Radius of Curvature : " + str(BestOutput[1]))
|
|
# print("Best Shape Factor : " + str(BestOutput[2]))
|
|
# print("Best Initial X Coordinate : " + str(BestOutput[3]))
|
|
# print("Best Initial Angle : " + str(BestOutput[4]))
|
|
# print("Lapsed Time in Seconds : " + str(time.time() -starttime))
|