@@ -131,98 +131,86 @@ namespace {
131131 template <typename T> PixelDirection& operator =(T&&) = delete ;
132132 };
133133
134- void findStarShape (const CGrayBitmap& bitmap, CStar& star, const double backgroundNoiseLevel)
134+ void findStarShape (const CGrayBitmap& bitmap, CStar& star, double backgroundNoiseLevel)
135135 {
136136 constexpr int AngleResolution = 10 ; // Test the axis every 10 degrees.
137+ constexpr double GradRadFactor = 3.14159265358979323846 / 180.0 ;
138+
137139 std::array<CStarAxisInfo, 360 / AngleResolution> starAxes;
140+ std::array<int , 4 > xcoords, ycoords;
138141
139- // Preallocate the vector for the inner loop.
140- PIXELDISPATCHVECTOR vPixels;
141- vPixels.reserve (10 );
142+ const auto StarAxisRadius = [&starAxes, AngleResolution](const int angle) -> double
143+ {
144+ const int index = ((angle + 360 ) % 360 ) / AngleResolution;
145+ return starAxes[index].m_fRadius ;
146+ };
147+
148+ std::ranges::transform (std::views::iota (0 , static_cast <int >(starAxes.size ())), starAxes.begin (), [](const int n) { return CStarAxisInfo{ 0.0 , 0.0 , n * AngleResolution }; });
142149
143150 const int width = bitmap.Width ();
144151 const int height = bitmap.Height ();
152+ const double bitmapMultiplier = bitmap.GetMultiplier ();
153+ backgroundNoiseLevel *= bitmapMultiplier;
145154
146- for (int lAngle = 0 ; lAngle < 360 ; lAngle += AngleResolution )
155+ for (CStarAxisInfo& axisInfo : starAxes )
147156 {
148- double fSquareSum = 0.0 ;
149- double fSum = 0.0 ;
157+ double squareSum = 0.0 ;
150158
151159 for (double fPos = 0.0 ; fPos <= star.m_fMeanRadius * 2.0 ; fPos += 0.10 )
152160 {
153- constexpr double GradRadFactor = 3.14159265358979323846 / 180.0 ;
154- const double fX = star.m_fX + std::cos (lAngle * GradRadFactor) * fPos ;
155- const double fY = star.m_fY + std::sin (lAngle * GradRadFactor) * fPos ;
156- double fLuminance = 0 ;
157-
158- // Compute luminance at fX, fY
159- vPixels.resize (0 );
160- ComputePixelDispatch (QPointF (fX , fY ), vPixels);
161+ const double fX = star.m_fX + std::cos (axisInfo.m_lAngle * GradRadFactor) * fPos ;
162+ const double fY = star.m_fY + std::sin (axisInfo.m_lAngle * GradRadFactor) * fPos ;
163+ double luminanceOfPixel = 0 ;
161164
162- for (const CPixelDispatch& pixel : vPixels)
165+ std::array<double , 4 > proportions = ComputeAll4PixelDispatches (QPointF{ fX , fY }, xcoords, ycoords);
166+
167+ for (const int n : { 0 , 1 , 2 , 3 })
163168 {
164- if (pixel.m_lX < 0 || pixel.m_lX >= width || pixel.m_lY < 0 || pixel.m_lY >= height)
165- continue ;
166-
167- double pixelBrightness;
168- bitmap.GetPixel (static_cast <size_t >(pixel.m_lX ), static_cast <size_t >(pixel.m_lY ), pixelBrightness);
169- pixelBrightness = std::max (pixelBrightness - backgroundNoiseLevel, 0.0 ); // Exclude negative values.
170- fLuminance += pixelBrightness * pixel.m_fPercentage ;
169+ if (const size_t x = xcoords[n], y = ycoords[n]; x >= 0 && x < width && y >= 0 && y < height)
170+ {
171+ luminanceOfPixel += proportions[n] * std::max (bitmap.getUncheckedValue (x, y) - backgroundNoiseLevel, 0.0 );
172+ }
171173 }
172- fSquareSum += fPos * fPos * fLuminance ;
173- fSum += fLuminance ;
174- // fNrValues += fLuminance * 2;
174+ squareSum += fPos * fPos * luminanceOfPixel;
175+ axisInfo.m_fSum += luminanceOfPixel;
175176 }
176177
177- const double fStdDev = fSum > 0.0 ? std::sqrt (fSquareSum / fSum ) : 0.0 ;
178- CStarAxisInfo ai{ .m_fRadius = fStdDev * 1.5 , .m_fSum = fSum , .m_lAngle = lAngle };
179-
180- starAxes[lAngle / AngleResolution] = std::move (ai);
178+ axisInfo.m_fRadius = axisInfo.m_fSum > 0.0 ? 1.5 * std::sqrt (squareSum / axisInfo.m_fSum ) : 0.0 ;
179+ axisInfo.m_fSum /= bitmapMultiplier; // Correct the sum, because bitmap.getUncheckedValue(x, y) did not divide by the bitmap-multiplier.
181180 }
182181
183- const auto StarAxisRadius = [&starAxes](const int angle) -> double
184- {
185- const int index = ((angle + 360 ) % 360 ) / AngleResolution;
186- return starAxes[index].m_fRadius ;
187- };
188-
189182 //
190183 // Do a search over the first 180 degrees calculating the diameter, setting majorAxis to the largest found.
191184 //
192185 double majorAxis{ 0.0 };
193- double majorAxisAngle{ 0.0 } ;
194- for (int i = 0 ; i < starAxes. size () / 2 ; ++i )
186+ int majorAxisAngle = 0 ;
187+ for (int angle = 0 ; angle < 180 ; angle += AngleResolution )
195188 {
196- auto a{ starAxes[i].m_fRadius };
197- auto b{ starAxes[i + (starAxes.size () / 2 )].m_fRadius };
198- auto diameter{ a + b };
189+ const double diameter = StarAxisRadius (angle) + StarAxisRadius (angle + 180 );
199190 if (diameter > majorAxis)
200191 {
201192 majorAxis = diameter;
202- if (a >= b) majorAxisAngle = starAxes[i].m_lAngle ;
203- else majorAxisAngle = starAxes[i + (starAxes.size () / 2 )].m_lAngle ;
193+ majorAxisAngle = angle;
204194 }
205195 }
206196
207- star.m_fMajorAxisAngle = majorAxisAngle;
208- auto a{ StarAxisRadius (majorAxisAngle) };
209- auto b{ StarAxisRadius (majorAxisAngle + 180 ) };
210- star.m_fLargeMajorAxis = std::max (a, b);
211- star.m_fSmallMajorAxis = std::min (a, b);
197+ double a = StarAxisRadius (majorAxisAngle);
198+ double b = StarAxisRadius (majorAxisAngle + 180 );
199+ star.m_fMajorAxisAngle = a >= b ? majorAxisAngle : (majorAxisAngle + 180 + 360 ) % 360 ;
200+ std::tie (star.m_fSmallMajorAxis , star.m_fLargeMajorAxis ) = std::minmax (a, b);
212201 a = StarAxisRadius (majorAxisAngle + 90 );
213202 b = StarAxisRadius (majorAxisAngle + 270 );
214- star.m_fLargeMinorAxis = std::max (a, b);
215- star.m_fSmallMinorAxis = std::min (a, b);
203+ std::tie (star.m_fSmallMinorAxis , star.m_fLargeMinorAxis ) = std::minmax (a, b);
216204
217- double minorAxis{ StarAxisRadius (majorAxisAngle + 90 ) + StarAxisRadius (majorAxisAngle + 270 ) } ;
205+ const double minorAxis = a + b ;
218206
219207 //
220208 // Compute eccentricity - definition of the ecccentricity formula says to use the semi-minor and
221209 // semi-major axis values, but using the minor and major axis values gives the same result.
222210 //
223211 if (0 != majorAxis) // Just to avoid a division by zero, should never happen.
224212 {
225- star.eccentricity = std::sqrt (1.0 - ((minorAxis * minorAxis)/ (majorAxis * majorAxis)));
213+ star.eccentricity = std::sqrt (1.0 - ((minorAxis * minorAxis) / (majorAxis * majorAxis)));
226214 }
227215 }
228216}
0 commit comments