File indexing completed on 2024-04-14 03:42:21
0001 /* 0002 SPDX-FileCopyrightText: 2020 Patrick Molenaar <pr_molenaar@hotmail.com> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 0006 Some code fragments were adapted from Peter Kirchgessner's FITS plugin 0007 SPDX-FileCopyrightText: Peter Kirchgessner <http://members.aol.com/pkirchg> 0008 */ 0009 0010 #include <math.h> 0011 #include <cmath> 0012 0013 #include "fits_debug.h" 0014 #include "fitsbahtinovdetector.h" 0015 #include "hough/houghline.h" 0016 #include "fitsdata.h" 0017 0018 #include <QElapsedTimer> 0019 #include <QtConcurrent> 0020 0021 //void FITSBahtinovDetector::configure(const QString &setting, const QVariant &value) 0022 //{ 0023 // if (!setting.compare("NUMBER_OF_AVERAGE_ROWS", Qt::CaseInsensitive)) 0024 // { 0025 // if (value.canConvert <int> ()) 0026 // { 0027 // NUMBER_OF_AVERAGE_ROWS = value.value <int> (); 0028 0029 // // Validate number of average rows value 0030 // if (NUMBER_OF_AVERAGE_ROWS % 2 == 0) 0031 // { 0032 // NUMBER_OF_AVERAGE_ROWS--; 0033 // qCInfo(KSTARS_FITS) << "Warning, number of rows must be an odd number, correcting number of rows to " 0034 // << NUMBER_OF_AVERAGE_ROWS; 0035 // } 0036 // // Rows must be a positive number! 0037 // if (NUMBER_OF_AVERAGE_ROWS < 1) 0038 // { 0039 // NUMBER_OF_AVERAGE_ROWS = 1; 0040 // qCInfo(KSTARS_FITS) << "Warning, number of rows must be positive correcting number of rows to " 0041 // << NUMBER_OF_AVERAGE_ROWS; 0042 // } 0043 // } 0044 // } 0045 //} 0046 0047 QFuture<bool> FITSBahtinovDetector::findSources(QRect const &boundary) 0048 { 0049 FITSImage::Statistic const &stats = m_ImageData->getStatistics(); 0050 switch (stats.dataType) 0051 { 0052 case TSHORT: 0053 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<int16_t>, boundary); 0054 0055 case TUSHORT: 0056 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<uint16_t>, boundary); 0057 0058 case TLONG: 0059 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<int32_t>, boundary); 0060 0061 case TULONG: 0062 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<uint32_t>, boundary); 0063 0064 case TFLOAT: 0065 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<float>, boundary); 0066 0067 case TLONGLONG: 0068 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<int64_t>, boundary); 0069 0070 case TDOUBLE: 0071 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<double>, boundary); 0072 0073 default: 0074 case TBYTE: 0075 return QtConcurrent::run(this, &FITSBahtinovDetector::findBahtinovStar<uint8_t>, boundary); 0076 0077 } 0078 } 0079 0080 template <typename T> 0081 bool FITSBahtinovDetector::findBahtinovStar(const QRect &boundary) 0082 { 0083 if (boundary.isEmpty()) 0084 return false; 0085 0086 QList<Edge*> starCenters; 0087 int subX = qMax(0, boundary.isNull() ? 0 : boundary.x()); 0088 int subY = qMax(0, boundary.isNull() ? 0 : boundary.y()); 0089 int subW = (boundary.isNull() ? m_ImageData->width() : boundary.width()); 0090 int subH = (boundary.isNull() ? m_ImageData->height() : boundary.height()); 0091 0092 int BBP = m_ImageData->getBytesPerPixel(); 0093 uint16_t dataWidth = m_ImageData->width(); 0094 0095 // #1 Find offsets 0096 uint32_t size = subW * subH; 0097 uint32_t offset = subX + subY * dataWidth; 0098 0099 // #2 Create new buffer 0100 auto * buffer = new uint8_t[size * BBP]; 0101 // If there is no offset, copy whole buffer in one go 0102 if (offset == 0) 0103 { 0104 memcpy(buffer, m_ImageData->getImageBuffer(), size * BBP); 0105 } 0106 else 0107 { 0108 uint8_t * dataPtr = buffer; 0109 const uint8_t * origDataPtr = m_ImageData->getImageBuffer(); 0110 // Copy data line by line 0111 for (int height = subY; height < (subY + subH); height++) 0112 { 0113 uint32_t lineOffset = (subX + height * dataWidth) * BBP; 0114 memcpy(dataPtr, origDataPtr + lineOffset, subW * BBP); 0115 dataPtr += (subW * BBP); 0116 } 0117 } 0118 0119 // #3 Create new FITSData to hold it 0120 FITSImage::Statistic stats; 0121 stats.width = subW; 0122 stats.height = subH; 0123 stats.dataType = m_ImageData->getStatistics().dataType; 0124 stats.bytesPerPixel = BBP; 0125 stats.samples_per_channel = size; 0126 FITSData* boundedImage = new FITSData(); 0127 boundedImage->restoreStatistics(stats); 0128 boundedImage->setProperty("dataType", stats.dataType); 0129 0130 // #4 Set image buffer 0131 boundedImage->setImageBuffer(buffer); 0132 0133 boundedImage->calculateStats(true); 0134 0135 QElapsedTimer timer1; 0136 0137 // Rotate image 180 degrees in steps of 1 degree 0138 QMap<int, BahtinovLineAverage> lineAveragesPerAngle; 0139 const int steps = 180; 0140 double radPerStep = M_PI / steps; 0141 0142 timer1.start(); 0143 0144 for (int angle = 0; angle < steps; angle++) 0145 { 0146 // TODO Apply multi threading to speed up calculation 0147 BahtinovLineAverage lineAverage = calculateMaxAverage<T>(boundedImage, angle); 0148 // Store line average in map 0149 lineAveragesPerAngle.insert(angle, lineAverage); 0150 } 0151 0152 qCDebug(KSTARS_FITS) << "Getting max average for all 180 rotations took" << timer1.elapsed() << "milliseconds"; 0153 0154 // Not needed anymore 0155 delete boundedImage; 0156 0157 // Calculate Bahtinov angles 0158 QVector<HoughLine*> bahtinov_angles; 0159 0160 // For all three Bahtinov angles 0161 for (int index1 = 0; index1 < 3; index1++) 0162 { 0163 double maxAverage = 0.0; 0164 double maxAngle = 0.0; 0165 int maxAverageOffset = 0; 0166 for (int angle = 0; angle < steps; angle++) 0167 { 0168 BahtinovLineAverage lineAverage = lineAveragesPerAngle[angle]; 0169 if (lineAverage.average > maxAverage) 0170 { 0171 maxAverage = lineAverage.average; 0172 maxAverageOffset = lineAverage.offset; 0173 maxAngle = angle; 0174 } 0175 } 0176 HoughLine* pHoughLine = new HoughLine(maxAngle * radPerStep, maxAverageOffset, subW, subH, maxAverage); 0177 if (pHoughLine != nullptr) 0178 { 0179 bahtinov_angles.append(pHoughLine); 0180 } 0181 0182 // Remove data around peak to prevent it from being detected again 0183 int minBahtinovAngleOffset = 18; 0184 for (int subAngle = maxAngle - minBahtinovAngleOffset; subAngle < maxAngle + minBahtinovAngleOffset; subAngle++) 0185 { 0186 int angleInRange = subAngle; 0187 if (angleInRange < 0) 0188 { 0189 angleInRange += 180; 0190 } 0191 if (angleInRange >= 180) 0192 { 0193 angleInRange -= 180; 0194 } 0195 lineAveragesPerAngle.remove(angleInRange); 0196 } 0197 } 0198 0199 // Proceed with focus offset calculation, but only when at least 3 lines have been detected 0200 QVector<HoughLine*> top3Lines; 0201 if (bahtinov_angles.size() >= 3) 0202 { 0203 HoughLine::getSortedTopThreeLines(bahtinov_angles, top3Lines); 0204 0205 // Debug output 0206 qCDebug(KSTARS_FITS) << "Sorted bahtinov angles:"; 0207 foreach (HoughLine* ln, top3Lines) 0208 { 0209 ln->printHoughLine(); 0210 } 0211 0212 // Determine intersection between outer lines 0213 HoughLine* oneLine = top3Lines[0]; 0214 HoughLine* otherLine = top3Lines[2]; 0215 QPointF intersection; 0216 HoughLine::IntersectResult result = oneLine->Intersect(*otherLine, intersection); 0217 if (result == HoughLine::INTERESECTING) 0218 { 0219 0220 qCDebug(KSTARS_FITS) << "Intersection: " << intersection.x() << ", " << intersection.y(); 0221 0222 // Determine offset between intersection and middle line 0223 HoughLine* midLine = top3Lines[1]; 0224 QPointF intersectionOnMidLine; 0225 double distance; 0226 if (midLine->DistancePointLine(intersection, intersectionOnMidLine, distance)) 0227 { 0228 qCDebug(KSTARS_FITS) << "Distance between intersection and midline is " << distance 0229 << " at mid line point " << intersectionOnMidLine.x() << ", " 0230 << intersectionOnMidLine.y(); 0231 0232 // Add star center to selected stars 0233 // Maximum Radius 0234 int maxR = qMin(subW - 1, subH - 1) / 2; 0235 BahtinovEdge* center = new BahtinovEdge(); 0236 center->width = maxR / 3; 0237 center->x = subX + intersection.x(); 0238 center->y = subY + intersection.y(); 0239 // Set distance value in HFR 0240 center->HFR = distance; 0241 0242 center->offset.setX(subX + intersectionOnMidLine.x()); 0243 center->offset.setY(subY + intersectionOnMidLine.y()); 0244 oneLine->Offset(subX, subY); 0245 midLine->Offset(subX, subY); 0246 otherLine->Offset(subX, subY); 0247 center->line.append(*oneLine); 0248 center->line.append(*midLine); 0249 center->line.append(*otherLine); 0250 starCenters.append(center); 0251 } 0252 else 0253 { 0254 qCWarning(KSTARS_FITS) << "Closest point does not fall within the line segment."; 0255 } 0256 } 0257 else 0258 { 0259 qCWarning(KSTARS_FITS) << "Lines are not intersecting (result: " << result << ")"; 0260 } 0261 } 0262 0263 // Clean up Bahtinov line array (of pointers) as they are no longer needed 0264 for (int index = 0; index < bahtinov_angles.size(); index++) 0265 { 0266 HoughLine* pLineAverage = bahtinov_angles[index]; 0267 if (pLineAverage != nullptr) 0268 { 0269 delete pLineAverage; 0270 } 0271 } 0272 bahtinov_angles.clear(); 0273 0274 // Clean up line averages array as they are no longer needed 0275 lineAveragesPerAngle.clear(); 0276 0277 top3Lines.clear(); 0278 0279 m_ImageData->setStarCenters(starCenters); 0280 0281 return true; 0282 } 0283 0284 template <typename T> 0285 BahtinovLineAverage FITSBahtinovDetector::calculateMaxAverage(const FITSData *data, int angle) 0286 { 0287 int BBP = data->getBytesPerPixel(); 0288 int size = data->getStatistics().samples_per_channel; 0289 int width = data->width(); 0290 int height = data->height(); 0291 int numChannels = data->channels(); 0292 0293 BahtinovLineAverage lineAverage; 0294 auto * rotimage = new T[size * BBP]; 0295 0296 rotateImage(data, angle, rotimage); 0297 0298 // Calculate average pixel value for each row 0299 auto * rotBuffer = reinterpret_cast<T *>(rotimage); 0300 0301 // printf("Angle;%d;Width;%d;Height;%d;Rows;%d;;RowSum;", angle, width, height, NUMBER_OF_AVERAGE_ROWS); 0302 0303 int NUMBER_OF_AVERAGE_ROWS = getValue("NUMBER_OF_AVERAGE_ROWS", 1).toInt(); 0304 if (NUMBER_OF_AVERAGE_ROWS % 2 == 0) 0305 { 0306 NUMBER_OF_AVERAGE_ROWS--; 0307 qCWarning(KSTARS_FITS) << "Warning, number of rows must be an odd number, correcting number of rows to " 0308 << NUMBER_OF_AVERAGE_ROWS; 0309 } 0310 // Rows must be a positive number! 0311 if (NUMBER_OF_AVERAGE_ROWS < 1) 0312 { 0313 NUMBER_OF_AVERAGE_ROWS = 1; 0314 qCWarning(KSTARS_FITS) << "Warning, number of rows must be positive correcting number of rows to " 0315 << NUMBER_OF_AVERAGE_ROWS; 0316 } 0317 0318 for (int y = 0; y < height; y++) 0319 { 0320 int yMin = y - ((NUMBER_OF_AVERAGE_ROWS - 1) / 2); 0321 int yMax = y + ((NUMBER_OF_AVERAGE_ROWS - 1) / 2); 0322 0323 unsigned long multiRowSum = 0; 0324 // Calculate average over multiple rows 0325 for (int y1 = yMin; y1 <= yMax; y1++) 0326 { 0327 int y2 = y1; 0328 if (y2 < 0) 0329 { 0330 y2 += height; 0331 } 0332 if (y2 >= height) 0333 { 0334 y2 -= height; 0335 } 0336 if (y2 < 0 || y2 >= height) 0337 { 0338 qCWarning(KSTARS_FITS) << "Y still out of bounds: 0 <=" << y2 << "<" << height; 0339 } 0340 0341 for (int x1 = 0; x1 < width; x1++) 0342 { 0343 int index = y2 * width + x1; 0344 unsigned long channelAverage = 0; 0345 for (int i = 0; i < numChannels; i++) 0346 { 0347 int offset = size * i; 0348 channelAverage += rotBuffer[index + offset]; 0349 } 0350 multiRowSum += qRound(channelAverage / static_cast<double>(numChannels)); 0351 } 0352 } 0353 // printf("%lu;", multiRowSum); 0354 0355 double average = multiRowSum / static_cast<double>(width * NUMBER_OF_AVERAGE_ROWS); 0356 if (average > lineAverage.average) 0357 { 0358 lineAverage.average = average; 0359 lineAverage.offset = y; 0360 } 0361 } 0362 // printf(";;MaxAverage;%.3f;MaxAverageIndex;%lu\r\n", lineAverage.average, lineAverage.offset); 0363 // fflush(stdout); 0364 0365 rotBuffer = nullptr; 0366 delete[] rotimage; 0367 rotimage = nullptr; 0368 0369 return lineAverage; 0370 } 0371 0372 /** Rotate an image by angle degrees. 0373 * verbose generates extra info on stdout. 0374 * return true if successful and rotated image. 0375 * @param angle The angle over which the image needs to be rotated 0376 * @param rotImage The image that needs to be rotated, also the rotated image return value 0377 */ 0378 template <typename T> 0379 bool FITSBahtinovDetector::rotateImage(const FITSData *data, int angle, T * rotImage) 0380 { 0381 int BBP = data->getBytesPerPixel(); 0382 int size = data->getStatistics().samples_per_channel; 0383 int width = data->width(); 0384 int height = data->height(); 0385 int numChannels = data->channels(); 0386 0387 int hx, hy; 0388 size_t bufferSize; 0389 0390 /* Check allocation buffer for rotated image */ 0391 if (rotImage == nullptr) 0392 { 0393 qWarning() << "No memory allocated for rotated image buffer!"; 0394 return false; 0395 } 0396 0397 while (angle < 0) 0398 { 0399 angle = angle + 360; 0400 } 0401 while (angle >= 360) 0402 { 0403 angle = angle - 360; 0404 } 0405 0406 hx = qFloor((width + 1) / 2.0); 0407 hy = qFloor((height + 1) / 2.0); 0408 0409 bufferSize = size * numChannels * BBP; 0410 memset(rotImage, 0, bufferSize); 0411 0412 auto * rotBuffer = reinterpret_cast<T *>(rotImage); 0413 auto * buffer = reinterpret_cast<const T *>(data->getImageBuffer()); 0414 0415 double innerCircleRadius = (0.5 * qSqrt(2.0) * qMin(hx, hy)); 0416 double angleInRad = angle * M_PI / 180.0; 0417 double sinAngle = qSin(angleInRad); 0418 double cosAngle = qCos(angleInRad); 0419 int leftEdge = qCeil(hx - innerCircleRadius); 0420 int rightEdge = qFloor(hx + innerCircleRadius); 0421 int topEdge = qCeil(hy - innerCircleRadius); 0422 int bottomEdge = qFloor(hy + innerCircleRadius); 0423 0424 for (int i = 0; i < numChannels; i++) 0425 { 0426 int offset = size * i; 0427 for (int x1 = leftEdge; x1 < rightEdge; x1++) 0428 { 0429 for (int y1 = topEdge; y1 < bottomEdge; y1++) 0430 { 0431 // translate point back to origin: 0432 double x2 = x1 - hx; 0433 double y2 = y1 - hy; 0434 0435 // rotate point 0436 double xnew = x2 * cosAngle - y2 * sinAngle; 0437 double ynew = x2 * sinAngle + y2 * cosAngle; 0438 0439 // translate point back: 0440 x2 = xnew + hx; 0441 y2 = ynew + hy; 0442 0443 int orgIndex = y1 * height + x1; 0444 int newIndex = qRound(y2) * height + qRound(x2); 0445 0446 if (newIndex >= 0 && newIndex < static_cast<int>(bufferSize)) 0447 { 0448 rotBuffer[newIndex + offset] = buffer[orgIndex + offset]; 0449 } // else index out of bounds, do not update pixel 0450 } 0451 } 0452 } 0453 0454 return true; 0455 }