OTSImageProcess.cpp 61 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377
  1. #pragma once
  2. #include "stdafx.h"
  3. #include <opencv2/core/core.hpp>
  4. #include <opencv2/highgui/highgui.hpp>
  5. #include <opencv2/opencv.hpp>
  6. #include "OTSImageProcess.h"
  7. #include "OTSImageProcessParam.h"
  8. #include <OTSFieldData.h>
  9. #include "OTSMorphology.h"
  10. #include "../OTSLog/COTSUtilityDllFunExport.h"
  11. #include "FieldMgr.h"
  12. #include "BaseFunction.h"
  13. #include "GreyPeak.h"
  14. namespace OTSIMGPROC
  15. {
  16. using namespace cv;
  17. using namespace std;
  18. const int nBlackColor = 255;
  19. namespace
  20. {
  21. cv::Mat GetMatDataFromParticles(COTSParticleList parts, int width, int height)
  22. {
  23. int rows, cols;
  24. cols = width;
  25. rows = height;
  26. BYTE* pPixel = new BYTE[width * height];
  27. Mat cvcopyImg = Mat(rows, cols, CV_8UC1, pPixel);
  28. cvcopyImg = Scalar::all(0);
  29. for (auto p : parts)
  30. {
  31. auto area = p->GetPixelArea();
  32. if (area > 100)
  33. {
  34. auto fea = p->GetFeature();
  35. auto segs = fea->GetSegmentsList();
  36. for (auto seg : segs)
  37. {
  38. auto row = seg->GetHeight();
  39. for (int i = 0; i < seg->GetLength(); i++)
  40. {
  41. auto col = seg->GetStart() + i;
  42. cvcopyImg.at<char>(row, col) = 255;
  43. }
  44. }
  45. }
  46. }
  47. return cvcopyImg;
  48. }
  49. }
  50. COTSImageProcess::COTSImageProcess(COTSImageProcessParamPtr a_pImageProcessParam)
  51. {
  52. m_imageProcessParam = a_pImageProcessParam;
  53. }
  54. COTSImageProcess::~COTSImageProcess()
  55. {
  56. }
  57. BOOL COTSImageProcess::RemoveBGByCVconnectivities(CBSEImgPtr inBSEImg, double a_pixelSize, COTSFieldDataPtr m_pFieldData)
  58. {
  59. ASSERT(m_pFieldData);
  60. ASSERT(inBSEImg);
  61. ASSERT(m_imageProcessParam);
  62. int nWidthImg = inBSEImg->GetWidth();
  63. int nHeightImg = inBSEImg->GetHeight();
  64. m_pFieldData->Width = nWidthImg;
  65. m_pFieldData->Height = nHeightImg;
  66. long nImgSize = nWidthImg * nHeightImg;
  67. BYTE* pSrcImg = inBSEImg->GetImageDataPointer();
  68. BYTE* pTempImg = new BYTE[nImgSize];
  69. CRect r = CRect(0, 0, nWidthImg, nHeightImg);
  70. CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
  71. BinaryProcess(inBSEImg, m_imageProcessParam, imgNoBGBinary/*, nNumParticle*/);
  72. BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();
  73. long nEigenGrayStart = m_imageProcessParam->GetParticleGray().GetStart();
  74. long nEigenGrayEnd = m_imageProcessParam->GetParticleGray().GetEnd();
  75. // get the area image
  76. Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
  77. Mat labels = Mat::zeros(cvcopyImg.size(), CV_32S);
  78. Mat stats, centroids;
  79. int number = connectedComponentsWithStats(cvcopyImg, labels, stats, centroids, 8, CV_32S);
  80. double rMin = m_imageProcessParam->GetIncArea().GetStart()/2.0;
  81. double rMax = m_imageProcessParam->GetIncArea().GetEnd()/2.0;
  82. double partAreaMin = rMin * rMin * 3.14159;
  83. double partAreaMax = rMax * rMax * 3.14159;
  84. COTSParticleList listParticleOut;
  85. for (size_t i = 1; i < number; i++)
  86. {
  87. int center_x = centroids.at<double>(i, 0);
  88. int center_y = centroids.at<double>(i, 1);
  89. //矩形边框
  90. int x = stats.at<int>(i, CC_STAT_LEFT);
  91. int y = stats.at<int>(i, CC_STAT_TOP);
  92. int w = stats.at<int>(i, CC_STAT_WIDTH);
  93. int h = stats.at<int>(i, CC_STAT_HEIGHT);
  94. int area = stats.at<int>(i, CC_STAT_AREA);
  95. double actualArea = area * a_pixelSize * a_pixelSize;
  96. if (actualArea >= partAreaMin && actualArea < partAreaMax)
  97. {
  98. Rect rectMax = Rect(x, y, w, h);
  99. Mat rectROI = labels(rectMax).clone();
  100. Mat imageROI = Mat::zeros(rectMax.size(), cvcopyImg.type());
  101. //exclude the point which intersect into this bounding box but is not in this group.
  102. int label = i;
  103. for (int row = 0; row < rectROI.rows; row++)
  104. {
  105. for (int col = 0; col < rectROI.cols; col++)
  106. {
  107. int v = rectROI.at<int>(row, col);
  108. if (v == label)
  109. {
  110. imageROI.at<uchar>(row, col) = 255;//set the value to 255,so we won't consider other pixel when we find segment and feature in this ROI.
  111. }
  112. }
  113. }
  114. COTSParticleList roiParts;
  115. if (GetOneParticleFromROI(rectMax.x, rectMax.y, rectMax.width, rectMax.height, imageROI.data, roiParts))
  116. {
  117. if (roiParts.size() > 0)
  118. {
  119. COTSParticlePtr roiPart = roiParts[0];//we will find only one part in the roi.
  120. roiPart->SetXRayPos(CPoint(center_x, center_y));
  121. CRect r = CRect(x, y, x + w, y + h);
  122. roiPart->SetParticleRect(r);
  123. roiPart->SetActualArea(actualArea);
  124. roiPart->SetPixelArea(area);
  125. listParticleOut.push_back(roiPart);
  126. }
  127. }
  128. }
  129. }
  130. int nTagId;
  131. COTSParticleList listParticleSearched;
  132. for (auto pParticle : listParticleOut)
  133. {
  134. COTSFeaturePtr pFeature = pParticle->GetFeature();
  135. COTSSegmentsList listSegment = pFeature->GetSegmentsList();
  136. long nPixelNum = 0;
  137. long eigenPixelNum = 0;
  138. long nPixelAll = 0;
  139. int nStartS = 0;
  140. int nHeightS = 0;
  141. int nLengthS = 0;
  142. for (auto pSegment : listSegment)
  143. {
  144. // get particle average gray
  145. nStartS = pSegment->GetStart();
  146. nHeightS = pSegment->GetHeight();
  147. nLengthS = pSegment->GetLength();
  148. nPixelNum += (long)nLengthS;
  149. for (unsigned int i = 0; i < nLengthS; i++)
  150. {
  151. long nValueTemp = (long)*(pSrcImg + nHeightS * nWidthImg + nStartS + i);
  152. nPixelAll += nValueTemp;
  153. if (nValueTemp >= nEigenGrayStart && nValueTemp <= nEigenGrayEnd)
  154. {
  155. eigenPixelNum += 1;
  156. }
  157. }
  158. }
  159. double r = m_imageProcessParam->GetIncArea().GetStart() / 2;
  160. double minPartArea = 3.14159 * r * r;
  161. if (eigenPixelNum* a_pixelSize*a_pixelSize > minPartArea)
  162. {
  163. BYTE nAveGray = (BYTE)(nPixelAll / nPixelNum);
  164. pParticle->SetAveGray(nAveGray);
  165. auto fieldOTSRect = m_pFieldData->GetOTSRect();
  166. CPoint leftTop = fieldOTSRect.GetTopLeft();
  167. CRect rectInSinglefld = pParticle->GetParticleRect();
  168. CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_pixelSize, leftTop.y - rectInSinglefld.top * a_pixelSize);
  169. CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_pixelSize, leftTop.y - rectInSinglefld.bottom * a_pixelSize);
  170. COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
  171. pParticle->SetOTSRect(recInOTSCord);
  172. listParticleSearched.push_back(pParticle);
  173. }
  174. }
  175. m_pFieldData->SetParticleList(listParticleSearched);
  176. delete[]pTempImg;
  177. return TRUE;
  178. }
  179. BOOL COTSImageProcess::GetParticlesBySpecialGrayRange(CBSEImgPtr a_pBSEImg, CIntRangePtr a_grayRange,CDoubleRangePtr a_diameterRange,double a_pixelSize, COTSFieldDataPtr m_pFieldData)
  180. {
  181. ASSERT(m_pFieldData);
  182. ASSERT(a_pBSEImg);
  183. ASSERT(a_grayRange);
  184. int nWidthImg = a_pBSEImg->GetWidth();
  185. int nHeightImg = a_pBSEImg->GetHeight();
  186. m_pFieldData->Width = nWidthImg;
  187. m_pFieldData->Height = nHeightImg;
  188. long nImgSize = nWidthImg * nHeightImg;
  189. BYTE* pSrcImg = a_pBSEImg->GetImageDataPointer();
  190. BYTE* pTempImg = new BYTE[nImgSize];
  191. CRect r = CRect(0, 0, nWidthImg, nHeightImg);
  192. CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
  193. long nNumParticle = 0;
  194. GetSpecialGrayRangeImage(a_pBSEImg, *a_grayRange, imgNoBGBinary, nNumParticle);
  195. BYTE* pPixel = imgNoBGBinary->GetImageDataPointer();
  196. if (nNumParticle == 0)
  197. {
  198. COTSParticleList listParticleEmpty;
  199. listParticleEmpty.clear();
  200. m_pFieldData->SetParticleList(listParticleEmpty);
  201. }
  202. else
  203. {
  204. // get the area image
  205. Mat cvcopyImg = Mat(nHeightImg, nWidthImg, CV_8UC1, pPixel);
  206. Mat labels = Mat::zeros(cvcopyImg.size(), CV_32S);
  207. Mat stats, centroids;
  208. int number = connectedComponentsWithStats(cvcopyImg, labels, stats, centroids, 8, CV_32S);
  209. double rStart = a_diameterRange->GetStart() / 2.0;
  210. double rEnd = a_diameterRange->GetEnd() / 2.0;
  211. double areaStart = rStart * rStart * 3.14159;
  212. double areaEnd = rEnd * rEnd * 3.14159;
  213. COTSParticleList listParticleOut;
  214. for (size_t i = 1; i < number; i++)
  215. {
  216. int center_x = centroids.at<double>(i, 0);
  217. int center_y = centroids.at<double>(i, 1);
  218. //矩形边框
  219. int x = stats.at<int>(i, CC_STAT_LEFT);
  220. int y = stats.at<int>(i, CC_STAT_TOP);
  221. int w = stats.at<int>(i, CC_STAT_WIDTH);
  222. int h = stats.at<int>(i, CC_STAT_HEIGHT);
  223. int area = stats.at<int>(i, CC_STAT_AREA);
  224. double actualArea = area * a_pixelSize * a_pixelSize;
  225. if (actualArea >= areaStart && actualArea < areaEnd)
  226. {
  227. Rect rectMax = Rect(x, y, w, h);
  228. Mat rectROI = labels(rectMax).clone();
  229. Mat imageROI = Mat::zeros(rectMax.size(), cvcopyImg.type());
  230. //exclude the point which intersect into this bounding box but is not in this group.
  231. int label = i;
  232. for (int row = 0; row < rectROI.rows; row++)
  233. {
  234. for (int col = 0; col < rectROI.cols; col++)
  235. {
  236. int v = rectROI.at<int>(row, col);
  237. if (v == label)
  238. {
  239. imageROI.at<uchar>(row, col) = 255;
  240. }
  241. }
  242. }
  243. COTSParticleList roiParts;
  244. if (!GetOneParticleFromROI(rectMax.x, rectMax.y, rectMax.width, rectMax.height, imageROI.data, roiParts))
  245. {
  246. continue;
  247. }
  248. if (roiParts.size() > 0)
  249. {
  250. COTSParticlePtr roiPart = roiParts[0];
  251. roiPart->SetXRayPos(CPoint(center_x, center_y));
  252. CRect r = CRect(x, y, x + w, y + h);
  253. roiPart->SetParticleRect(r);
  254. roiPart->SetActualArea(actualArea);
  255. roiPart->SetPixelArea(area);
  256. listParticleOut.push_back(roiPart);
  257. }
  258. }
  259. }
  260. // form a image only have particles on
  261. //COTSSegmentsList listImage;
  262. for (auto pParticle : listParticleOut)
  263. {
  264. int area = pParticle->GetActualArea();
  265. double pActualArea = area ;
  266. COTSFeaturePtr pFeature = pParticle->GetFeature();
  267. COTSSegmentsList listSegment = pFeature->GetSegmentsList();
  268. long nPixelNum = 0;
  269. long nPixelAll = 0;
  270. int nStartS = 0;
  271. int nHeightS = 0;
  272. int nLengthS = 0;
  273. for (auto pSegment : listSegment)
  274. {
  275. // get particle average gray
  276. nStartS = pSegment->GetStart();
  277. nHeightS = pSegment->GetHeight();
  278. nLengthS = pSegment->GetLength();
  279. nPixelNum += (long)nLengthS;
  280. for (unsigned int i = 0; i < nLengthS; i++)
  281. {
  282. long nValueTemp = (long)*(pSrcImg + nHeightS * nWidthImg + nStartS + i);
  283. nPixelAll += nValueTemp;
  284. }
  285. }
  286. BYTE nAveGray = (BYTE)(nPixelAll / nPixelNum);
  287. pParticle->SetAveGray(nAveGray);
  288. auto fieldOTSRect = m_pFieldData->GetOTSRect();
  289. CPoint leftTop = fieldOTSRect.GetTopLeft();
  290. CRect rectInSinglefld = pParticle->GetParticleRect();
  291. CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_pixelSize, leftTop.y - rectInSinglefld.top * a_pixelSize);
  292. CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_pixelSize, leftTop.y - rectInSinglefld.bottom * a_pixelSize);
  293. COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
  294. pParticle->SetOTSRect(recInOTSCord);
  295. }
  296. m_pFieldData->SetParticleList(listParticleOut);
  297. }
  298. delete[]pTempImg;
  299. return TRUE;
  300. }
  301. BOOL COTSImageProcess::SplitFieldImageIntoMatricsParticle(CBSEImgPtr fieldImg, double a_PixelSize, COTSFieldDataPtr outFldData)
  302. {
  303. int nWidthImg = fieldImg->GetWidth();
  304. int nHeightImg = fieldImg->GetHeight();
  305. outFldData->Width = nWidthImg;
  306. outFldData->Height = nHeightImg;
  307. CRect r = CRect(0, 0, nWidthImg, nHeightImg);
  308. CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
  309. double xrayStep = m_imageProcessParam->GetXrayStep();
  310. xrayStep = xrayStep / a_PixelSize;
  311. auto rect = fieldImg->GetImageRect();
  312. std::vector<CPoint> matrixPs;
  313. GetMatrixPointsFromRect(rect, xrayStep, matrixPs);
  314. //int colnum = ceil((double)rect.Width() / xrayStep + 0.5);
  315. //if (colnum % 2 == 0) colnum += 1;//let the number to be an odd number.Then we can make the middle point in the center of the particle exactly.
  316. //int rownum = ceil((double)rect.Height() / xrayStep + 0.5);
  317. //if (rownum % 2 == 0) rownum += 1;
  318. //CPoint theFirst = CPoint(rect.left - (colnum * xrayStep - rect.Width()) / 2 + xrayStep / 2, rect.top - (rownum * xrayStep - rect.Height()) / 2 + xrayStep / 2);
  319. //for (int i = 0; i < rownum; i++)
  320. //{
  321. // for (int j = 0; j < colnum; j++)
  322. // {
  323. // double x = (double)theFirst.x + (double)j * xrayStep;
  324. // double y = (double)theFirst.y + (double)i * xrayStep;
  325. // CPoint thePoint = CPoint(x, y);
  326. // matrixPs.push_back(thePoint);
  327. // }
  328. //}
  329. int i = 0;
  330. COTSParticleList matrixParts;
  331. for (auto point : matrixPs)
  332. {
  333. COTSParticlePtr part = COTSParticlePtr(new COTSParticle());
  334. COTSSegmentsList segs;
  335. for (int i = 0; i < xrayStep; i++)
  336. {
  337. COTSSegmentPtr seg = COTSSegmentPtr(new COTSSegment());
  338. seg->SetStart(point.x - xrayStep / 2);
  339. seg->SetLength(xrayStep);
  340. seg->SetHeight(point.y - xrayStep / 2 + i);
  341. //auto originalSegs = segsOnTheSameHeight[seg->GetHeight()];
  342. int currentH = seg->GetHeight();
  343. int segStart = seg->GetStart();
  344. int segEnd = seg->GetEnd();
  345. if (segStart < 0) segStart = 0;
  346. if (segEnd >= rect.Width()) segEnd = rect.Width()-1;
  347. if (currentH >= rect.Height()) continue;
  348. if (currentH < 0) continue;
  349. seg->SetStart(segStart);
  350. seg->SetEnd(segEnd);
  351. segs.push_back(seg);
  352. }
  353. if (segs.size() > 0)
  354. {
  355. part->GetFeature()->SetSegmentsList(segs);
  356. part->SetActualArea(xrayStep * xrayStep);
  357. part->CalXRayPos();
  358. part->SetFieldId(outFldData->GetId());
  359. part->SetAnalysisId(i);
  360. matrixParts.push_back(part);
  361. i++;
  362. }
  363. }
  364. outFldData->SetParticleList(matrixParts);
  365. return true;
  366. //-----------
  367. }
  368. BOOL COTSImageProcess::SplitRawImageIntoParticlesByWaterShed(CBSEImgPtr fieldImg, double a_PixelSize, COTSFieldDataPtr outFldData)
  369. {
  370. int nWidthImg = fieldImg->GetWidth();
  371. int nHeightImg = fieldImg->GetHeight();
  372. outFldData->Width = nWidthImg;
  373. outFldData->Height = nHeightImg;
  374. CRect r = CRect(0, 0, nWidthImg, nHeightImg);
  375. CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
  376. COTSParticleList listMainParticle;
  377. BinaryProcess(fieldImg, m_imageProcessParam, imgNoBGBinary);
  378. CDoubleRange areaRng = m_imageProcessParam->GetIncArea();
  379. GetParticlesFromBinaryImage(imgNoBGBinary, areaRng, a_PixelSize, listMainParticle);
  380. outFldData->SetParticleList(listMainParticle);
  381. CBSEImgPtr noBgImg = CBSEImgPtr(new CBSEImg(fieldImg->GetImageRect()));
  382. RemoveBackGround(fieldImg, noBgImg);
  383. BlurImage(noBgImg);
  384. std::vector<CPoint> seeds;
  385. seeds.clear();
  386. FindSeedsByGrayScale(noBgImg, a_PixelSize, seeds);
  387. std::vector<CPoint> matrixseeds;
  388. matrixseeds.clear();
  389. if (seeds.size() <20)
  390. {
  391. auto step = areaRng.GetStart() * 4;
  392. GetMatrixPointsFromRect(fieldImg->GetImageRect(), step / a_PixelSize, matrixseeds);
  393. seeds.insert(seeds.end(), matrixseeds.begin(), matrixseeds.end());
  394. }
  395. Mat matImg = GetMatDataFromBseImg(noBgImg);
  396. Mat marks(matImg.size(), CV_32S);
  397. marks = Scalar::all(0);
  398. for (int index = 0; index < seeds.size(); index++)
  399. {
  400. auto p = seeds[index];
  401. marks.at<INT32>(p.y, p.x) = index + 10;
  402. }
  403. Mat imageGray3;
  404. cv::cvtColor(matImg, imageGray3, cv::COLOR_GRAY2RGB);//灰度转换 opencv4
  405. cv::watershed(imageGray3, marks);
  406. auto matBinNobg = GetMatDataFromParticles(listMainParticle, nWidthImg, nHeightImg);
  407. //the first "long" is the index of particle,the second "long" is the row number of the segment, the third "long" is the sequence of the pixel.
  408. std::map<long, std::map<long, std::set<long>>> pixelMap;
  409. for (int i = 0; i < marks.rows; i++)
  410. {
  411. for (int j = 0; j < marks.cols; j++)
  412. {
  413. if (matBinNobg.at<char>(i, j) == 0)
  414. {
  415. marks.at<int>(i, j) = 0;
  416. continue;
  417. }
  418. long index = marks.at<int>(i, j);
  419. if (index == -1)//if this pixel is on the border,then it will be absorbed by it's neighbor.
  420. {
  421. if (i + 1 < marks.rows && marks.at<int>(i + 1, j) != -1 && matBinNobg.at<char>(i + 1, j) != 0)
  422. {
  423. index = marks.at<int>(i + 1, j);
  424. }
  425. else if (j + 1 < marks.cols && marks.at<int>(i, j + 1) != -1 && matBinNobg.at<char>(i, j + 1) != 0)
  426. {
  427. index = marks.at<int>(i, j + 1);
  428. }
  429. else if (i > 0 && marks.at<int>(i - 1, j) != -1 && matBinNobg.at<char>(i - 1, j) != 0)
  430. {
  431. index = marks.at<int>(i - 1, j);
  432. }
  433. else if (j > 0 && marks.at<int>(i, j - 1) != -1 && matBinNobg.at<char>(i, j - 1) != 0)
  434. {
  435. index = marks.at<int>(i, j - 1);
  436. }
  437. marks.at<int>(i, j) = index;
  438. }
  439. pixelMap[index][i].insert(j);
  440. }
  441. }
  442. COTSParticleList listParticleOut;
  443. for (auto it = pixelMap.begin(); it != pixelMap.end(); it++)
  444. {
  445. if (it->first == -1) continue;
  446. if (it->first == 0) continue;
  447. auto pixelsegs = it->second;
  448. COTSParticlePtr newsubpart = COTSParticlePtr(new COTSParticle());
  449. COTSFeaturePtr newFea = COTSFeaturePtr(new COTSFeature());
  450. for (auto pixelseg : pixelsegs)
  451. {
  452. auto start = (*pixelseg.second.begin());
  453. for (auto pix = pixelseg.second.begin(); pix != pixelseg.second.end(); pix++)
  454. {
  455. auto curNode = pix;
  456. ++curNode;
  457. if ((curNode) != pixelseg.second.end())
  458. {
  459. if (*curNode - *pix > 1)
  460. {
  461. COTSSegmentPtr seg = COTSSegmentPtr(new COTSSegment());
  462. seg->SetStart(start);
  463. seg->SetEnd(*pix);
  464. seg->SetHeight(pixelseg.first);
  465. newFea->AddSegment(seg);
  466. start = *curNode;
  467. }
  468. }
  469. else
  470. {
  471. COTSSegmentPtr seg = COTSSegmentPtr(new COTSSegment());
  472. seg->SetStart(start);
  473. seg->SetEnd(*pix);
  474. seg->SetHeight(pixelseg.first);
  475. newFea->AddSegment(seg);
  476. }
  477. }
  478. }
  479. newsubpart->SetFeature(newFea);
  480. newsubpart->CalXRayPos();
  481. newsubpart->CalActualArea(a_PixelSize);
  482. CPoint pos = newsubpart->GetXRayPos();
  483. for (auto mp : listMainParticle)
  484. {
  485. auto fieldOTSRect = outFldData->GetOTSRect();
  486. CPoint leftTop = fieldOTSRect.GetTopLeft();
  487. CRect rectInSinglefld = mp->GetParticleRect();
  488. CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_PixelSize, leftTop.y - rectInSinglefld.top * a_PixelSize);
  489. CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_PixelSize, leftTop.y - rectInSinglefld.bottom * a_PixelSize);
  490. COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
  491. mp->SetOTSRect(recInOTSCord);
  492. mp->SetFieldId(outFldData->GetId());
  493. if (mp->IfContain(pos))
  494. {
  495. CRect rectInSinglefld = newsubpart->GetParticleRect();
  496. CPoint OTSLeftTop = CPoint(leftTop.x + rectInSinglefld.left * a_PixelSize, leftTop.y - rectInSinglefld.top * a_PixelSize);
  497. CPoint OTSRightBottom = CPoint(leftTop.x + rectInSinglefld.right * a_PixelSize, leftTop.y - rectInSinglefld.bottom * a_PixelSize);
  498. COTSRect recInOTSCord = COTSRect(OTSLeftTop, OTSRightBottom);
  499. newsubpart->SetFieldId(mp->GetFieldId());
  500. newsubpart->SetOTSRect(recInOTSCord);
  501. mp->AddSubParticle(newsubpart);
  502. break;
  503. }
  504. }
  505. }
  506. return true;
  507. }
  508. CIntRangePtr COTSImageProcess::CalBackground(CBSEImgPtr m_pBSEImg)
  509. {
  510. auto ranges = CalcuGrayLevelRange(m_pBSEImg);
  511. return ranges[0];
  512. }
  513. std::vector<CIntRangePtr> COTSImageProcess::CalcuGrayLevelRange(CBSEImgPtr m_pBSEImg)
  514. {
  515. CIntRangePtr pBackground = CIntRangePtr(new CIntRange());
  516. WORD originChartData[MAXBYTE];
  517. WORD firstSmoothChart[MAXBYTE];
  518. WORD secondSmooth[MAXBYTE];
  519. //1. get chart data
  520. m_pBSEImg->SetChartData();
  521. memcpy(originChartData, m_pBSEImg->GetBSEChart(), sizeof(WORD) * MAXBYTE);
  522. originChartData[0] = 0;
  523. originChartData[254] = 0;
  524. linearSmooth5(originChartData, firstSmoothChart, MAXBYTE);
  525. linearSmooth5(firstSmoothChart, secondSmooth, MAXBYTE);
  526. //2. get down edge
  527. int nLengthEdge = MAXBYTE + 2;
  528. WORD n_aBSEChart[MAXBYTE + 2];
  529. memset(n_aBSEChart, 0, sizeof(WORD) * nLengthEdge);
  530. std::map<long, std::vector<int>> peakMap;// hold all the peaks in this spectrum which are sorted by there area.
  531. std::vector<int> currentUpSeries;
  532. std::vector<int> currentPeakSeries;
  533. // make sure the wave begin with up edge and end with down edge
  534. n_aBSEChart[0] = 0;
  535. n_aBSEChart[nLengthEdge - 1] = 0;
  536. memcpy(&n_aBSEChart[1], &secondSmooth, sizeof(WORD) * MAXBYTE);
  537. int nLengthCom = MAXBYTE + 1;
  538. // up edge
  539. for (int i = 0; i < nLengthCom; i++)
  540. {
  541. if (n_aBSEChart[i] <= n_aBSEChart[i + 1])//this is a upward edge
  542. {
  543. if (currentPeakSeries.size() > 0)
  544. {
  545. int seriesSize = currentPeakSeries.size();
  546. long area = 0;
  547. for (int i = 0; i < seriesSize; i++)
  548. {
  549. area = area + n_aBSEChart[currentPeakSeries[i]];
  550. }
  551. peakMap[area] = currentPeakSeries;
  552. currentPeakSeries.clear();
  553. }
  554. currentUpSeries.push_back(i + 1);// save all the continuous up edge
  555. }
  556. else//this is a downward edge
  557. {
  558. // encounter a downward edge means upward edge series end,
  559. if (currentUpSeries.size() > 0)
  560. {
  561. currentPeakSeries = currentUpSeries;
  562. currentUpSeries.clear();
  563. }
  564. currentPeakSeries.push_back(i + 1);
  565. }
  566. }
  567. if (currentPeakSeries.size() > 0)
  568. {
  569. int seriesSize = currentPeakSeries.size();
  570. long area = 0;
  571. for (int i = 0; i < seriesSize; i++)
  572. {
  573. area = area + n_aBSEChart[currentPeakSeries[i]];
  574. }
  575. peakMap[area] = currentPeakSeries;
  576. currentPeakSeries.clear();
  577. }
  578. std::vector<CIntRangePtr> ranges;
  579. std::map<long, std::vector<int>>::reverse_iterator it;
  580. for (it=peakMap.rbegin();it!=peakMap.rend();it++)
  581. {
  582. CIntRangePtr pRange = CIntRangePtr(new CIntRange());
  583. pRange->SetStart(it->second[0]);
  584. pRange->SetEnd(it->second[it->second.size()-1]);
  585. ranges.push_back(pRange);
  586. }
  587. return ranges;
  588. }
  589. void COTSImageProcess::GetSpecialGrayRangeImage(CBSEImgPtr a_pImgIn, CIntRange a_SpecialGrayRange, CBSEImgPtr a_pBinImgOut, long& foundedPixelNum)
  590. {
  591. // the background pixel will be 0,and the other part will be 255.
  592. ASSERT(a_pImgIn);
  593. int nWidthImg = a_pImgIn->GetWidth();
  594. int nHeightImg = a_pImgIn->GetHeight();
  595. long nImgSize = nWidthImg * nHeightImg;
  596. BYTE* pTempImg = new BYTE[nImgSize];
  597. BYTE* pSrcImg = a_pImgIn->GetImageDataPointer();
  598. BYTE* pPixel = new byte[nImgSize];
  599. long nBGStart;
  600. long nBGEnd;
  601. long nNumParticle = 0;
  602. nBGStart = a_SpecialGrayRange.GetStart();
  603. nBGEnd = a_SpecialGrayRange.GetEnd();
  604. // delete background
  605. for (unsigned int i = 0; i < nImgSize; i++)
  606. {
  607. if (pSrcImg[i] >= nBGStart && pSrcImg[i] <= nBGEnd)
  608. {
  609. pPixel[i] = 255;
  610. nNumParticle++;
  611. }
  612. else
  613. {
  614. pPixel[i] = 0;
  615. }
  616. }
  617. BErode3(pPixel, pTempImg, 5, nHeightImg, nWidthImg);
  618. BDilate3(pTempImg, pPixel, 5, nHeightImg, nWidthImg);
  619. a_pBinImgOut->SetImageData(pPixel, nWidthImg, nHeightImg);
  620. foundedPixelNum = nNumParticle;
  621. delete[] pTempImg;
  622. return;
  623. }
  624. void COTSImageProcess::BinaryProcess(CBSEImgPtr a_pImgIn, COTSImageProcessParamPtr a_pImageProcessParam, CBSEImgPtr a_pBinImgOut/*,long& foundedPixelNum*/)
  625. {
  626. // the background pixel will be 0,and the other part will be 255.
  627. ASSERT(a_pImgIn);
  628. ASSERT(a_pImageProcessParam);
  629. int nWidthImg = a_pImgIn->GetWidth();
  630. int nHeightImg = a_pImgIn->GetHeight();
  631. long nImgSize = nWidthImg * nHeightImg;
  632. BYTE* pSrcImg = a_pImgIn->GetImageDataPointer();
  633. BYTE* pPixel= new BYTE[nImgSize];
  634. Mat srcImgMat = GetMatDataFromBseImg(a_pImgIn);
  635. Mat rstMat;
  636. long nBGStart;
  637. long nBGEnd;
  638. long nPartStart;
  639. long nPartEnd;
  640. long nNumParticle = 0;
  641. if (a_pImageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::MANUAL || a_pImageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::WaterShed)
  642. {
  643. nBGStart = a_pImageProcessParam->GetBGGray().GetStart();
  644. nBGEnd = a_pImageProcessParam->GetBGGray().GetEnd();
  645. nPartStart = a_pImageProcessParam->GetParticleGray().GetStart();
  646. nPartEnd = a_pImageProcessParam->GetParticleGray().GetEnd();
  647. //CVRemoveBG(srcImgMat, rstMat, nBGStart, nBGEnd);
  648. RemoveBG_old(srcImgMat, rstMat, nBGStart, nBGEnd, nNumParticle);
  649. pPixel = rstMat.data;
  650. }
  651. else
  652. {
  653. switch (a_pImageProcessParam->GetAutoBGRemoveType())
  654. {
  655. case OTS_AUTOBGREMOVE_TYPE::DOWNWARD:
  656. //RemoveBG_old(srcImgMat, rstMat, 0, nBGEnd, nNumParticle);
  657. //CVRemoveBG(srcImgMat, rstMat, 0, nBGEnd, nNumParticle);
  658. AutoRemove_background_OTS(srcImgMat, rstMat, 1);
  659. break;
  660. case OTS_AUTOBGREMOVE_TYPE::UPWARD:
  661. //RemoveBG_old(srcImgMat, rstMat, nBGStart, 255, nNumParticle);
  662. //CVRemoveBG(srcImgMat, rstMat, nBGStart, 255, nNumParticle);
  663. AutoRemove_background_OTS(srcImgMat, rstMat, 0);
  664. break;
  665. case OTS_AUTOBGREMOVE_TYPE::MIDDLE:
  666. //RemoveBG_old(srcImgMat, rstMat, nBGStart, nBGEnd, nNumParticle);
  667. //CVRemoveBG(srcImgMat, rstMat, nBGStart, nBGEnd, nNumParticle);
  668. AutoRemove_background_OTS(srcImgMat, rstMat, 2);
  669. break;
  670. default:
  671. break;
  672. }
  673. pPixel = rstMat.data;
  674. }
  675. a_pBinImgOut->SetImageData(pPixel,nWidthImg,nHeightImg);
  676. return ;
  677. }
  678. void COTSImageProcess::RemoveBackGround(CBSEImgPtr a_pImgIn, CBSEImgPtr a_pBinImgOut)
  679. {
  680. // the background pixel will be 0,and the other part will be 255.
  681. ASSERT(a_pImgIn);
  682. ASSERT(m_imageProcessParam);
  683. int nWidthImg = a_pImgIn->GetWidth();
  684. int nHeightImg = a_pImgIn->GetHeight();
  685. long nImgSize = nWidthImg * nHeightImg;
  686. BYTE* pTempImg = new BYTE[nImgSize];
  687. BYTE* pSrcImg = a_pImgIn->GetImageDataPointer();
  688. BYTE* pPixel = new byte[nImgSize];
  689. long nBGStart;
  690. long nBGEnd;
  691. long nPartStart;
  692. long nPartEnd;
  693. long nNumParticle = 0;
  694. if (m_imageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::MANUAL|| m_imageProcessParam->GetBGRemoveType() == OTS_BGREMOVE_TYPE::WaterShed)
  695. {
  696. nBGStart = m_imageProcessParam->GetBGGray().GetStart();
  697. nBGEnd = m_imageProcessParam->GetBGGray().GetEnd();
  698. nPartStart = m_imageProcessParam->GetParticleGray().GetStart();
  699. nPartEnd = m_imageProcessParam->GetParticleGray().GetEnd();
  700. // delete background
  701. for (unsigned int i = 0; i < nImgSize; i++)
  702. {
  703. if (pSrcImg[i] >= nBGStart && pSrcImg[i] <= nBGEnd)
  704. {
  705. pPixel[i] = 0;
  706. }
  707. else
  708. {
  709. pPixel[i] = pSrcImg[i];
  710. nNumParticle++;
  711. }
  712. if (pSrcImg[i]<nPartStart || pSrcImg[i]>nPartEnd)
  713. {
  714. pPixel[i] = 0;
  715. }
  716. }
  717. }
  718. a_pBinImgOut->SetImageData(pPixel, nWidthImg, nHeightImg);
  719. delete[] pTempImg;
  720. return;
  721. }
  722. BOOL COTSImageProcess::GetParticles(long left, long top, long a_nWidth, long a_nHeight, const BYTE* a_pPixel, COTSParticleList& a_listParticles)
  723. {
  724. ASSERT(a_pPixel);
  725. if (!a_pPixel)
  726. {
  727. return FALSE;
  728. }
  729. //a_listParticles.clear();
  730. COTSParticleList findedParts;
  731. COTSSegmentsList listSegment;
  732. listSegment.clear();
  733. //1. get segment line by line
  734. if (!GetSegmentList(left, top, a_nWidth, a_nHeight, a_pPixel, listSegment))
  735. {
  736. return FALSE;
  737. }
  738. if ((int)listSegment.size() == 0)
  739. {
  740. return FALSE;
  741. }
  742. //2. save the temp feature
  743. COTSFeatureList listFeature;
  744. listFeature.clear();
  745. if (!GetFeatureList(listSegment, listFeature))//get every feature for all the particle,the complete feature.
  746. {
  747. return FALSE;
  748. }
  749. if ((int)listFeature.size() == 0)
  750. {
  751. return FALSE;
  752. }
  753. /*COTSParticleList listParticles;
  754. listParticles.clear();*/
  755. if (!ChangeFeaturelist(listFeature, findedParts))
  756. {
  757. return FALSE;
  758. }
  759. for (auto f : findedParts)
  760. {
  761. a_listParticles.push_back(f);
  762. }
  763. return TRUE;
  764. }
  765. BOOL COTSImageProcess::GetOneParticleFromROI(long left, long top, long a_nWidth, long a_nHeight, const BYTE* a_pPixel, COTSParticleList& a_listParticles)
  766. {
  767. ASSERT(a_pPixel);
  768. if (!a_pPixel)
  769. {
  770. return FALSE;
  771. }
  772. //a_listParticles.clear();
  773. COTSParticleList findedParts;
  774. COTSSegmentsList listSegment;
  775. listSegment.clear();
  776. //1. get segment line by line
  777. if (!GetSegmentList(left, top, a_nWidth, a_nHeight, a_pPixel, listSegment))
  778. {
  779. return FALSE;
  780. }
  781. if ((int)listSegment.size() == 0)
  782. {
  783. return FALSE;
  784. }
  785. //2. save the temp feature
  786. COTSFeatureList listFeature;
  787. listFeature.clear();
  788. COTSFeaturePtr fea = COTSFeaturePtr(new COTSFeature());
  789. fea->SetSegmentsList(listSegment);
  790. listFeature.push_back(fea);
  791. if ((int)listFeature.size() == 0)
  792. {
  793. return FALSE;
  794. }
  795. if (!ChangeFeaturelist(listFeature, findedParts))
  796. {
  797. return FALSE;
  798. }
  799. for (auto f : findedParts)
  800. {
  801. a_listParticles.push_back(f);
  802. }
  803. return TRUE;
  804. }
  805. BOOL COTSImageProcess::GetSegmentList(long left, long top, long a_nWidth, long a_nHeight, const BYTE* a_pPixel, COTSSegmentsList& a_listSegments)
  806. {
  807. ASSERT(a_pPixel);
  808. long nImgSize = a_nWidth * a_nHeight;
  809. a_listSegments.clear();
  810. //1. get segment line by line
  811. long nLine, nm, nn;
  812. long nStart = 0, nLength = 0;
  813. for (nLine = 0; nLine < a_nHeight; nLine++)
  814. {
  815. for (nm = 0; nm < a_nWidth; nm += (nLength + 1))
  816. {
  817. nLength = 0;
  818. // get start
  819. if (*(a_pPixel + nLine * a_nWidth + nm) != 0)
  820. {
  821. nStart = nm;
  822. nLength++;
  823. //get length
  824. for (nn = nm + 1; nn < a_nWidth; nn++)
  825. {
  826. // check if segment is over, break
  827. if (nLength != 0)
  828. {
  829. if (*(a_pPixel + nLine * a_nWidth + nn) == 0)
  830. break;
  831. }
  832. if (*(a_pPixel + nLine * a_nWidth + nn) != 0)
  833. {
  834. nLength++;
  835. }
  836. }
  837. // generate segment
  838. COTSSegmentPtr pSegment = COTSSegmentPtr(new COTSSegment(nLine + top, nStart + left, nLength));
  839. a_listSegments.push_back(pSegment);
  840. }
  841. else
  842. {
  843. continue;
  844. }
  845. }
  846. }
  847. if ((int)a_listSegments.size() == 0)
  848. {
  849. //LogErrorTrace(__FILE__, __LINE__, _T("no particle is found."));
  850. return FALSE;
  851. }
  852. return TRUE;
  853. }
  854. BOOL COTSImageProcess::GetFeatureList(COTSSegmentsList& a_listSegments, COTSFeatureList& a_listFeatures)
  855. {
  856. COTSSegmentsList listSegmentNew;
  857. std::map<long, COTSSegmentsList > mapOneLineSegments;
  858. for each (auto s in a_listSegments)
  859. {
  860. mapOneLineSegments[s->GetHeight()].push_back(s);//sorting all the segments base on the line number.
  861. }
  862. std::map<long, COTSSegmentsList >::iterator lineItr = mapOneLineSegments.begin();//find the highest line
  863. while (lineItr != mapOneLineSegments.end())
  864. {
  865. for (auto s = lineItr->second.begin(); s < lineItr->second.end(); )//find one segment of this line.
  866. {
  867. COTSSegmentPtr bottomSeg = *s;
  868. listSegmentNew.clear();
  869. listSegmentNew.push_back(*s);
  870. s = lineItr->second.erase(s);
  871. std::map<long, COTSSegmentsList >::iterator tempItr = lineItr;
  872. tempItr++;
  873. for (; tempItr != mapOneLineSegments.end(); tempItr++)//find all other lines of segments
  874. {
  875. if (tempItr->first - bottomSeg->GetHeight() > 1)
  876. {
  877. break;
  878. }
  879. for (auto nextLineSegment = tempItr->second.begin(); nextLineSegment < tempItr->second.end();)//find next line's all segments
  880. {
  881. if (((*nextLineSegment)->GetStart() - (bottomSeg->GetStart() + bottomSeg->GetLength())) > 1)
  882. {
  883. break;
  884. }
  885. if (bottomSeg->UpDownConection(**nextLineSegment))
  886. {
  887. listSegmentNew.push_back(*nextLineSegment);
  888. bottomSeg = *nextLineSegment;
  889. nextLineSegment = tempItr->second.erase(nextLineSegment);
  890. break;
  891. }
  892. if (tempItr->second.size() > 0)
  893. {
  894. nextLineSegment++;
  895. }
  896. else
  897. {
  898. break;
  899. }
  900. }
  901. }
  902. COTSFeaturePtr pFeature = COTSFeaturePtr(new COTSFeature());
  903. pFeature->SetSegmentsList(listSegmentNew);
  904. //check if this new feature is connected with other found feature.
  905. COTSSegmentPtr topSeg = listSegmentNew[0];//find the toppest segment of this new feature.
  906. COTSSegmentPtr bottomSegment = listSegmentNew[listSegmentNew.size() - 1];//find the lowest segment of this new feature.
  907. bool haveMerged = false;
  908. for each (auto f in a_listFeatures)
  909. {
  910. for (auto seg : f->GetSegmentsList())
  911. {
  912. if (bottomSegment->UpDownConection(*seg) || topSeg->UpDownConection(*seg))
  913. {
  914. COTSSegmentsList segs = f->GetSegmentsList();
  915. for (auto s : listSegmentNew)
  916. {
  917. segs.push_back(s);
  918. }
  919. f->SetSegmentsList(segs);
  920. haveMerged = true;
  921. break;
  922. }
  923. }
  924. if (haveMerged)
  925. {
  926. break;
  927. }
  928. }
  929. if (!haveMerged)
  930. {
  931. a_listFeatures.push_back(pFeature);
  932. }
  933. if (lineItr->second.size() == 0)
  934. {
  935. break;
  936. }
  937. }
  938. lineItr++;
  939. }
  940. return true;
  941. }
  942. BOOL COTSImageProcess::ChangeFeaturelist(COTSFeatureList& a_listFeatures, COTSParticleList& a_listParticle)
  943. {
  944. for (auto pFeature : a_listFeatures)
  945. {
  946. COTSParticlePtr pParticle = COTSParticlePtr(new COTSParticle());
  947. pParticle->SetFeature(pFeature);
  948. a_listParticle.push_back(pParticle);
  949. }
  950. if ((int)a_listParticle.size() == 0)
  951. {
  952. return FALSE;
  953. }
  954. return TRUE;
  955. }
  956. BOOL COTSImageProcess::GetParticlesFromBinaryImage(CBSEImgPtr m_pBSEImg, CDoubleRange a_diameterRange, double a_pixelSize, COTSParticleList& listParticleOut)
  957. {
  958. int height, width;
  959. height = m_pBSEImg->GetHeight();
  960. width = m_pBSEImg->GetWidth();
  961. BYTE* pPixel = m_pBSEImg->GetImageDataPointer();
  962. Mat nobgMat = Mat(height, width, CV_8UC1, pPixel);
  963. Mat labels = Mat::zeros(nobgMat.size(), CV_32S);
  964. Mat stats, centroids;
  965. int number = connectedComponentsWithStats(nobgMat, labels, stats, centroids, 8, CV_32S);
  966. double rStart = a_diameterRange.GetStart()/2.0;
  967. double rEnd = a_diameterRange.GetEnd()/2.0;
  968. double areaStart = rStart* rStart * 3.14159;
  969. double areaEnd = rEnd* rEnd * 3.14159;
  970. //COTSParticleList listParticleOut;
  971. for (size_t i = 1; i < number; i++)
  972. {
  973. int center_x = centroids.at<double>(i, 0);
  974. int center_y = centroids.at<double>(i, 1);
  975. //矩形边框
  976. int x = stats.at<int>(i, CC_STAT_LEFT);
  977. int y = stats.at<int>(i, CC_STAT_TOP);
  978. int w = stats.at<int>(i, CC_STAT_WIDTH);
  979. int h = stats.at<int>(i, CC_STAT_HEIGHT);
  980. int area = stats.at<int>(i, CC_STAT_AREA);
  981. double actualArea = area * a_pixelSize * a_pixelSize;
  982. if (actualArea >= areaStart && actualArea < areaEnd)
  983. {
  984. Rect rectMax = Rect(x, y, w, h);
  985. Mat rectROI = labels(rectMax).clone();
  986. Mat imageROI = Mat::zeros(rectMax.size(), nobgMat.type());
  987. //exclude the point which intersect into this bounding box but is not in this group.
  988. int label = i;
  989. for (int row = 0; row < rectROI.rows; row++)
  990. {
  991. for (int col = 0; col < rectROI.cols; col++)
  992. {
  993. int v = rectROI.at<int>(row, col);
  994. if (v == label)
  995. {
  996. imageROI.at<uchar>(row, col) = 255;
  997. }
  998. }
  999. }
  1000. COTSParticleList roiParts;
  1001. if (!GetOneParticleFromROI(rectMax.x, rectMax.y, rectMax.width, rectMax.height, imageROI.data, roiParts))
  1002. {
  1003. continue;
  1004. }
  1005. if (roiParts.size() > 0)
  1006. {
  1007. COTSParticlePtr roiPart = roiParts[0];
  1008. roiPart->SetXRayPos(CPoint(center_x, center_y));
  1009. CRect r = CRect(x, y, x + w, y + h);
  1010. roiPart->SetParticleRect(r);
  1011. roiPart->SetActualArea(actualArea);
  1012. roiPart->SetPixelArea(area);
  1013. listParticleOut.push_back(roiPart);
  1014. }
  1015. }
  1016. }
  1017. return TRUE;
  1018. }
  1019. void COTSImageProcess::findPeakAndValley(const vector<int>& v, vector<int>& peakPositions, vector<int>& valleyPositions)
  1020. {
  1021. vector<int> diff_v(v.size() - 1, 0);
  1022. // 计算V的一阶差分和符号函数trend
  1023. for (vector<int>::size_type i = 0; i != diff_v.size(); i++)
  1024. {
  1025. if (v[i + 1] - v[i] > 0)
  1026. diff_v[i] = 1;
  1027. else if (v[i + 1] - v[i] < 0)
  1028. diff_v[i] = -1;
  1029. else
  1030. diff_v[i] = 0;
  1031. }
  1032. // 对Trend作了一个遍历
  1033. for (int i = diff_v.size() - 1; i >= 0; i--)
  1034. {
  1035. if (diff_v[i] == 0 && i == diff_v.size() - 1)
  1036. {
  1037. diff_v[i] = 1;
  1038. }
  1039. else if (diff_v[i] == 0)
  1040. {
  1041. if (diff_v[i + 1] >= 0)
  1042. diff_v[i] = 1;
  1043. else
  1044. diff_v[i] = -1;
  1045. }
  1046. }
  1047. for (vector<int>::size_type i = 0; i != diff_v.size() - 1; i++)
  1048. {
  1049. if (diff_v[i + 1] - diff_v[i] == -2)
  1050. peakPositions.push_back(i + 1);
  1051. if (diff_v[i + 1] - diff_v[i] == 2)
  1052. valleyPositions.push_back(i + 1);
  1053. }
  1054. }
  1055. std::vector<CIntRange> COTSImageProcess::GetValidGreyLevelRanges(CBSEImgPtr pBSEImg)
  1056. {
  1057. WORD originChartData[MAXBYTE];
  1058. //1. get chart data
  1059. pBSEImg->SetChartData();
  1060. memcpy(originChartData, pBSEImg->GetBSEChart(), sizeof(WORD) * MAXBYTE);
  1061. originChartData[0] = 0;
  1062. std::vector<int> greyChart;
  1063. for (int i = 0; i < MAXBYTE; i++)
  1064. {
  1065. greyChart.push_back((int)originChartData[i]);
  1066. }
  1067. std::vector<int> peaks;
  1068. std::vector<int> valleys;
  1069. findPeakAndValley(greyChart, peaks, valleys);
  1070. if ((peaks.size() - valleys.size()) == 1)
  1071. {
  1072. valleys.push_back(254);
  1073. }
  1074. std::map<int, MyGreyPeak*> peakMap;// hold all the peaks in this spectrum which are sorted by there area.
  1075. MyGreyPeak* lastRng = NULL;
  1076. MyGreyPeak* firstRng = NULL;
  1077. for (int i = 0; i < peaks.size(); i++)
  1078. {
  1079. MyGreyPeak* gr = new MyGreyPeak();
  1080. if (i == 0)
  1081. {
  1082. gr->start = 0;
  1083. }
  1084. else
  1085. {
  1086. gr->start = valleys[i - 1];
  1087. }
  1088. gr->end = valleys[i];
  1089. for (int j = gr->start; j < gr->end; j++)
  1090. {
  1091. gr->area += originChartData[j];
  1092. }
  1093. gr->peakPos = peaks[i];
  1094. gr->peakValue = originChartData[peaks[i]];
  1095. gr->valleyValue = valleys[i];
  1096. if (firstRng == NULL)
  1097. {
  1098. firstRng = gr;
  1099. lastRng = gr;
  1100. }
  1101. else
  1102. {
  1103. gr->preRng = lastRng;
  1104. lastRng->nextRng = gr;
  1105. }
  1106. lastRng = gr;
  1107. peakMap[gr->start] = gr;
  1108. }
  1109. int validGreyLevelNum = 0;//statistic the none zero grey level number
  1110. long partArea = 0;
  1111. for (int i = 0; i < 255; i++)
  1112. {
  1113. if (originChartData[i] != 0)
  1114. {
  1115. validGreyLevelNum += 1;
  1116. partArea += originChartData[i];
  1117. }
  1118. }
  1119. auto curRange = firstRng;
  1120. int maxPeakValue = curRange->peakValue;
  1121. while (curRange != NULL)
  1122. {
  1123. if (curRange->peakValue > maxPeakValue)
  1124. {
  1125. maxPeakValue = curRange->peakValue;
  1126. }
  1127. curRange = curRange->nextRng;
  1128. }
  1129. curRange = firstRng;
  1130. while (curRange != NULL)
  1131. {
  1132. curRange->maxPeakvalue = maxPeakValue;
  1133. curRange->allNoneZeroGreyLevelNum = validGreyLevelNum;
  1134. curRange = curRange->nextRng;
  1135. }
  1136. MyGreyPeak* firstRng1=NULL;
  1137. MyGreyPeak::MergeGreyPeaks(firstRng,firstRng1);
  1138. std::vector<CIntRange> ranges;
  1139. auto curRange1 = firstRng1;
  1140. while (curRange1 != NULL)
  1141. {
  1142. CIntRange pRange;
  1143. pRange.SetStart(curRange1->start);
  1144. pRange.SetEnd(curRange1->end);
  1145. ranges.push_back(pRange);
  1146. curRange1 = curRange1->nextRng;
  1147. }
  1148. return ranges;
  1149. }
  1150. BOOL COTSImageProcess::FindSeedsByGrayScale(CBSEImgPtr fieldImg, double a_PixelSize, std::vector<CPoint>& OutSeeds)
  1151. {
  1152. ASSERT(fieldImg);
  1153. ASSERT(m_imageProcessParam);
  1154. int nWidthImg = fieldImg->GetWidth();
  1155. int nHeightImg = fieldImg->GetHeight();
  1156. long nImgSize = nWidthImg * nHeightImg;
  1157. std::vector<CIntRange> validPeaks = GetValidGreyLevelRanges(fieldImg);
  1158. CRect r = CRect(0, 0, nWidthImg, nHeightImg);
  1159. CBSEImgPtr imgNoBGBinary = CBSEImgPtr(new CBSEImg(r));
  1160. for (int i = 0; i < validPeaks.size(); i++)
  1161. {
  1162. CIntRange pRange = validPeaks[i];
  1163. if (pRange.GetStart() == 0)//prevent the 0 grey pixels being recognized as particle.
  1164. {
  1165. pRange.SetStart(3);
  1166. }
  1167. long nNumParticle = 0;
  1168. GetSpecialGrayRangeImage(fieldImg, pRange, imgNoBGBinary, nNumParticle);
  1169. if (nNumParticle == 0)
  1170. {
  1171. continue;
  1172. }
  1173. else
  1174. {
  1175. // get the area image
  1176. Mat matImg = GetMatDataFromBseImg(imgNoBGBinary);
  1177. Mat labels = Mat::zeros(matImg.size(), CV_32S);
  1178. Mat stats, centroids;
  1179. int number = connectedComponentsWithStats(matImg, labels, stats, centroids, 8, CV_32S);
  1180. auto incsize = m_imageProcessParam->GetIncArea();
  1181. double rStart = incsize.GetStart() / 2.0;
  1182. double rEnd = incsize.GetEnd() / 2.0;
  1183. double pixelarea = a_PixelSize * a_PixelSize;
  1184. double areaStart = rStart* rStart * 3.14159*pixelarea;
  1185. double areaEnd = rEnd* rEnd * 3.14159*pixelarea;
  1186. /*double areaStart = 0;
  1187. double areaEnd = 1000;*/
  1188. for (size_t i = 1; i < number; i++)
  1189. {
  1190. int center_x = centroids.at<double>(i, 0);
  1191. int center_y = centroids.at<double>(i, 1);
  1192. int area = stats.at<int>(i, CC_STAT_AREA);
  1193. double actualArea = area * a_PixelSize * a_PixelSize;
  1194. if (actualArea >= areaStart && actualArea < areaEnd)
  1195. {
  1196. CPoint seed = CPoint(center_x, center_y);
  1197. OutSeeds.push_back(seed);
  1198. }
  1199. }
  1200. }
  1201. }
  1202. return true;
  1203. }
  1204. BOOL COTSImageProcess::GetMatrixPointsFromRect(CRect rect, double xrayStep, std::vector<CPoint>& matrixPs)
  1205. {
  1206. int colnum = ceil((double)rect.Width() / xrayStep + 0.5);
  1207. if (colnum % 2 == 0) colnum += 1;//let the number to be an odd number.Then we can make the middle point in the center of the particle exactly.
  1208. int rownum = ceil((double)rect.Height() / xrayStep + 0.5);
  1209. if (rownum % 2 == 0) rownum += 1;
  1210. CPoint theFirst = CPoint(rect.left - (colnum * xrayStep - rect.Width()) / 2 + xrayStep / 2, rect.top - (rownum * xrayStep - rect.Height()) / 2 + xrayStep / 2);
  1211. for (int i = 0; i < rownum; i++)
  1212. {
  1213. for (int j = 0; j < colnum; j++)
  1214. {
  1215. double x = (double)theFirst.x + (double)j * xrayStep;
  1216. double y = (double)theFirst.y + (double)i * xrayStep;
  1217. CPoint thePoint = CPoint(x, y);
  1218. if (rect.PtInRect(thePoint))
  1219. {
  1220. matrixPs.push_back(thePoint);
  1221. }
  1222. }
  1223. }
  1224. return true;
  1225. }
  1226. BOOL COTSImageProcess::CalcuParticleImagePropertes(COTSParticlePtr a_pOTSPart, double a_PixelSize)
  1227. {
  1228. //--------- convert this particle data to image data,construct an image only with this particle.------
  1229. const int nExpand_Size = 3;
  1230. const int nWhiteColor = 0;
  1231. const int nThickness = 1;
  1232. // lineType Type of the line
  1233. const int nLineType = 8;
  1234. // get rectangle of the particle
  1235. CRect rect = a_pOTSPart->GetParticleRect();
  1236. if (a_pOTSPart->GetActualArea() < 30 * a_PixelSize)// the particle is too small that openCV can't calculate a width value of it. Then we take the upright rect of the particle as it's minArea rect.
  1237. {
  1238. double w = 0, h = 0;
  1239. w = (double)rect.Width()*a_PixelSize;
  1240. h = (double)rect.Height()*a_PixelSize;
  1241. a_pOTSPart->SetDMax(MAX(w, h));
  1242. a_pOTSPart->SetDMin(MIN(w, h));
  1243. a_pOTSPart->SetDMean((w + h) / 2);
  1244. a_pOTSPart->SetFeretDiameter((w + h) / 2);
  1245. a_pOTSPart->SetDElong(MAX(w, h));
  1246. a_pOTSPart->SetPerimeter((w+h)*2);
  1247. a_pOTSPart->SetDPerp(MIN(w, h));
  1248. a_pOTSPart->SetDInscr(MIN(w, h));
  1249. a_pOTSPart->SetMinWidth(w);
  1250. a_pOTSPart->SetMinHeight(h);
  1251. return true;
  1252. }
  1253. // calculate the particle image data size, expand 3 pixel at the edge
  1254. Mat particleImage = Mat::zeros(rect.Height() + nExpand_Size , rect.Width() + nExpand_Size , CV_8U);
  1255. // get the segment list
  1256. COTSSegmentsList listSegment = a_pOTSPart->GetFeature()->GetSegmentsList();
  1257. for (auto pSegment : listSegment)
  1258. {
  1259. int nStart = pSegment->GetStart() - rect.left + nExpand_Size;
  1260. int nEnd = pSegment->GetStart() + pSegment->GetLength() - rect.left - 1 + nExpand_Size;
  1261. int nHeight = pSegment->GetHeight() - rect.top + nExpand_Size;
  1262. line(particleImage, Point(nStart, nHeight), Point(nEnd, nHeight), Scalar(nBlackColor), nThickness, nLineType);
  1263. }
  1264. //--------abstract the contour of the particle.
  1265. vector<vector<Point>>contours;
  1266. findContours(particleImage, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
  1267. if (contours.size()==0)// the particle is too odd that openCV can't find a contour of it. Then we take the upright rect of the particle as it's minArea rect.
  1268. {
  1269. double w = 0, h = 0;
  1270. w = (double)rect.Width()*a_PixelSize;
  1271. h = (double)rect.Height()*a_PixelSize;
  1272. a_pOTSPart->SetDMax(MAX(w, h));
  1273. a_pOTSPart->SetDMin(MIN(w, h));
  1274. a_pOTSPart->SetDMean((w + h) / 2);
  1275. a_pOTSPart->SetFeretDiameter((w + h) / 2);
  1276. a_pOTSPart->SetDElong(MAX(w, h));
  1277. a_pOTSPart->SetPerimeter((w + h) * 2);
  1278. a_pOTSPart->SetDPerp(MIN(w, h));
  1279. a_pOTSPart->SetDInscr(MIN(w, h));
  1280. a_pOTSPart->SetMinWidth(w);
  1281. a_pOTSPart->SetMinHeight(h);
  1282. return true;
  1283. }
  1284. int imaxcontour = 0, imax = 0;
  1285. for (unsigned int i = 0; i < contours.size(); i++) {
  1286. int itmp = contourArea(contours[i]);
  1287. if (imaxcontour < itmp) {
  1288. imax = i;
  1289. imaxcontour = itmp;
  1290. }
  1291. }
  1292. vector<Point > listEdge = contours[imax];
  1293. vector<vector<Point>>Outcontours;
  1294. Outcontours.push_back(listEdge);
  1295. //---------calculate the minimium rectangle
  1296. auto rRect = cv::minAreaRect(listEdge);
  1297. Point2f p[4];
  1298. rRect.points(p);
  1299. int D_MIN =getDistance(p[0], p[1]);
  1300. int D_MinRecLen = 0;//minareaRect's length(the longger side).
  1301. for (int j = 0; j <= 2; j++)
  1302. {
  1303. //line(cvContourImg, p[j], p[(j + 1) % 4], Scalar(100, 100, 0), 2);
  1304. int d = getDistance(p[j], p[j + 1]);
  1305. if (d < D_MIN)
  1306. {
  1307. D_MIN = d;
  1308. }
  1309. if (d > D_MinRecLen)
  1310. {
  1311. D_MinRecLen = d;
  1312. }
  1313. }
  1314. a_pOTSPart->SetDMin(D_MIN*a_PixelSize);
  1315. a_pOTSPart->SetMinWidth(rRect.size.width);
  1316. a_pOTSPart->SetMinHeight(rRect.size.height);
  1317. //----------calculate the perimeter
  1318. double d = arcLength(listEdge, true);
  1319. a_pOTSPart->SetPerimeter(d*a_PixelSize);
  1320. //-----------calculate the Max Diameter. Find the min enclosing circle first ,then find the two farthest circle connected point.
  1321. Point2f center; float radius;
  1322. minEnclosingCircle(listEdge, center, radius);
  1323. //circle(cvContourImg, center, radius, Scalar(100), 2);
  1324. std::vector <Point> outContour = listEdge;
  1325. std::vector <Point> rst;
  1326. for (unsigned int k = 0; k < outContour.size(); k++)
  1327. {
  1328. Point p = outContour[k];
  1329. double d = sqrt(pow((p.x - center.x), 2) + pow((p.y - center.y), 2));
  1330. if (fabs(d - radius) < 0.01)
  1331. {
  1332. rst.push_back(p);
  1333. }
  1334. }
  1335. double D_MAX = 0;
  1336. Point lineDmax[2];
  1337. for (unsigned int m = 0; m < rst.size(); m++)
  1338. {
  1339. Point p = rst[m];
  1340. for (unsigned int n = m + 1; n < rst.size(); n++)
  1341. {
  1342. Point p1 = rst[n];
  1343. double d = sqrt(powf((p.x - p1.x), 2) + powf((p.y - p1.y), 2));
  1344. if (d > D_MAX)
  1345. {
  1346. D_MAX = d;
  1347. lineDmax[0] = p;
  1348. lineDmax[1] = p1;
  1349. }
  1350. }
  1351. }
  1352. a_pOTSPart->SetDMax(D_MAX*a_PixelSize);
  1353. //--------calculate the D_PERP property using the D_MAX's two endpoints.
  1354. std::vector<Point> curve1;
  1355. std::vector<Point> curve2;
  1356. for (unsigned int i = 0; i < outContour.size(); i++)
  1357. {
  1358. Point pt = outContour[i];
  1359. bool start = false;
  1360. int clockwise = Side(lineDmax[0], lineDmax[1], pt);// devide these points into two group ,separate into the two sides.
  1361. if (clockwise > 0)
  1362. {
  1363. curve1.push_back(pt);
  1364. }
  1365. else
  1366. {
  1367. curve2.push_back(pt);
  1368. }
  1369. }
  1370. double d_perp1 = 0, d_perp2 = 0;
  1371. for (unsigned int i = 0; i < curve1.size(); i++)
  1372. {
  1373. double d = getDist_P2L(curve1[i], lineDmax[0], lineDmax[1]);
  1374. if (d > d_perp1)
  1375. {
  1376. d_perp1 = d;
  1377. }
  1378. }
  1379. for (unsigned int i = 0; i < curve2.size(); i++)
  1380. {
  1381. double d = getDist_P2L(curve2[i], lineDmax[0], lineDmax[1]);
  1382. if (d > d_perp2)
  1383. {
  1384. d_perp2 = d;
  1385. }
  1386. }
  1387. a_pOTSPart->SetDPerp((d_perp1 + d_perp2)*a_PixelSize);
  1388. //----------find the diameter of max inscribed circle
  1389. int r;
  1390. Point inscribeCirclecenter;
  1391. FindInnerCircleInContour(outContour, inscribeCirclecenter, r);
  1392. //--------------------------------------------------------calculate the xraypos !
  1393. CRect rec = a_pOTSPart->GetParticleRect();
  1394. a_pOTSPart->SetXRayPos(CPoint(inscribeCirclecenter.x - nExpand_Size + rec.left - 1, inscribeCirclecenter.y - nExpand_Size + rec.top - 1));
  1395. a_pOTSPart->SetDInscr(r * 2 * a_PixelSize);
  1396. //---------------calculate the image other caracater: length/width realArea/minRectangeArea etc. we can use these propertes to do forward process.
  1397. double minRectArea = D_MIN * D_MinRecLen*a_PixelSize*a_PixelSize;//最小外接矩形面积
  1398. double fillRatio = a_pOTSPart->GetActualArea() / minRectArea;//实际面积与最小外接矩形面积比,that's the fill rate.
  1399. double lengthWidthRatio;
  1400. lengthWidthRatio = (double)D_MinRecLen / D_MIN;//长宽比
  1401. //decide if this shape is a strip shape :if the lenthWidthRatio>2 then it is. if the lengthWidthRatio<2 and the areaRatio<0.5 then it is.
  1402. bool isStripShape = false;
  1403. double curveLength = 0;
  1404. double D_MEAN=0;
  1405. Moments mu;
  1406. mu = moments(listEdge, false);
  1407. int nx = mu.m10 / mu.m00;
  1408. int ny = mu.m01 / mu.m00;
  1409. //circle(cvcopyImg, Point(nx, ny), 1, (255), 1);
  1410. Point ptCenter = Point((int)nx, (int)ny);
  1411. if (pointPolygonTest(listEdge, ptCenter, false) != 1)// the center point doesn't contain in the contour, we think it as curve shape.
  1412. {
  1413. isStripShape = true;
  1414. }
  1415. /*if (lengthWidthRatio >= 2 )// in PartA software this is true,but IncA because of the GB definition the everage feret diameter is always the mean value of all the chord.
  1416. {
  1417. isStripShape = true;
  1418. }*/
  1419. if (fillRatio <= 0.4)// only when the fill rate is very low,we think it as a curve shape,then we choose the mean width as the feret diameter.
  1420. {
  1421. isStripShape = true;
  1422. }
  1423. if (isStripShape)
  1424. {
  1425. curveLength = a_pOTSPart->GetPerimeter()/2 - a_pOTSPart->GetDInscr()/2;// thinking this particle as a strip rectangle.the width is the max inscribe circle diameter/2.
  1426. if (curveLength < D_MAX)
  1427. {
  1428. curveLength = D_MAX;
  1429. }
  1430. if (curveLength < MIN_DOUBLE_VALUE || a_pOTSPart->GetActualArea()<MIN_DOUBLE_VALUE)
  1431. {
  1432. D_MEAN = 0;
  1433. }
  1434. else
  1435. {
  1436. D_MEAN = a_pOTSPart->GetActualArea() / curveLength;
  1437. }
  1438. a_pOTSPart->SetDMean(D_MEAN*a_PixelSize);
  1439. a_pOTSPart->SetFeretDiameter(D_MEAN*a_PixelSize);
  1440. a_pOTSPart->SetDElong (curveLength*a_PixelSize);
  1441. }
  1442. else//it's a ball shape particle
  1443. {
  1444. curveLength = D_MAX;
  1445. double ftd = 0, maxD = 0, minD = 0, dratio = 0;
  1446. GetParticleAverageChord(outContour, a_PixelSize, ftd);
  1447. a_pOTSPart->SetDMean(ftd);
  1448. a_pOTSPart->SetFeretDiameter(ftd);
  1449. a_pOTSPart->SetDElong(curveLength*a_PixelSize);
  1450. }
  1451. if (a_pOTSPart->GetAspectRatio() > 2)
  1452. {
  1453. double angle;
  1454. if (rRect.size.width > rRect.size.height) // w > h
  1455. {
  1456. angle = abs(rRect.angle);
  1457. }
  1458. else
  1459. {
  1460. angle = 90.0 + abs(rRect.angle);
  1461. }
  1462. a_pOTSPart->SetOrientation(angle);
  1463. }
  1464. else
  1465. {
  1466. a_pOTSPart->SetOrientation(0);
  1467. }
  1468. return true;
  1469. }
  1470. void COTSImageProcess::ImshowImage(CBSEImgPtr img)
  1471. {
  1472. BYTE* data = img->GetImageDataPointer();
  1473. //Mat cvImg;
  1474. cv::Size s;
  1475. s.width = img->GetImageSize().cx;
  1476. s.height = img->GetImageSize().cy;
  1477. Mat cvImg=Mat::zeros(s, CV_8U);
  1478. cvImg.data = data;
  1479. cv::imshow("dd", cvImg);
  1480. cv::waitKey();
  1481. }
  1482. void COTSImageProcess::ImshowChartData(CBSEImgPtr img)
  1483. {
  1484. img->SetChartData();
  1485. WORD* data = img->GetBSEChart();
  1486. //Mat cvImg;
  1487. cv::Size s;
  1488. s.width = 255;
  1489. s.height = 100;
  1490. Mat cvImg = Mat::zeros(s, CV_8U);
  1491. //cvImg.data = data;
  1492. WORD nBSEChart[MAXBYTE];
  1493. //1. get chart data
  1494. linearSmooth5(data, nBSEChart, MAXBYTE);
  1495. for (int i=1;i<255;i++)
  1496. {
  1497. line(cvImg, Point(i, 100-nBSEChart[i]), Point(i+1, 100-nBSEChart[i+1]), Scalar(nBlackColor), 1, 8);
  1498. }
  1499. cv::imshow("chart", cvImg);
  1500. cv::waitKey();
  1501. }
  1502. BOOL COTSImageProcess::MergeBigBoundaryParticles(COTSFieldDataList allFields,double pixelSize,int scanFieldSize, CSize ResolutionSize, COTSParticleList& mergedParts)
  1503. {
  1504. class BorderPart
  1505. {
  1506. typedef std::shared_ptr<BorderPart> CBorderPartPtr;
  1507. BorderPart(COTSParticlePtr p)
  1508. {
  1509. myPart = p;
  1510. headerParticle = NULL;
  1511. }
  1512. public:
  1513. COTSParticlePtr myPart;
  1514. COTSParticle* headerParticle;//used to merge particles ,if this particle has been merged then this pointer will point to the first particle of these merged particles or else it's NULL.
  1515. static std::vector<CBorderPartPtr> ConvertPartToBorderPart(COTSParticleList parts)
  1516. {
  1517. static std::map<COTSParticle*, CBorderPartPtr> allborderPart;
  1518. std::vector<CBorderPartPtr> borderParts;
  1519. for (auto p : parts)
  1520. {
  1521. if (allborderPart.find(p.get()) == allborderPart.end())
  1522. {
  1523. auto borderp = CBorderPartPtr(new BorderPart(p));
  1524. borderParts.push_back(borderp);
  1525. allborderPart[p.get()] = borderp;
  1526. }
  1527. else
  1528. {
  1529. borderParts.push_back(allborderPart[p.get()]);
  1530. }
  1531. }
  1532. return borderParts;
  1533. }
  1534. };
  1535. auto FldMgr = new CFieldMgr(scanFieldSize, ResolutionSize);
  1536. std::map<COTSParticle*, COTSParticleList> mapMergeParticles;//hold up all the boundary connected particles. the pair's first is also the member of these particles.
  1537. std::map<COTSParticle*, COTSSegmentsList> mapMergedSegments;//hold up all the segment's corresponding clone in the connected particles.
  1538. for (auto centerfld : allFields)
  1539. {
  1540. // find neighbor field on the left.
  1541. auto leftFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::LEFT);
  1542. if (leftFld != nullptr)
  1543. {
  1544. auto lParts = centerfld->GetLeftBorderedBigParticles();
  1545. auto rParts = leftFld->GetRightBorderedBigParticles();
  1546. auto leftParts = BorderPart::ConvertPartToBorderPart(lParts);
  1547. auto rightParts = BorderPart::ConvertPartToBorderPart(rParts);
  1548. for (auto leftp : leftParts)
  1549. {
  1550. for (auto rightp : rightParts)
  1551. {
  1552. if (leftp->myPart->IsConnected(rightp->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::LEFT))
  1553. {
  1554. if (leftp->headerParticle != NULL)
  1555. {
  1556. if (rightp->headerParticle == NULL)
  1557. {
  1558. rightp->headerParticle = leftp->headerParticle;
  1559. mapMergeParticles[leftp->headerParticle].push_back(rightp->myPart);
  1560. }
  1561. }
  1562. else
  1563. {
  1564. if (rightp->headerParticle != NULL)
  1565. {
  1566. leftp->headerParticle = rightp->myPart.get();
  1567. mapMergeParticles[rightp->myPart.get()].push_back(leftp->myPart);
  1568. }
  1569. else
  1570. {
  1571. leftp->headerParticle = leftp->myPart.get();
  1572. rightp->headerParticle = leftp->myPart.get();
  1573. mapMergeParticles[leftp->myPart.get()].push_back(rightp->myPart);
  1574. }
  1575. }
  1576. }
  1577. }
  1578. }
  1579. }
  1580. //find neighbor field on the upward
  1581. auto upFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::UP);
  1582. if (upFld != nullptr)
  1583. {
  1584. auto topBorderParts = centerfld->GetTopBorderedBigParticles();
  1585. auto bottomBorderParts = upFld->GetBottomBorderedBigParticles();
  1586. auto upParts = BorderPart::ConvertPartToBorderPart(topBorderParts);
  1587. auto downParts = BorderPart::ConvertPartToBorderPart(bottomBorderParts);
  1588. for (auto upprt : upParts)
  1589. {
  1590. for (auto downprt : downParts)
  1591. {
  1592. if (upprt->myPart->IsConnected(downprt->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::UP))
  1593. {
  1594. if (upprt->headerParticle != NULL)
  1595. {
  1596. if (downprt->headerParticle == NULL)
  1597. {
  1598. downprt->headerParticle = upprt->headerParticle;
  1599. mapMergeParticles[upprt->headerParticle].push_back(downprt->myPart);
  1600. }
  1601. }
  1602. else
  1603. {
  1604. if (downprt->headerParticle != NULL)
  1605. {
  1606. upprt->headerParticle = downprt->headerParticle;
  1607. mapMergeParticles[downprt->myPart.get()].push_back(upprt->myPart);
  1608. }
  1609. else
  1610. {
  1611. upprt->headerParticle = upprt->myPart.get();
  1612. downprt->headerParticle = upprt->myPart.get();
  1613. mapMergeParticles[upprt->myPart.get()].push_back(downprt->myPart);
  1614. }
  1615. }
  1616. }
  1617. }
  1618. }
  1619. }
  1620. //find neighbor field on the downward.
  1621. auto downFld = FldMgr->FindNeighborField(allFields, centerfld,SORTING_DIRECTION::DOWN);
  1622. if (downFld != nullptr)
  1623. {
  1624. auto bottomParts = centerfld->GetBottomBorderedBigParticles();
  1625. auto topParts = downFld->GetTopBorderedBigParticles();
  1626. auto downParts = BorderPart::ConvertPartToBorderPart(bottomParts);
  1627. auto upParts= BorderPart::ConvertPartToBorderPart(topParts);
  1628. for (auto downprt : downParts)
  1629. {
  1630. for (auto upprt : upParts)
  1631. {
  1632. if (downprt->myPart->IsConnected(upprt->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::DOWN))
  1633. {
  1634. if (downprt->headerParticle != NULL)
  1635. {
  1636. if (upprt->headerParticle == NULL)
  1637. {
  1638. upprt->headerParticle = downprt->headerParticle;
  1639. mapMergeParticles[downprt->headerParticle].push_back(upprt->myPart);
  1640. }
  1641. }
  1642. else
  1643. {
  1644. if (upprt->headerParticle != NULL)
  1645. {
  1646. downprt->headerParticle = upprt->headerParticle;
  1647. mapMergeParticles[upprt->headerParticle].push_back(downprt->myPart);
  1648. }
  1649. else
  1650. {
  1651. downprt->headerParticle = downprt->myPart.get();
  1652. upprt->headerParticle = downprt->myPart.get();
  1653. mapMergeParticles[downprt->myPart.get()].push_back(upprt->myPart);
  1654. }
  1655. }
  1656. }
  1657. }
  1658. }
  1659. }
  1660. //find neighbor field on the right.
  1661. auto rightFld = FldMgr->FindNeighborField(allFields, centerfld, SORTING_DIRECTION::RIGHT);
  1662. if (rightFld != nullptr)
  1663. {
  1664. auto rParts = centerfld->GetRightBorderedBigParticles();
  1665. auto lParts = rightFld->GetLeftBorderedBigParticles();
  1666. auto rightParts = BorderPart::ConvertPartToBorderPart(rParts);
  1667. auto leftParts = BorderPart::ConvertPartToBorderPart(lParts);
  1668. for (auto rightprt : rightParts)
  1669. {
  1670. for (auto leftprt : leftParts)
  1671. {
  1672. if (rightprt->myPart->IsConnected(leftprt->myPart.get(), centerfld->Width, centerfld->Height, (int)SORTING_DIRECTION::RIGHT))
  1673. {
  1674. if (rightprt->headerParticle != NULL)
  1675. {
  1676. if (leftprt->headerParticle == NULL)
  1677. {
  1678. leftprt->headerParticle = rightprt->headerParticle;
  1679. mapMergeParticles[rightprt->headerParticle].push_back(leftprt->myPart);
  1680. }
  1681. }
  1682. else
  1683. {
  1684. if (leftprt->headerParticle != NULL)
  1685. {
  1686. rightprt->headerParticle = leftprt->headerParticle;
  1687. mapMergeParticles[leftprt->headerParticle].push_back(rightprt->myPart);
  1688. }
  1689. else
  1690. {
  1691. rightprt->headerParticle = rightprt->myPart.get();
  1692. leftprt->headerParticle = rightprt->myPart.get();
  1693. mapMergeParticles[rightprt->myPart.get()].push_back(leftprt->myPart);
  1694. }
  1695. }
  1696. }
  1697. }
  1698. }
  1699. }
  1700. }
  1701. /*for (auto particle : mapMergeParticles)
  1702. {
  1703. }*/
  1704. static int partTagId;
  1705. for (auto pair : mapMergeParticles)
  1706. {
  1707. struct EleAreaPercentage
  1708. {
  1709. EleAreaPercentage(double p, CElementChemistryPtr e)
  1710. {
  1711. areaPercentage = p;
  1712. eleData = e;
  1713. }
  1714. double areaPercentage;
  1715. CElementChemistryPtr eleData;
  1716. };
  1717. auto newPart = COTSParticlePtr(new COTSParticle());
  1718. COTSSegmentsList newSegs;
  1719. auto p = pair.first;
  1720. newPart->SetSEMPos(p->GetSEMPos());
  1721. //firstly,we sum up all the merged particles's area and get the represent string.
  1722. std::string partsStr = std::to_string(p->GetFieldId()) + ":" + std::to_string(p->GetAnalysisId());
  1723. double allPartArea = p->GetActualArea();//Get the first particle's area.
  1724. for (auto other : pair.second)// Get the total area of all these particles for the use of ele calcu.
  1725. {
  1726. partsStr += "," + std::to_string(other->GetFieldId()) + ":" + std::to_string(other->GetAnalysisId());//Get the subparticles string such as "1:1,2:1" etc.
  1727. allPartArea += other->GetActualArea();//Get other particle's area
  1728. }
  1729. // calculate all the new segment's position.
  1730. std::vector <COTSParticle*> allSubParts;
  1731. allSubParts.push_back(p);
  1732. for (auto other : pair.second)// Get the total area of all these particles for the use of ele calcu.
  1733. {
  1734. allSubParts.push_back(other.get());
  1735. }
  1736. for (auto subp : allSubParts)
  1737. {
  1738. int fid = subp->GetFieldId();
  1739. CPoint myFldPos;
  1740. for (auto f : allFields)//find this particle's filed.
  1741. {
  1742. if (f->GetId() == fid)
  1743. {
  1744. myFldPos = f->GetPosition();
  1745. }
  1746. }
  1747. int fldWidth = allFields[0]->Width;
  1748. int fldHeight = allFields[0]->Height;
  1749. CPoint fldLeftUpPos = CPoint(myFldPos.x + fldWidth / 2 , myFldPos.y + fldHeight / 2 );
  1750. for (auto s : subp->GetFeature()->GetSegmentsList())
  1751. {
  1752. COTSSegmentPtr newseg = COTSSegmentPtr(new COTSSegment());
  1753. newseg->SetStart(s->GetStart() + fldLeftUpPos.x);
  1754. newseg->SetHeight((0 - s->GetHeight()) + fldLeftUpPos.y);//the coordinate system of segment in a field is different with the OTS coordinate system.OTS system's y axis is upward positive ,yet the field is downward positive.
  1755. newseg->SetLength(s->GetLength());
  1756. newSegs.push_back(newseg);
  1757. }
  1758. }
  1759. COTSFeaturePtr newFeature = COTSFeaturePtr(new COTSFeature());
  1760. newFeature->SetSegmentsList(newSegs);
  1761. newPart->SetFeature(newFeature);
  1762. newPart->CalCoverRectFromSegment();
  1763. //second, we get all the element data and their area percentage .
  1764. std::map<std::string, std::vector<EleAreaPercentage>> mapEleData;
  1765. CPosXrayPtr pXray1 = p->GetXrayInfo();
  1766. if (pXray1 != nullptr)
  1767. {
  1768. for (auto ele : pXray1->GetElementQuantifyData())
  1769. {
  1770. mapEleData[ele->GetName().GetBuffer()].push_back(EleAreaPercentage(p->GetActualArea() / allPartArea, ele));
  1771. }
  1772. }
  1773. for (auto other : pair.second)
  1774. {
  1775. auto otherXray = other->GetXrayInfo();
  1776. if (otherXray != nullptr)
  1777. {
  1778. for (auto eledata : otherXray->GetElementQuantifyData())
  1779. {
  1780. mapEleData[eledata->GetName().GetBuffer()].push_back(EleAreaPercentage(other->GetActualArea() / allPartArea, eledata));
  1781. }
  1782. }
  1783. }
  1784. // third,we calculate all the element's new percentage data and get a new element chemistry list.
  1785. CElementChemistriesList newCheList;
  1786. for (auto eledata : mapEleData)
  1787. {
  1788. CElementChemistryPtr newEleche = CElementChemistryPtr(new CElementChemistry());
  1789. newEleche->SetName(CString(eledata.first.c_str()));
  1790. double newPercentage = 0;
  1791. for (auto d : eledata.second)
  1792. {
  1793. newPercentage += d.areaPercentage * d.eleData->GetPercentage();
  1794. }
  1795. newEleche->SetPercentage(newPercentage);
  1796. newCheList.push_back(newEleche);
  1797. }
  1798. CPosXrayPtr xray(new CPosXray());
  1799. xray->SetElementQuantifyData(newCheList);
  1800. newPart->SetXrayInfo(xray);
  1801. newPart->SetConnectedParticlesSequentialString(partsStr);
  1802. newPart->SetActualArea(allPartArea);
  1803. partTagId++;
  1804. newPart->SetParticleId(partTagId);
  1805. newPart->SetAnalysisId(partTagId);
  1806. std::string name = p->GetClassifyName();
  1807. newPart->SetClassifyName(name);
  1808. newPart->SetColor(p->GetColor());
  1809. mergedParts.push_back(newPart);
  1810. }
  1811. return true;
  1812. }
  1813. }