#!/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") # ========== MAIN ========== if __name__ == '__main__': 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 if (len(sys.argv) <= 1): PrintCorrectUsage() sys.exit() ImageFlag = 0 PosFileFlag = 0 NegFileFlag = 0 OutFlag = 0 Counter = 1 SolveMethod = 0 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))