SCRIPTS - added script to look at images offline to improve setttings for alc

This commit is contained in:
leonarski_f committed 2025-09-11 14:53:33 +02:00
1 parent 5f94ccd9d2
commit 206bfcf2db
1 file changed
+804
+804
View File
@@ -0,0 +1,804 @@
{
"cells": [
{
"cell_type": "code",
"id": "initial_id",
"metadata": {
"collapsed": true,
"ExecuteTime": {
"end_time": "2025-09-09T06:54:20.171830792Z",
"start_time": "2025-09-03T09:37:36.977231Z"
}
},
"source": [
"import cv2\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns"
],
"outputs": [],
"execution_count": 47
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-03T09:20:49.172543Z",
"start_time": "2025-09-03T09:20:49.169238Z"
}
},
"cell_type": "code",
"source": "threshold = 5",
"id": "19360e6dda22efb5",
"outputs": [],
"execution_count": 30
},
{
"metadata": {
"jupyter": {
"is_executing": true
}
},
"cell_type": "code",
"source": [
"bkg = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg1_50_50.jpg\")\n",
"curr_img = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/1002_0_1_50_50_curr_image.tiff\")"
],
"id": "9fac6cc561e58b70",
"outputs": [],
"execution_count": null
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-05T11:27:26.501898Z",
"start_time": "2025-09-05T11:27:26.496995Z"
}
},
"cell_type": "code",
"source": [
"def flatfield_correction(raw_image, flat_image, dark_image=None):\n",
" \"\"\"\n",
" Apply flat-field correction to remove uneven illumination\n",
"\n",
" Args:\n",
" raw_image: Your actual image with crystal\n",
" flat_image: Image of uniform illumination (no sample)\n",
" dark_image: Dark frame (camera with no light), optional\n",
" \"\"\"\n",
" # Convert to float for calculations\n",
" raw = raw_image.astype(np.float32)\n",
" flat = flat_image.astype(np.float32)\n",
"\n",
" if dark_image is not None:\n",
" dark = dark_image.astype(np.float32)\n",
" # Subtract dark frame from both\n",
" raw_corrected = raw - dark\n",
" flat_corrected = flat - dark\n",
" else:\n",
" raw_corrected = raw\n",
" flat_corrected = flat\n",
"\n",
" # Avoid division by zero\n",
" flat_corrected[flat_corrected == 0] = 1\n",
"\n",
" # Apply correction\n",
" mean_flat = np.mean(flat_corrected)\n",
" corrected = (raw_corrected / flat_corrected) * mean_flat\n",
"\n",
" # Convert back to original dtype\n",
" return np.clip(corrected, 0, 255).astype(raw_image.dtype)"
],
"id": "114345da7d31b16d",
"outputs": [],
"execution_count": 139
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-05T11:43:24.969687Z",
"start_time": "2025-09-05T11:43:24.954049Z"
}
},
"cell_type": "code",
"source": [
"def create_ring_mask(image, center=None, inner_radius=50, outer_radius=100):\n",
" \"\"\"\n",
" Create a mask to exclude a ring region\n",
"\n",
" Args:\n",
" image_shape: (height, width) of your image\n",
" center: (cx, cy) center of ring, None for image center\n",
" inner_radius: Inner radius of ring to mask out\n",
" outer_radius: Outer radius of ring to mask out\n",
"\n",
" Returns:\n",
" mask: Binary mask (0 = exclude ring, 255 = use for correction)\n",
" \"\"\"\n",
"\n",
" if len(image.shape) == 3:\n",
" h, w, _ = image.shape\n",
" else:\n",
" h, w = image.shape\n",
"\n",
"\n",
" if center is None:\n",
" center = (w // 2, h // 2)\n",
"\n",
" # Create coordinate arrays\n",
" cx, cy = center\n",
"\n",
" w = int(w)\n",
" h = int(h)\n",
"\n",
"\n",
" x = np.arange(w)\n",
" y = np.arange(h)\n",
" X, Y = np.meshgrid(x, y)\n",
"\n",
" # Calculate distances from center\n",
" distances = np.sqrt((X - cx)**2 + (Y - cy)**2)\n",
"\n",
"\n",
" # Create mask: exclude the ring region\n",
" mask = np.ones((h, w), dtype=np.uint8) * 255\n",
" ring_region = (distances >= inner_radius) & (distances <= outer_radius)\n",
" mask[ring_region] = 0\n",
"\n",
" return mask\n",
"\n",
"def flatfield_correction_with_mask(raw_image, flat_image, mask, dark_image=None):\n",
" \"\"\"\n",
" Apply flat-field correction with mask to exclude certain regions\n",
"\n",
" Args:\n",
" raw_image: Your camera image\n",
" flat_image: Flat-field reference\n",
" mask: Binary mask (0 = ignore, 255 = use)\n",
" dark_image: Dark frame (optional)\n",
" \"\"\"\n",
" # Convert to float\n",
" raw = raw_image.astype(np.float32)\n",
" flat = flat_image.astype(np.float32)\n",
"\n",
" if dark_image is not None:\n",
" dark = dark_image.astype(np.float32)\n",
" raw_corrected = raw - dark\n",
" flat_corrected = flat - dark\n",
" else:\n",
" raw_corrected = raw\n",
" flat_corrected = flat\n",
"\n",
" # Apply mask to flat field - only use unmasked regions for correction\n",
" mask_float = mask.astype(np.float32) / 255.0\n",
"\n",
" # Calculate mean only from unmasked regions\n",
" if len(flat_corrected.shape) == 3:\n",
" # Color image\n",
" mean_flat = []\n",
" for channel in range(flat_corrected.shape[2]):\n",
" masked_flat = flat_corrected[:,:,channel] * mask_float\n",
" valid_pixels = masked_flat[mask > 0]\n",
" mean_flat.append(np.mean(valid_pixels) if len(valid_pixels) > 0 else 1.0)\n",
" mean_flat = np.array(mean_flat).reshape(1, 1, -1)\n",
" else:\n",
" # Grayscale\n",
" masked_flat = flat_corrected * mask_float\n",
" valid_pixels = masked_flat[mask > 0]\n",
" mean_flat = np.mean(valid_pixels) if len(valid_pixels) > 0 else 1.0\n",
"\n",
" # Avoid division by zero\n",
" flat_corrected[flat_corrected <= 0] = 1\n",
"\n",
" # Apply correction\n",
" corrected = (raw_corrected / flat_corrected) * mean_flat\n",
"\n",
" # In masked regions, keep original values\n",
" mask_3d = mask_float\n",
" if len(corrected.shape) == 3:\n",
" mask_3d = np.stack([mask_float] * corrected.shape[2], axis=2)\n",
"\n",
" # Blend: use corrected where mask=1, original where mask=0\n",
" final_result = corrected * mask_3d + raw_corrected * (1 - mask_3d)\n",
"\n",
" return np.clip(final_result, 0, 255).astype(np.uint8)\n",
"\n",
"def create_interactive_ring_mask(image):\n",
" \"\"\"\n",
" Interactively create a ring mask by clicking on the image\n",
" Click center, then inner radius point, then outer radius point\n",
" \"\"\"\n",
" points = []\n",
"\n",
" def mouse_callback(event, x, y, flags, param):\n",
" if event == cv2.EVENT_LBUTTONDOWN:\n",
" points.append((x, y))\n",
" cv2.circle(image, (x, y), 3, (0, 255, 0), -1)\n",
" cv2.imshow('Select Ring', image)\n",
"\n",
" if len(points) == 1:\n",
" print(\"Click on inner radius of ring\")\n",
" elif len(points) == 2:\n",
" print(\"Click on outer radius of ring\")\n",
"\n",
" cv2.imshow('Select Ring', image)\n",
" print(\"Click on center of ring\")\n",
" cv2.setMouseCallback('Select Ring', mouse_callback)\n",
" cv2.waitKey(0)\n",
" cv2.destroyAllWindows()\n",
"\n",
" if len(points) == 3:\n",
" center = points[0]\n",
" inner_radius = np.sqrt((points[1][0] - center[0])**2 + (points[1][1] - center[1])**2)\n",
" outer_radius = np.sqrt((points[2][0] - center[0])**2 + (points[2][1] - center[1])**2)\n",
"\n",
" mask = create_ring_mask(image.shape, center, inner_radius, outer_radius)\n",
" return mask\n",
" else:\n",
" print(\"Need 3 points!\")\n",
" return None\n"
],
"id": "a1ae8c16d1cb31f7",
"outputs": [],
"execution_count": 154
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-08T14:23:06.438502Z",
"start_time": "2025-09-08T14:23:02.773985Z"
}
},
"cell_type": "code",
"source": [
"min_contour_area = 4\n",
"mask = create_ring_mask(curr_img, inner_radius=1000,outer_radius=2000)\n",
"diff_image = cv2.absdiff(curr_img, bkg)\n",
"diff_image = cv2.cvtColor(diff_image, cv2.COLOR_RGB2GRAY)\n",
"thresh_value, binary_image = cv2.threshold(diff_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n",
"thresh_value, thresh_normal = cv2.threshold(diff_image, 50, 255, cv2.THRESH_BINARY)\n",
"ff = flatfield_correction(curr_img, bkg, dark_image=None)\n",
"ff = cv2.cvtColor(ff, cv2.COLOR_RGB2GRAY)\n",
"ff_masked = flatfield_correction_with_mask(curr_img, bkg, mask, dark_image=None)\n",
"ff_masked = cv2.cvtColor(ff_masked, cv2.COLOR_RGB2GRAY)\n",
"thresh_value, thresh_mask = cv2.threshold(ff_masked, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)\n",
"kernel = np.ones((3,3), np.uint8)\n",
"thresh = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel)\n",
"thresh = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel)\n",
"\n",
"contours, _ = cv2.findContours(thresh_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n",
"valid_contours = [cnt for cnt in contours if cv2.contourArea(cnt) > min_contour_area]\n",
"\n",
"largest_contour = max(valid_contours, key=cv2.contourArea)\n",
"\n",
"# Step 6: Find leftmost point\n",
"leftmost_point = tuple(largest_contour[largest_contour[:,:,0].argmin()][0])\n",
"print(leftmost_point)\n",
"\n",
"leftmost_point = None\n",
"\n",
"for contour in contours:\n",
" # Loop through each point in the current contour\n",
" for point in contour:\n",
" x, y = point[0] # Point is a nested array [ [x, y] ]\n",
"\n",
" # Check if this is the leftmost point\n",
" if leftmost_point is None or x < leftmost_point[0]:\n",
" leftmost_point = (x, y)\n",
" plot_contour = contour\n",
"print(leftmost_point)\n",
"\n",
"c_img = cv2.drawContours(thresh, contours, -1, (0, 255, 0), 1)\n",
"plt.subplot(2, 2, 1)\n",
"plt.imshow(diff_image)\n",
"plt.subplot(2, 2, 2)\n",
"plt.imshow(thresh_normal)\n",
"#plt.subplot(2, 2, 4)\n",
"#plt.imshow(ff_masked)\n",
"\n",
"\n",
"result_image = cv2.cvtColor(thresh, cv2.COLOR_GRAY2BGR) if len(thresh.shape) == 2 else thresh.copy()\n",
"\n",
"# Draw contour\n",
"result_image = cv2.drawContours(result_image, [largest_contour], -1, (0, 255, 0), 2)\n",
"\n",
"# Mark leftmost point\n",
"result_image = cv2.circle(result_image, leftmost_point, 100, (0, 0, 255), -1)\n",
"result_image = cv2.putText(result_image, f'Leftmost: {leftmost_point}',\n",
" (leftmost_point[0] + 10, leftmost_point[1] - 10),\n",
" cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 0, 255), 2)\n",
"\n",
"print(f'Leftmost point: {leftmost_point}')\n",
"\n",
"plt.subplot(2, 2, 3)\n",
"plt.imshow(binary_image)\n"
],
"id": "a5363e8d3c71448",
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(np.int32(24), np.int32(979))\n",
"(np.int32(24), np.int32(979))\n",
"Leftmost point: (np.int32(24), np.int32(979))\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7f97e92802c0>"
]
},
"execution_count": 354,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"<Figure size 640x480 with 3 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgUAAAGiCAYAAAB3bbXGAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAbc9JREFUeJzt3Xl8HNWd7/3PqW51a1+8SLLANibse9iMJsATJr4WxGEg8TyvQHwTkiHhCbGZC04I8U0CZDIzTmBuZkIGyOQ1MzjzPElguPcCwxInjo3xAMKAgwMY8AVisMGWvMhSa+2tzvNHd5eqqlvGNpIly983r7a6q05VnWqkOr86WxlrrUVERESOeM54Z0BEREQmBgUFIiIiAigoEBERkTwFBSIiIgIoKBAREZE8BQUiIiICKCgQERGRPAUFIiIiAigoEBERkTwFBSIiIgJM8KDg7rvv5phjjqG8vJy5c+fy/PPPj3eWROQwoGuHyMGZsEHBAw88wNKlS7ntttv4/e9/z5lnnklbWxs7d+4c76yJyASma4fIwTMT9YFIc+fO5bzzzuMf//EfAXBdl5kzZ3LDDTfwrW99a5xzJyITla4dIgcvOt4ZKCWVSrFhwwaWLVvmLXMch3nz5tHe3l5ym2QySTKZ9D67rktXVxdTp07FGDPmeRY5ENZaent7aWlpwXEmbIXdYedArx26bsjh5FBcNyZkULB7926y2SxNTU2B5U1NTbzxxhslt1m+fDnf+973DkX2REbNtm3bOProo8c7G5PGgV47dN2Qw9FYXjcmZFBwMJYtW8bSpUu9zz09PcyaNYsLWUDUlI28ob/1RHcGcohkbJqneZyamprxzsoRbeTrxieJso/rhsg4yJDmaZ4Y0+vGhAwKpk2bRiQSobOzM7C8s7OT5ubmktvE43Hi8XjR8qgp23dQoDhAxotFVdSj7ECvHSNeN/iA64bIeMjfw47ldWNCNmbGYjHOOeccVq9e7S1zXZfVq1fT2to6jjkrwdpgbcPBbC8io+KwunaITEATsqYAYOnSpVxzzTWce+65nH/++fzDP/wD/f39fOlLXxrvrAV92IhNd4oio+qwuXaITEATNij47Gc/y65du7j11lvp6OjgrLPOYuXKlUUdiD6QtSM0EVgOfdvBeBxT5MgyatcOkSPQhJ2n4MNKJBLU1dXxcXOl2gZlwsnYNGvtw/T09FBbWzve2ZE877rBFbpuyISTsWnW8siYXjcmZJ8CEREROfQUFIiIiAigoEBERETyFBSIiIgIoKBARERE8hQUiIiICKCgQERERPIUFIiIiAigoEBERETyFBSIiIgIoKBARERE8hQUiIiICKCgQERERPIUFIiIiAigoEBERETyFBSIiIgIoKBARERE8hQUiIiICKCgQERERPIUFIiIiAigoEBERETyFBSIiIgIoKBARERE8hQUyH4wofe+z8aEkobWi4jIYSM63hmQw4Dx/vEv8H0MfwasGU5m7ZhlTURERo+CAskJF+y5hYDN/fS9BXKFvj8Z+fXeMl/iUrtWoCAiMuEoKDgSGRMslL2AIFTy52sICktDiwMxQi5wsL64IL+msKG3E1+QUWCL3oiIyDhQUHDEKZTsJZoDQs0EpujfUPJCVGCHtzXhct2AFyr4mxW8mgbr27lRDYKIyDhSUHDEMCOU7OG7fUJVAN4tfr6GIbQLL11+I4d8wT68ncG3ncmlHY4JSrU9iIjIeFBQMKmFA4HhZeH6gNxbX41BqQCisKEvTgDygUAoeBhe4Qsy8m8smHxaGwg+Cj/z+/GaG0RE5FAY9SGJt99+O8aYwOukk07y1g8NDbF48WKmTp1KdXU1CxcupLOzM7CPrVu3smDBAiorK2lsbOTmm28mk8kcZI6O1EKlRM2AMRhvsRlOYwzeikBAYLztAuu9VSYUSJjQi9D24eMZ/P95+/EfV8MbjwgT77ohcmQak5qCU089ld/97nfDB4kOH+amm27i8ccf58EHH6Suro4lS5bwmc98hmeeeQaAbDbLggULaG5u5tlnn2XHjh184QtfoKysjL/92789iNwcYYVKyXkD8m8L/xYFC8NrCbwd4bsLL7ahffo7DnrHD3ZLLGYx+dqG3I9Cc0S+6UF9DSa9iXXdEDkyjUlQEI1GaW5uLlre09PDv/zLv/DLX/6SP/3TPwXgvvvu4+STT+a5557jggsu4Le//S2vvfYav/vd72hqauKss87i+9//Prfccgu33347sVhsLLI8CZQu7APdBI0/Xb6avmgoYqngwL/AFqXxuhp4/QYCHQhyijoWmmCzgi+wMMOtDL7j5FcqOJi0dN0QGX9jMqPhm2++SUtLC8ceeyyLFi1i69atAGzYsIF0Os28efO8tCeddBKzZs2ivb0dgPb2dk4//XSampq8NG1tbSQSCTZt2jTiMZPJJIlEIvA6ovjv7r1mAn8TAQQCgkI1vrdxKG2+6r5QlUugapdANe9wOorXmRJ5KHWsQNNCLp3XqGD826hJYbLSdUNk/I16UDB37lxWrFjBypUruffee9myZQsXXXQRvb29dHR0EIvFqK+vD2zT1NRER0cHAB0dHYE/7ML6wrqRLF++nLq6Ou81c+bM0T2xicrfru/VDvgL3VC/gEBNQXFhHCz4h/dvRup/gP9n6DiFgj2/nTEG44SOOVJeneFl3rF951h6siU5XOm6ITIxjHrzwWWXXea9P+OMM5g7dy6zZ8/m3//936moqBjtw3mWLVvG0qVLvc+JRGLy/4GXmGvAECoszfA677MviAiWrYX6+xJNCqZUU8O+BCcxKLQSgM3vxgyPXCw0H0Cw7cAM78vYXB6CzRRGzQmThK4bIhPDmD8Qqb6+nhNOOIG33nqL5uZmUqkU3d3dgTSdnZ1eW2Jzc3NRr+LC51LtjQXxeJza2trAa1IL1A4M3+V7d/COAeME7+rDtQrDO/Pd5TsllrHvgMC/zhuNUGK7QB58tREj9IcIVEV4TQowXMMQPieZLHTdEBkfYx4U9PX18fbbbzNjxgzOOeccysrKWL16tbd+8+bNbN26ldbWVgBaW1t55ZVX2Llzp5dm1apV1NbWcsoppxxEDibhnaSvmn24zPcVpMYhWND6CneC/QKKqvJ9TQkjHbPoVVgffj/SNqE0gf4HJYOSUP5LBg0KDCaT8b9uiByZRr354Bvf+AaXX345s2fPZvv27dx2221EIhGuvvpq6urquPbaa1m6dClTpkyhtraWG264gdbWVi644AIA5s+fzymnnMLnP/957rjjDjo6OvjOd77D4sWLicfjB5GjyVRYlCisA3fM4TTG98PXVFDqzt6XPLh7U7zOFJoF8isKVfilRgj4C+9wusK2XmGfn9AI3/69VghfPqwdXmQDGWNiBYEfNAxTCibedUPkyDTqQcF7773H1VdfzZ49e5g+fToXXnghzz33HNOnTwfg7//+73Ech4ULF5JMJmlra+Oee+7xto9EIjz22GNcf/31tLa2UlVVxTXXXMNf/dVfjXZWDz++Qt0UFgRaA0J30OFgIFCb4N/v/jYNlMiMCa0YaV/WguMEA4bAUxfzfRYsGGOxgUcvQ4kHLWAsWFMqeJgogYECgv2l64bIxGCsnZw9tRKJBHV1dXzcXEnUlI13dj6k4N1/oB0+PKIgVEtQ3Ga/j4Dgg6rgTbBQ9g4TeMphiTSBQCD06xaeyrjw2bfcuta/MrcPOxwMWP/yUseYgDI2zVr7MD09PWrHnkC86wZXTILrhkw2GZtmLY+M6XVDzz44HJSsIQgHAsNVBvuuGcgX2uEAoFRA4JXvwSp6k/9ciCe9z4Ukdj/2ndvBcP4KQYB3LoVRCHY4GLG+5d7IA5tvePAFJBqVICJyUBQUTHRFowUKNQS+975C35tboFT/g/C+ipoGShTejgmW0773JrS9CU93zHDgUBRg2HBg4gsCCnkrlP/5N9ab+TA8Y6LNNyUwHDgoMBAROWBjPvpAPoTAsEOCve79NQX5NOFCOvCm1CiBwGdfHwGDN7wxdIgcO5zOU6pPXWA/oeOOdK6BEQq+c/UyFjg532df8IIvvYiI7DfVFExY4YCAUCE5vA4Kiwq1COECOBwIFC/zNwGY/A69G/LCbX6opiDwPrzMf+h8c4IB7wFHNrCf0IbWd6fvq6YwFCY88h3QqxnIpfM6HwaaJlRjICKyP1RTMBGF+gv4ilTfDbA/IPDVHJQcFljyFt632vjXFG9fqnOg/1iFqn7/eq+A9+3YFAICm5ubIBDtlKi9KKrRYHiaZO88CvkYXh74vorei4jISFRTMJF5ZWa4an24MC3qQzBSB0JftftwrYANFszhwt9/9+9fVqrA9Y7rq1GwJbYxXkV//gY/d3xb2M7rVOjbZ6CfgfVtU+DvZIivS0FwmYiI7JtqCiYU311wqU6B/rZ/f5JCslJ9BcKBRCCJE9yn8e1nnzfalhJVA8Xp/bUawYyGd5jj+Gcz9G8S+g78AVGgrwHF51sUtIiIyEgUFExgw2WeKSq8jVfwhdf5Cs5CdbpXrR5+4mCpW/rCbXWwSj5YyI5UsNtg0wHh96VOsLDr/LMdi4KhwkczXHtgTP7ZDqFaEBPMmwn+s++8iIiImg8mjtDdbuFnUdu6r+AM1A6U6P7vH/YXLsf3t3z0qu7D9e8jfDb59/6mAP9BvQ6AvnyEKxqMyQ1l9A1LxFooi+BWxXGr4rmgwDGYrAUXjGvBtRjXBddCxsW4uRcZF7JZyGax2SxkXaybVQdEkUnPEolC/bQMU5vSRMsskTJLcsAhkzEkBw2ppEM6aejvjZBJmXyz45FLQcFEYYJvih6BnA8ajL9qvOTQQorWFTr2jRwJjHDXH7hbz78PzS9gfW38gUcam3ATQ+jOP9COwXDaUv0QgExLHYmzp5NsLCcbNxgXTAYiKYgkIZK2OClLJGWJDro4aZfIUBYnmcUkM5hkGjOUxCZTMJTEplLYVDoXHIjIJGQ5/sxLine truncated
},
"metadata": {},
"output_type": "display_data"
}
],
"execution_count": 354
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-08T14:02:01.581369Z",
"start_time": "2025-09-08T14:01:56.983788Z"
}
},
"cell_type": "code",
"source": [
"# Load image and convert\n",
"img = cv2.imread('/home/leonarski_f/aaredaq/daq/src/aaredaq/991_0_1_50_50_curr_image.tiff', cv2.IMREAD_COLOR)\n",
"#img = cv2.imread('/home/leonarski_f/aaredaq/daq/src/aaredaq/0_280_100_50_curr_image_colour.jpg', cv2.IMREAD_COLOR)\n",
"\n",
"bkg = cv2.imread('/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg1_50_50.jpg')\n",
"# Convert to grayscale & enhance\n",
"gray_bkg = cv2.cvtColor(bkg, cv2.COLOR_BGR2GRAY)\n",
"gray = cv2.subtract(gray, gray_bkg)\n",
"gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n",
"gray = cv2.GaussianBlur(gray, (5, 5), 0)\n",
"gray = cv2.equalizeHist(gray)\n",
"\n",
"# Mask out camera-case ring\n",
"#gray = cv2.subtract(gray, gray_bkg)\n",
"mask = np.zeros_like(gray, dtype=np.uint8)\n",
"h, w = gray.shape\n",
"radius = min(h, w)//2 - 20\n",
"center = (w//2, h//2)\n",
"\n",
"norm=gray\n",
"print(radius)\n",
"\n",
"#norm = flatfield_correction_with_mask(gray, gray_bkg, mask)\n",
"\n",
"#radius = 500#radius#-factor\n",
"print(gray.shape)\n",
"# # gray = cv2.bitwise_and(gray, gray, mask=mask)\n",
"# # x, y = center[0] - radius, center[1] - radius\n",
"# gray_cropped = gray[y:y+2*radius, x:x+2*radius]\n",
"# #gray_bkg_cropped = gray[y:y+2*radius, x:x+2*radius]\n",
"# #cv2.circle(mask, center, radius, 93, -1)\n",
"# print(gray_cropped.shape)\n",
"#\n",
"# #masked = cv2.bitwise_and(gray, gray, mask=mask)\n",
"# #norm = gray_cropped\n",
"# masked=gray_cropped\n",
"# background= cv2.medianBlur(masked, 51)\n",
"#\n",
"# norm = cv2.subtract(masked, background)\n",
"#norm = cv2.absdiff(masked,gray_bkg)\n",
"kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (25,25))\n",
"tophat = cv2.morphologyEx(norm, cv2.MORPH_TOPHAT, kernel)\n",
"\n",
"enhanced = cv2.addWeighted(norm, 0.7, tophat, 0.3, 0)\n",
"# Adaptive thresholding\n",
"thresh = cv2.adaptiveThreshold(enhanced, 255, cv2.ADAPTIVE_THRESH_MEAN_C,\n",
" cv2.THRESH_BINARY_INV, 41, 2)\n",
"\n",
"# Find contours\n",
"contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n",
"\n",
"best_point = None\n",
"min_x = w + 1\n",
"\n",
"print(len(contours))\n",
"for cnt in contours:\n",
" x,y,wc,hc = cv2.boundingRect(cnt)\n",
"\n",
" # Skip very large contours (likely mask boundary)\n",
" if cv2.contourArea(cnt) > 0.2 * h * w:\n",
" print(\"Skipping contour with area\", cv2.contourArea(cnt))\n",
" continue\n",
"\n",
" # Compute distance of contour centroid from image center\n",
" M = cv2.moments(cnt)\n",
" if M[\"m00\"] != 0:\n",
" cx = int(M[\"m10\"]/M[\"m00\"])\n",
" cy = int(M[\"m01\"]/M[\"m00\"])\n",
" dist = np.sqrt((cx-center[0])**2 + (cy-center[1])**2)\n",
" # If centroid is near mask radius, skip (mask edge)\n",
" #if abs(dist - radius) < 60:\n",
" # print(\"Skipping contour near mask radius\")\n",
" # continue\n",
"\n",
" # Skip vertical pin\n",
" if hc > 3*wc:\n",
" print(\"Skipping vertical pin\")\n",
" continue\n",
"\n",
" leftmost = tuple(cnt[cnt[:, :, 0].argmin()][0])\n",
" if leftmost[0] < min_x:\n",
" min_x = leftmost[0]\n",
" best_point = leftmost\n",
"\n",
"# Draw result\n",
"if best_point is not None:\n",
" print(f\"best_point = {best_point}\")\n",
" #cv2.circle(thresh, best_point, 100, (0, 0, 255), -1)\n",
" print(radius)\n",
" corrected_point = (best_point[0]+factor, best_point[1] +factor)\n",
" print(f\"corrected_point {corrected_point}\")\n",
" cv2.circle(img, corrected_point, 100, (0, 0, 255), -1)\n",
" print(\"Loop tip (leftmost point):\", best_point)\n",
"else:\n",
" print(\"No suitable contour found\")\n",
"\n",
"plt.figure(figsize=(10,8))\n",
"plt.figure(figsize=(14,10))\n",
"plt.subplot(2,4,1); plt.title(\"Original\"); plt.imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))#; plt.axis(\"off\")\n",
"plt.subplot(2,4,2); plt.title(\"Cropped\"); plt.imshow(gray_cropped, cmap=\"gray\"); plt.axis(\"off\")\n",
"plt.subplot(2,4,3); plt.title(\"Background\"); plt.imshow(norm, cmap=\"gray\"); plt.axis(\"off\")\n",
"plt.subplot(2,4,4); plt.title(\"Top-hat\"); plt.imshow(tophat, cmap=\"gray\"); plt.axis(\"off\")\n",
"plt.subplot(2,4,5); plt.title(\"Enhanced\"); plt.imshow(enhanced, cmap=\"gray\"); plt.axis(\"off\")\n",
"plt.subplot(2,4,6); plt.title(\"Threshold\"); plt.imshow(thresh, cmap=\"gray\")#; plt.axis(\"off\")\n",
"plt.subplot(2,4,7); plt.title(\"image with contours\"); #plt.imshow(img)#; plt.axis(\"off\")\n",
"contour_vis = cv2.cvtColor(masked, cv2.COLOR_GRAY2BGR)\n",
"cv2.drawContours(contour_vis, contours, -1, (0,255,0), 1)\n",
"cv2.circle(contour_vis, best_point, 100, (0, 0, 255), -1)\n",
"plt.imshow(contour_vis)\n",
"plt.show()\n",
"\n",
"cv2.imwrite('masked.tiff', enhanced)\n",
"cv2.imwrite('thresh.tiff', thresh)"
],
"id": "d98c070865b11f36",
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1003\n",
"(2046, 2046)\n",
"10291\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping contour with area 2054203.5\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"Skipping vertical pin\n",
"best_point = (np.int32(0), np.int32(1421))\n",
"1003\n",
"corrected_point (np.int32(250), np.int32(1671))\n",
"Loop tip (leftmost point): (np.int32(0), np.int32(1421))\n"
]
},
{
"data": {
"text/plain": [
"<Figure size 1000x800 with 0 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<Figure size 1400x1000 with 7 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAABH0AAALRCAYAAADGGN8bAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsvXeUXNWVPbwr5+rKVV2dg7rViq2ABJIQEkk2JjobgxGDxwkcZjzO8xkbe/CMzTgMY/jBGINtjLFxwAZMFlFIAhRoZXUO1aGqK+f4vj96neNXpepWSwYj2W+vpaWuF+67L91377777CMTBEGABAkSJEiQIEGCBAkSJEiQIEGChL8ryN/uCkiQIEGCBAkSJEiQIEGCBAkSJEh48yGRPhIkSJAgQYIECRIkSJAgQYIECX+HkEgfCRIkSJAgQYIECRIkSJAgQYKEv0NIpI8ECRIkSJAgQYIECRIkSJAgQcLfISTSR4IECRIkSJAgQYIECRIkSJAg4e8QEukjQYIECRIkSJAgQYIECRIkSJDwdwiJ9JEgQYIECRIkSJAgQYIECRIkSPg7hET6SJAgQYIECRIkSJAgQYIECRIk/B1CIn0kSJAgQYIECRIkSJAgQYIECRL+DiGRPhLednzjG9+ATCY7pX3vu+8+yGQyDA0NvbmVEmFoaAgymQz33XffW3YMCRIkSBBDanckSDj90dzcjEsvvfTtrsbbBplMhm984xtvdzUkSJBwiti6dSuMRuPbXQ0JfwNIpI+EvwoHDx7ENddcg7q6Omg0Gni9Xnz4wx/GwYMH3+6qSZAg4e8A/f39+PjHP47W1lZotVqYzWasX78eP/rRj5BOp9/u6kmQIOE0AE0Aif+5XC5s3rwZjz/++NtdPQkSJEhgVLZVs/17/vnn3+6qnhJuvfVWPPzww293NSRUQPl2V0DCmYvf//73+NCHPgSbzYYbbrgBLS0tGBoawj333IPf/va3ePDBB3HVVVedsJx///d/x5e//OVTqsO1116LD37wg9BoNKe0vwQJEk5fPPbYY3jf+94HjUaDj3zkI1iyZAlyuRxefvllfOELX8DBgwdx9913v93VlCBBwmmCW265BS0tLRAEAVNTU7jvvvtwySWX4JFHHvmHVuRIkCDh9MEvfvGLst8///nP8fTTTx+3vKur629ZrTcNt956K9773vfiyiuvfLurIkEEifSRcEro7+/Htddei9bWVrz44otwOp287rOf/SzOPfdcXHvttejp6UFra2vVMpLJJAwGA5RKJZTKU3sUFQoFFArFKe0rQYKE0xeDg4P44Ac/iKamJmzbtg21tbW87sYbb0RfXx8ee+yxqvuWSiXkcjlotdq/VXUlSJBwGuCd73wnVq9ezb9vuOEGuN1u/OpXvzpjSR/qK0mQIOHvA9dcc03Z7507d+Lpp58+brkECW8mpPAuCaeE733ve0ilUrj77rvLCB8AcDgcuOuuu5BMJvHd734XwF98ew4dOoSrr74aVqsVGzZsKFsnRjqdxmc+8xk4HA6YTCZcfvnl8Pl8x8WPV/P0oRj7l19+GWvWrIFWq0Vrayt+/vOflx0jFArh3/7t37B06VIYjUaYzWa8853vxBtvvPEmXikJEiScCr773e8ikUjgnnvuKSN8CO3t7fjsZz8LYEYqfdNNN+GXv/wlFi9eDI1GgyeeeAIAsHfvXrzzne+E2WyG0WjEBRdcgJ07d5aVRe3Iiy++iI9//OOw2+0wm834yEc+gnA4XLYttS9PPfUUuru7odVqsWjRIvz+978/ro6RSASf+9zn0NDQAI1Gg/b2dvzXf/0XSqXScdtt3boVNTU1sFgsuO666xCJRP6ayydBggQAFosFOp2ubGLptttuw7p162C326HT6bBq1Sr89re/rbr//fffjzVr1kCv18NqtWLjxo146qmn5jzmz372MyiVSnzhC1/gZcFgENdeey3MZjO/42+88cZxvl3kr9Hf349LLrkEJpMJH/7whwHMkD+f//znuT3p7OzEbbfdBkEQeP+5vMAq+0/U9+rr68PWrVthsVhQU1OD66+/HqlUqmzfbDaLf/mXf4HT6eQ+2djY2JzXQYIECaeO+bzvQHn/p7OzE1qtFqtWrcKLL754Usfz+Xy48sorYTQa4XQ68W//9m8oFotl28yn7ZTJZEgmk/jZz37GYWpbt249pWsg4c2FpPSRcEp45JFH0NzcjHPPPbfq+o0bN6K5ufm4mfj3ve99WLBgAW699dbjGi4xtm7dit/85je49tprcfbZZ+OFF17Au971rnnXr6+vD+9973txww034LrrrsNPf/pTbN26FatWrcLixYsBAAMDA3j44Yfxvve9Dy0tLZiamsJdd92F8847D4cOHYLX65338SRIkPDm4pFHHkFrayvWrVs3r+23bduG3/zmN7jpppvgcDjQ3NyMgwcP4txzz4XZbMYXv/hFqFQq3HXXXdi0aRNeeOEFrF27tqyMm266CRaLBd/4xjdw9OhR3HnnnRgeHsbzzz9fRkz39vbiAx/4AD7xiU/guuuuw7333ov3ve99eOKJJ3DRRRcBAFKpFM477zz4fD58/OMfR2NjI1555RV85StfwcTEBH74wx8CAARBwBVXXIGXX34Zn/jEJ9DV1YU//OEPuO66696cCylBwj8QotEopqenIQgC/H4/br/9diQSibIZ9B/96Ee4/PLL8eEPfxi5XA4PPvgg3ve+9+HRRx8t62d885vfxDe+8Q2sW7cOt9xyC9RqNXbt2oVt27bh4osvrnr8u+++G5/4xCfw1a9+Fd/+9rcBzCgPL7vsMrz66qv45Cc/iYULF+KPf/zjrO94oVDAli1bsGHDBtx2223Q6/UQBAGXX345nnvuOdxwww3o7u7Gk08+iS984Qvw+Xz4wQ9+cMrX7P3vfz9aWlrwne98B3v27MFPfvITuFwu/Nd//Rdv89GPfhT3338/rr76aqxbtw7btm07qT6ZBAkS5o+Tfd9feOEF/PrXv8ZnPvMZaDQa3HHHHXjHO96BV199FUuWLDnh8YrFIrZs2YK1a9fitttuwzPPPIP//u//RltbGz75yU/ydvNpO3/xi1/gox/9KNasWYOPfexjAIC2trY38epIOGUIEiScJCKRiABAuOKKK+bc7vLLLxcACLFYTLj55psFAMKHPvSh47ajdYTdu3cLAITPfe5zZdtt3bpVACDcfPPNvOzee+8VAAiDg4O8rKmpSQAgvPjii7zM7/cLGo1G+PznP8/LMpmMUCwWy44xODgoaDQa4ZZbbilbBkC499575zxfCRIkvDmIRqPzamMIAAS5XC4cPHiwbPmVV14pqNVqob+/n5eNj48LJpNJ2LhxIy+jdmTVqlVCLpfj5d/97ncFAMIf//hHXkbty+9+97uy+tbW1gorVqzgZd/61rcEg8EgHDt2rKxOX/7ylwWFQiGMjIwIgiAIDz/8sABA+O53v8vbFAoF4dxzz5XaHQkS5gl6hyv/aTQa4b777ivbNpVKlf3O5XLCkiVLhPPPP5+X9fb2CnK5XLjqqquO6yeUSiX+u6mpSXjXu94lCIIg/OhHPxJkMpnwrW99q2z73/3udwIA4Yc//CEvKxaLwvnnn3/cO37dddcJAIQvf/nLZWVQO/Htb3+7bPl73/teQSaTCX19fYIgzN1fqew/Ud/rn/7pn8q2u+qqqwS73c6/9+3bJwAQPvWpT5Vtd/XVVx9XpgQJEk4eN954Y9k4aL7vuyAI3Na9/vrrvGx4eFjQarXCVVdddcJjU5sjHvcIgiCsWLFCWLVqVdmy+bSdgiAIBoNBuO666054bAl/W0jhXRJOGvF4HABgMpnm3I7Wx2IxXvaJT3zihOVTWManPvWpsuWf/vSn513HRYsWlamQnE4nOjs7MTAwwMs0Gg3k8plXoFgsIhgMwmg0orOzE3v27Jn3sSRIkPDmgtqME7UxYpx33nlYtGgR/y4Wi3jqqadw5ZVXlvmK1dbW4uqrr8bLL79c1jYBwMc+9jGoVCr+/clPfhJKpRJ//vOfy7bzer1lJvUUCrZ3715MTk4CAB566CGce+65sFqtmJ6e5n8XXnghisUiS6///Oc/Q6lUls2mKRSKk2rvJEiQMIMf//jHePrpp/H000/j/vvvx+bNm/HRj360LPxSp9Px3+FwGNFoFOeee27Zd//hhx9GqVTC17/+de4nECrD0YGZcNTPfvaz+K//+i/8+7//e9m6J554AiqVCv/8z//My+RyOW688cZZz0PcHgAz7YRCocBnPvOZsuWf//znIQjCX5WhrLJfdu655yIYDHL7SO1f5bE/97nPnfIxJUiQMDtO9n0/55xzsGrVKv7d2NiIK664Ak8++eRxIVqzoVo7IB4zAfNrOyWcvpDCuyScNGggRuTPbKhGDrW0tJyw/OHhYcjl8uO2bW9vn3cdGxsbj1tmtVrL/DlKpRJ+9KMf4Y477sDg4GBZw2i32+d9LAkSJLy5MJvNAE7cxohR2V4EAgGkUil0dnYet21XVxdKpRJGR0c53BMAFixYULad0WhEbW1tmWcYMNMWVQ78Ojo6AMx4ang8HvT29qKnp+c4zzOC3+8HMNPe1dbWwmg0lq2vVm8JEiTMjTVr1pQZOX/oQx/CihUrcNNNN+HSSy+FWq3Go48+im9/+9vYt28fstksbyt+p/v7+yGXy8uI5Nnwwgsv4LHHHsOXvvSlMh8fAr3jer2+bPlsfRqlUon6+vrjyvB6vccR4ZTdZ3h4+IT1nA2V/SWr1QpgZlBnNpu5T1YZoiG1URIkvDU42fe9su8CzPRJUqkUAoEAbDYbQqFQ2Xqn08mJcLRa7XF9lcoxE4B5tZ0STl9IpI+Ek0ZNTQ1qa2vR09Mz53Y9PT2oq6vjARxQzhK/lZgto5cg8hG69dZb8f/9f/8f/umf/gnf+ta3YLPZIJfL8bnPfe44o1UJEiT87WA2m+H1enHgwIF57/OLine truncated
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 350,
"metadata": {},
"output_type": "execute_result"
}
],
"execution_count": 350
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-08T11:21:42.726663Z",
"start_time": "2025-09-08T11:21:41.711795Z"
}
},
"cell_type": "code",
"source": [
"x#bkg_old = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_020925.jpg\")\n",
"#bkg_new = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_03092025.jpg\")\n",
"bkg_old = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg280_50_50_old.jpg\")\n",
"bkg_new = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg280_50_50.jpg\")\n",
"#$bkg_newest = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/bkg280.jpg\")\n",
"grey_old = cv2.cvtColor(bkg_old, cv2.COLOR_RGB2GRAY)\n",
"grey_new = cv2.cvtColor(bkg_new, cv2.COLOR_RGB2GRAY)\n",
"#grey_newest = cv2.cvtColor(bkg_newest, cv2.COLOR_RGB2GRAY)\n",
"diff_image = cv2.absdiff(grey_old, grey_new)\n",
"_, thresh_diff = cv2.threshold(diff_image, threshold, 255, cv2.THRESH_BINARY)\n",
"\n",
"\n",
"# Simple subtraction (may result in negative values)\n",
"diff_simple = cv2.subtract(bkg_new, bkg_old)\n",
"\n",
"# Weighted difference for motion detection\n",
"diff_weighted = cv2.addWeighted(bkg_new, 0.5, bkg_old, -0.5, 128)\n",
"\n",
"# Using numpy for difference\n",
"diff_numpy = np.abs(bkg_new.astype(np.int16) - bkg_old.astype(np.int16)).astype(np.uint8)\n",
"\n",
"plt.subplot(1,2,1)\n",
"plt.imshow(bkg_old)\n",
"plt.subplot(1,2,2)\n",
"plt.imshow(bkg_new)"
],
"id": "5c26fa11f85e6483",
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.image.AxesImage at 0x7f97bc9dc7a0>"
]
},
"execution_count": 275,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"text/plain": [
"<Figure size 640x480 with 2 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAESCAYAAAAWg2R/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs/W/Mtl12FoQfa9/DTKEwb+XfDBNL1Q8KNQpJ0WESMFEnjEpIjP0A2EhjGklMS6IjKo2m1D+xSf2gEhFjYuwXTdAPEoJJQ6kfmsBQtKaRVCFoSArRdwrUztv2F6ad51y/D3sdxzrWvq77ed+3fe7p+zzvvZ/cz3Vd57nP/WfttY691tpr7zMyM/GcntNzek7P6Tk9p+f0hqf1S92A5/ScntNzek7P6Tk9p69GelZ6ntNzek7P6Tk9p+f0oUjPSs9zek7P6Tk9p+f0nD4U6VnpeU7P6Tk9p+f0nJ7ThyI9Kz3P6Tk9p+f0nJ7Tc/pQpGel5zk9p+f0nJ7Tc3pOH4r0rPQ8p+f0nJ7Tc3pOz+lDkZ6Vnuf0nJ7Tc3pOz+k5fSjSs9LznJ7Tc3pOz+k5PacPRXpWep7Tc3pOz+k5Pafn9KFIH2il54//8T+Ov+/v+/vwNV/zNfj0pz+Nv/SX/tIvdZOe03N6Tq9BesaO5/ScntO99IFVev7kn/yT+PznP48/+kf/KP63/+1/w2/5Lb8Fn/vc5/ATP/ETv9RNe07P6Tl9gNMzdjyn5/ScHkvxQX3h6Kc//Wn8Y//YP4b//D//zwEA13Xh67/+6/GH/tAfwh/5I3/kl7h1z+k5PacPanrGjuf0nJ7TY+kjv9QNuJd+7ud+Dj/yIz+C7/zO79S1tRY++9nP4gtf+MLdZ7785S/jy1/+sn5f14Wf/MmfxK/5Nb8GEfHkbX5Oz+k5zZSZ+Omf/ml86lOfwlpfHafy+8WOZ9x4Ts/pg5eeEjs+kErP3/7bfxsvXrzAJz7xiXH9E5/4BP7KX/krd5/5nu/5Hvx7/96/99Vo3nN6Ts/pfaS/8Tf+Bv7ev/fv/arU9X6x4xk3ntNz+uCmp8COD6TS8wtJ3/md34nPf/7z+v2lL30Jv/E3/kYEAhGBBBAIAGnf76VztS+O6yzl3nP3793P/15uHlkPy9N/v9/vLysXADICC+zVow3yHwBy5H3UUg4AGTcF9+VAZACR2IuvqbpYy67pEv32r9S4Bp/DffKmXY0M/T5b7Ku/zBEs9yUrw6RbIkcD7j1zXuPvUfcj370vg5WeYNH6cb73HH0tM5FI/Kpf9atefWNeUXocN1bz7/uRUfjY37/nv/eXADLvylqJwSFrv8BE2WKdqv5xjHjp7+je3rTufK7+TyQyApE5+3+Td2MQqv+egxjhaEMJXvXrigRySVZvkSwHlkR9yy4BmbwPLFxIAtQYxRq3ZMsPvMhsYLt9ytpX198DPpzX3u27MLR+OCXuYmPmnjM5Rva9Z8Hura4F6+o8o9w7dXny8jOvJ8OOD6TS82t/7a/Fw8MDvvjFL47rX/ziF/HJT37y7jMf+9jH8LGPfez2RgQiVg1AD/d9pafFBniMMWLcO6/fv/cyxcE+D0S8ByknUNwDrJcqOnG/78wnChneUzRfBmxbrk/aOECynLT646a+pXve80BGbkUITpe1SVb1PkQauPD7JeXN27KxqIQ5Ehk7V6Se3lnvKE/93D1lEVLCdiWhMpXnUKTiACqNxV0A4/0UwDi9WVbeNq1zVRvvgpAIcHL/Oba8fs/13Mz8GJ2eKr1f7HgMNyLimOAfUf8PMpGXw655Nn7PM8+N7HQ7mqHu8FuN97vR+L3jRAwdZNLgHnbETduEBWlGSAPXRoAwaQqb7KIxp8vqcjtP6Ipjyi4r8NBP1+e1Sy15zpoTsgrskhpbWimCfiWssdOasfrOu8QivxLqzCYVGagpnMoZU14PjHgv38EyIhubHkk3/JGJqCWmoNKk8b3RWXe+4/fQ2x+pOgaJHsfYX2z6QO7e+uhHP4pv+qZvwg/+4A/q2nVd+MEf/EF85jOf+QWW2rb9Heh41zswIXvPiYJ9znr3SnyMERxMCojf7c/z3nsO76Uc+4fA3eeX/UVsmFjg/W0px9p/iL6OCGAFEKt+15Nl2WFh318BrIVYu1zlrXJRZasu1Vt/KxDrAVhrrwuvKOGtz6jPZfWvKAuT7e2+e1+jns110BNh/Wh66c/HxvK/1/GM4KQ069QEEHOSuMe1mQfHO396AYP77/ktRqmPM/FXMT0NdgD3ZH+C9E7UkU9K9GRqeXJSbR0AL1Vbc88jOPIeFZ73xmMYecdzeIQfS3GhbEhuiv/XOuvj8wuL2LCWZEz4QBmPJSxYi/Jt+LNKLtfGh4lvq9rS8o1F3Ol2ImIbPuqDYc3AP6heyveKjSEoTBn4d9Nn0qxkM2IrcQe+djmlAqoMqN739QcIx2sAX8ozJx/K6HJ+ag4bjEwed57vgsaDs57TMn2i9IH09ADA5z//eXzrt34rfttv+234x//xfxz/6X/6n+Jnf/Zn8S//y//yL7hMxyhq9mY/4CT0Lx7CN7Olz0KHRSjvpz8lNI3jGgaz3mjkR77RG/8du2XDyDxcLjcsF5uxg8/Gaetvb0z3WhWBdk3o99mV0G/RCg3A25uzLZ7Isn1ckgJALi9qzC5sz7KbgUCui4XXc+qJGxxdWFV1iWaBB+s3qSp/YbWVSoWqustYWY6h4kx3MZ/XrJspb1K7rWlVsTyfLG2IG2TMLTSlYS4O3MrI0wHTLya9Gux4CTozxyOeHD7uWCNL93AItAjE8FqOSUZZTmU0brBiivLj+PBu10788HtrLToxR7qdDPuHeNjLKrcOF5gi1vQIUFYW0NK78y64cr8KH0SVWbmIbb/luQsNjmiXQyQQyXvl94l9bQmf+tkteFFX28fcdsSWKtLgQhoudXIMEQ2bgDWvXEbXxrmXYYbujQGcXuMb8+be/KK+Pm7Qe18033mGQ0Du1f0U6QOr9Pze3/t78bf+1t/Cd33Xd+Htt9/Gb/2tvxXf//3ffxOg+K4piyFE8AN1gOP3qQh5nvP+u9d987skKl+SjfXVvCpmPsHpfQFXdJn3fjvQeP1qYxCgSoDiDm0Eeh4NFHBXthClwE7AAwBYaodPzgNKBxhPREjSt8pyKxl2i5Akj4wVlxVEwTanSymFcwjpCfQ97fRSF0HmyF/Nn4r4BK27SMDnTDHqhwmId7jqHlvfZEnFLrw3ADql6lR8v/rplWHH+5F1e+Ix+Bj3bh50DeqOLL80zXxD+XgEMx6byO4aTdGGhmQoTJFx3pMisHmActxSsR9ObG9JmhKxjaImXgBaixCmDPlnPw8hIXA+Qp+b34ZlbZgx0qfvhzpLQizlAy61ZscsXWY4xOj/lPN1gwPqe6ZwlzjScGRKjtHgNIjejYe49BUVP/meZP6O4uR4rUbdPMe2+/1Dw/wqpA/sOT2/2PTOO+/grbfeQmAvc9zMLLiZ1uF37ui7R77HByqOXIxHyfOG58lbkHs35ebdlZ79X9j3eSsUKGm41TkCWOnXzzpOc8+B8lAJpHXW9aFl3U4uQa1ANONK7Mvon3d+ZV+/s2yT5zinXc+wS1cDnsq0JdPsx0eZ2SXcxPYgt9eIeazu9Gd53eJ6znifRJYXDGCw4kkF/yIHwaDGHd7FuyU3IuYTVyYyL3zpS1/Cxz/+8Xct6YOQiBuLS6/vMTln3aWbdNjD8/bIpPTyCeu+DLwfJee9Kj2MsUkpA8BiHAr7Kk9CXZESA/PG7P93TMmk68AIKkG6a8aL4Mbl0PGuZTacwUUzx5B5f8rHUB30IVqYTDdOmNqfhDuTSeYdlfX3rTb5xfpjvY4vVb/i9wDJfALy+uLAiHf7FLIdnuEbXlYHT/w/6GVlnN1+t5R54Xoi7PjAenpeWWoZPtzRJ7z773MSTtwO4btX60+pmLOU3I271VVeDl6/MG/PVHq0UypqWSaaEhug9mcIbQRP/d9od1s+rjypBfRJuzLGMgRWRh3zCg2iNRLP5CgLz3LNMu1ODI/GOf13/1OKW4zc7sXxHg80y3ndIWUdz2Y1kV1Ra4/lLn6XW7szDmbnhOVLdg3ct8mp835A6mZA8v4QvT7pvcn5Y0+p70bMfW/ugrk/cdxXhpyzb1j/JThwz1B62ZLFxB9K8G5nDnnsDQZTblvqvegsT+4Q5yB/dpk47u+l69zxOGzLAPaus3d2YeaNy2T3HKXOl8QHEbhxo1vpyt216aBOru01wVbufKQmT9zuF70Z8SzMwdVKlOWU3lHtXlR8IsbmiHdLY8kLGN/vJlOMX5akd8ZU9EcetuGrCBZvvNIz2DRudw/dTnQtUDY1vq86h7wUl49abjSc23vv/bPdvoYcU6FgNeaOZnuWfg9x0qO0yiIaHKQMcWY+3ekFgtaco3/150td1jaCrBSNMRIx67ckQAE9H8yz+yCAQIIxQBm9pu7LRLuftOI2HdKJWWvzPQlY+zpYpvLvWW9bbGZtciltfD84kuz4GCgM/fxU1k/CmNpzblOzHO9f2Tm/v/Yaz07vow+PUH7ovjdGxyMKz8vqUFlHmeP7e1B66kIpJgMkpMB4a6OWpLBmjFxiIdbmmigsyFHmnrSJJC3f3ZFwBg+vWSCEoepElCyvcX/fSmg3Fr1VSAAPVk5KTjsmJqp6tkMHU5glMOm5l+Ye6jnlLnptJcsWydCly1xSt3tXVdNii20WZiTgWOx8Zd/YMtHuJfx7s8OL9DLDaZqIdwuZPMjLs1nvngwynho63nilB7iB+zt3mSapbwl/QJsVPNy9NufZfDjB4GzFI+D1qAUnwPI+cALusrQm7FdHEyYA6cpZZKELlZZ+JmaeoBDYivMpFASQIEi2hdHVDsjtsq0fjwXfuWUKtkE/CzyioSdO5hgD3yDU9bAtx1gSjGjiDK8Q2wSBWMfjXAUw0V4eI9MtwE3A6mvkwwZ1jFbTks3m0/eLine truncated
},
"metadata": {},
"output_type": "display_data"
}
],
"execution_count": 275
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-08T08:13:17.148953Z",
"start_time": "2025-09-08T08:13:17.046977Z"
}
},
"cell_type": "code",
"source": [
"bkg_old_np = bkg_old.astype(np.float32)\n",
"bkg_old_np_normalized = bkg_old.astype(np.float32) / 255.0\n",
"bkg_new_np = bkg_old.astype(np.float32)\n",
"bkg_new_np_normalized = bkg_old.astype(np.float32) / 255.0\n",
"diff_np=bkg_old_np-bkg_new_np\n",
"diff_np_normalized=bkg_old_np_normalized-bkg_new_np_normalized"
],
"id": "95a0ac65d0013e8",
"outputs": [],
"execution_count": 196
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-08T08:14:30.738697Z",
"start_time": "2025-09-08T08:14:30.671257Z"
}
},
"cell_type": "code",
"source": [
"print(\"mean: \", np.mean(grey_old.astype(np.float32)-grey_new.astype(np.float32)))\n",
"print(\"std: \", np.std(grey_old.astype(np.float32)-grey_new.astype(np.float32)))\n",
"diff_abs = np.abs(grey_old.astype(np.float32)-grey_new.astype(np.float32))\n",
"print(np.percentile(diff_abs, 99))"
],
"id": "83aae96e9df6c807",
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mean: -6.184164\n",
"std: 5.93317\n",
"21.0\n"
]
}
],
"execution_count": 200
},
{
"metadata": {},
"cell_type": "code",
"outputs": [],
"execution_count": null,
"source": [
"threshold_old = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_020925.jpg\")\n",
"threshold_new = cv2.imread(\"/home/leonarski_f/aaredaq/daq/src/aaredaq/alc_test_images/bkg280_03092025.jpg\")"
],
"id": "258a3d679516f15c"
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-05T13:33:03.863166Z",
"start_time": "2025-09-05T13:33:03.839938Z"
}
},
"cell_type": "code",
"source": [
"import numpy as np\n",
"from scipy import ndimage\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# Example showing the difference\n",
"image = np.random.rand(50, 50) * 100\n",
"# Add a bright spot\n",
"image[20:30, 35:45] = 200\n",
"\n",
"# Without thresholding - all pixels contribute\n",
"com_no_threshold = ndimage.center_of_mass(image)\n",
"print(f\"No threshold COM: ({com_no_threshold[1]:.2f}, {com_no_threshold[0]:.2f})\")\n",
"\n",
"# With thresholding - only bright pixels contribute\n",
"threshold = 150\n",
"thresholded = image.copy()\n",
"thresholded[thresholded < threshold] = 0\n",
"com_with_threshold = ndimage.center_of_mass(thresholded)\n",
"print(f\"With threshold COM: ({com_with_threshold[1]:.2f}, {com_with_threshold[0]:.2f})\")\n"
],
"id": "d70801a606268e76",
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"No threshold COM: (26.11, 24.33)\n",
"With threshold COM: (39.50, 24.50)\n"
]
}
],
"execution_count": 189
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-05T13:33:18.254639Z",
"start_time": "2025-09-05T13:33:18.248074Z"
}
},
"cell_type": "code",
"source": [
"def center_of_mass_with_threshold(image, threshold):\n",
" \"\"\"\n",
" Apply threshold before center of mass calculation\n",
" \"\"\"\n",
" # Create binary mask\n",
" binary = (image >= threshold).astype(float)\n",
"\n",
" # Calculate center of mass on thresholded image\n",
" com = ndimage.center_of_mass(binary)\n",
"\n",
" # Return as (x, y) coordinates\n",
" return (com[1], com[0]) if com[0] is not np.nan else None\n"
],
"id": "1ba37750a7b66b23",
"outputs": [],
"execution_count": 190
},
{
"metadata": {
"ExecuteTime": {
"end_time": "2025-09-05T13:33:27.236834Z",
"start_time": "2025-09-05T13:33:27.231958Z"
}
},
"cell_type": "code",
"source": [
"def weighted_center_of_mass_threshold(image, threshold):\n",
" \"\"\"\n",
" Keep original intensities above threshold, zero below\n",
" \"\"\"\n",
" weighted_image = image.copy()\n",
" weighted_image[weighted_image < threshold] = 0\n",
"\n",
" com = ndimage.center_of_mass(weighted_image)\n",
" return (com[1], com[0]) if com[0] is not np.nan else None\n"
],
"id": "59452341d64bc800",
"outputs": [],
"execution_count": 191
},
{
"metadata": {},
"cell_type": "code",
"outputs": [],
"execution_count": null,
"source": "",
"id": "8cf8041e736e5a43"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}