Обновить
-1

Пользователь

0,1
Рейтинг
1
Подписчики
Отправить сообщение

сеточки уходят в рандомо-генерацию мусора. вот и весь 'психоз'.

Мыши тоже, помню, плакали и кололись :)

Мне кажется абсолютно логично было бы вернуть в этой ситуации одно перечисление с тремя состояниями. Взаимоисключение в этом случае было бы невозможно, а 'лишние' значения ловились бы элементарно приемником (что лучше всегда и делать)

Думается, что все претензии из-за неумения ставить задачу нейронкам. Не нравятся полотна кода? Просите сети делать максимально компактный код. Получается неоптимально по скорости? Просите делать оптимально. Это не 'высосано из пальца'. Это реальные рекомендации по сотням чат-сессий.

Ну и да, актуальные LLM'ки фундаментально недетерминированные. Ожидать от них противного не стоит.

В любом случае буквально всегда есть вариант: не нравится - не используйте.

думается что это всего лишь вопрос конкретных реализаций а не самого принципа.

Я вот одного не понимаю. Чего они рвутся с максом на айфоны. Сервера не локализованы, всё идет через инфраструктуру США. Должны же радоваться по идее, но нет...

Лучше бы падишах. Но осел тоже норм.

По поводу такси - видимо, придётся вводить стеклянные перегородки между водителем и пассажирами...

Зачем такие сложности? Вверху был кадр из Кин-дза-дза.

По-хорошему можно заодно и колеса блокировать. Для уверенности, что точно никуда не поедешь.

То сейчас мы вышли на приличное такое плато роста производительности.

Даже на CPU это не так. А на GPU это прям совсем не так.

Абсолютно так. Я вот не понимаю массового нытья. Почему-то все всех с пеной у рта убеждают, что ИИшечка - прям отстой. Какие-то постоянные хохмочки, лулзы, подколки. Блин, все взрослые люди вроде бы. Вроде все сами в состоянии зайти и за буквально полчаса-час, даже абсолютно бесплатно, определить для себя: есть польза или нет. Ну вот какой практический смысл всех убеждать и отговаривать? Не нравится - не юзайте! Никто же не заставляет. Вроде.

К самой ИИшечке претензии тоже сильно надуманные. Такое чувство, что все пишут прямо идеальный код, абсолютно без ошибок, полностью покрытый тестами, максимально поддерживаемый, без лапши и т п. В общем, прямо призываю очень внимательно посмотреться в зеркало, прежде чем огульно критиковать чаты.

Всё уже когда-то было :) И, да, и это пройдет.

Скрытый текст

под эту модель 'Qwen3.6-27B' сколько нужно минимально видеопамяти для более-менее комфортной работы?

Все восторги обычно разбиваются об вопрос «а откуда это известно» (в смысле что «отлично»)? Как Вы можете быть уверены, что он не правдоподобную муму гонит (особенно с учётом того, что у Вас знаний отличить муму от правильного ответа очевидно недостаточно)

Откуда известно? Взяли глазами и посмотрели на результаты :) Программисты вместе с докторами. На множестве примеров. До этого кода было с пяток попыток запилить своё руками, или найти готовое, с весьма увы спорными результатами. Что должно получится известно и просто: контуры изображений должны максимально точно совпасть по всему набору, в результатах нет рокет саенса. Но сами сделать это увы долго не могли. Чатик вот справился. Знаний для решения недостаточно. Но постановку знаем хорошо и результаты можем оценить.

Это только один пример, можно сказать что уникальный, но он - один из множества. Я так то вообще из чатов (использую 3-4) не вылезаю фактически. Большую часть нудной текучки отдал чатам. Сам, конечно, тоже пишу куски. Результаты буквально превосходные.

Однако на 99% уверен, что специалист в этом языке Вам сходу парочку проблемных мест укажет.

Скорее нет :) Так то я как бы и сам специалист, почти 25 лет кодинга на этом вот языке и несколько проектов с нуля, несколько миллионов строк кода.

Я же написал в самом начале. Это код алгоритма регистрации. Гуглится. Ок, если сложно, вот ссылка:

https://3dqlab.stanford.edu/image-registration/

На входе пачка КТ изображений. Пациент во время скана бывает дергается, и это проблема для работы последующих алгоритмов (для примера: удаление костей по маске или по диапазону HU, расчет перфузионных карт CBV/CBF/MTT и т п). Нужна предобработка, что код и делает. Делает, как писал, отлично. Сомневающимся стучите в личку вышлю пачку скриншотов. Сюда тянуть не вижу смысла.

Вы так не беспокойтесь и не волнуйтесь, всё очень плотно тестируется :)

Не вопрос. Для примера, процедура регистрации (гуглится что это). почти 100% нейрокод, работает отлично:

procedure RegisterSeriesToFirstWordHU_RigidPyr_InPlace(const Frames: TArray<PDAWord>; // Frames[i] length = W*H
 W, H: integer; RescaleSlope, RescaleIntercept: double; FillHU: integer = -1024; MaxShiftPx: integer = 0;
 UseSubpixel: boolean = True; UseEdges: boolean = True; RoiX: integer = 0; RoiY: integer = 0;
 RoiW: integer = 0; RoiH: integer = 0; MaxAngleDeg: double = 5.0; AngleStepDeg: double = 1.0);

procedure RegisterSeriesToFirstWordFill_RigidPyr_InPlace(const Frames: TArray<PDAWord>;
 W, H: integer; RescaleSlope, RescaleIntercept: double; FillValueWord: word; MaxShiftPx: integer = 0;
 UseSubpixel: boolean = True; UseEdges: boolean = True; RoiX: integer = 0; RoiY: integer = 0;
 RoiW: integer = 0; RoiH: integer = 0; MaxAngleDeg: double = 5.0; AngleStepDeg: double = 1.0);

implementation

uses
 ProgressFormUnit,
 System.Threading,
 System.Diagnostics;

const
 EPS = 1e-12;

type
 TComplex = record
  Re, Im: double;
 end;

 TRigid = record
  AngleRad: double;
  Dx, Dy: double; // apply to Mov after rotation
  Peak: double;
 end;

function CConj(const A: TComplex): TComplex; inline;
begin
 Result.Re := A.Re;
 Result.Im := -A.Im;
end;

function CMul(const A, B: TComplex): TComplex; inline;
begin
 Result.Re := A.Re * B.Re - A.Im * B.Im;
 Result.Im := A.Re * B.Im + A.Im * B.Re;
end;

procedure CScale(var A: TComplex; S: double); inline;
begin
 A.Re := A.Re * S;
 A.Im := A.Im * S;
end;

function NextPow2(N: integer): integer;
begin
 Result := 1;
 while Result < N do Result := Result shl 1;
end;

procedure FFT1D(var Data: TArray<TComplex>; Invert: boolean);
var
 N, I, J, Bit, Len, HalfLen, K: integer;
 Ang, WRe, WIm, TmpRe, WlenRe, WlenIm: double;
 URe, UIm, VRe, VIm: double;
 V: TComplex;
begin
 N := Length(Data);
 if N <= 1 then Exit;

 // bit reversal
 J := 0;
 for I := 1 to N - 1 do
 begin
  Bit := N shr 1;
  while (J and Bit) <> 0 do
  begin
   J := J xor Bit;
   Bit := Bit shr 1;
  end;
  J := J xor Bit;

  if I < J then
  begin
   V := Data[I];
   Data[I] := Data[J];
   Data[J] := V;
  end;
 end;

 Len := 2;
 while Len <= N do
 begin
  HalfLen := Len shr 1;
  Ang := 2 * Pi / Len;
  if Invert then Ang := -Ang;

  SinCos(Ang, WlenIm, WlenRe);

  {WlenRe := Cos(Ang);
  WlenIm := Sin(Ang);}

  I := 0;
  while I < N do
  begin
   WRe := 1.0;
   WIm := 0.0;

   for K := 0 to HalfLen - 1 do
   begin
    URe := Data[I + K].Re;
    UIm := Data[I + K].Im;

    V := Data[I + K + HalfLen];

    // V * W
    VRe := V.Re * WRe - V.Im * WIm;
    VIm := V.Re * WIm + V.Im * WRe;

    Data[I + K].Re := URe + VRe;
    Data[I + K].Im := UIm + VIm;

    Data[I + K + HalfLen].Re := URe - VRe;
    Data[I + K + HalfLen].Im := UIm - VIm;

    // W *= Wlen
    TmpRe := WRe * WlenRe - WIm * WlenIm;
    WIm := WRe * WlenIm + WIm * WlenRe;
    WRe := TmpRe;
   end;

   Inc(I, Len);
  end;

  Len := Len shl 1;
 end;

 if Invert then
  for I := 0 to N - 1 do
   CScale(Data[I], 1.0 / N);
end;

procedure FFT2D(const Data: TArray<TComplex>; W, H: integer; Invert: boolean);
var
 X, Y: integer;
 Row:  TArray<TComplex>;
 Col:  TArray<TComplex>;
begin
 // rows
 SetLength(Row, W);
 for Y := 0 to H - 1 do
 begin
  for X := 0 to W - 1 do
   Row[X] := Data[Y * W + X];
  FFT1D(Row, Invert);
  for X := 0 to W - 1 do
   Data[Y * W + X] := Row[X];
 end;

 // cols
 SetLength(Col, H);
 for X := 0 to W - 1 do
 begin
  for Y := 0 to H - 1 do
   Col[Y] := Data[Y * W + X];
  FFT1D(Col, Invert);
  for Y := 0 to H - 1 do
   Data[Y * W + X] := Col[Y];
 end;
end;

procedure BuildHann1D(N: integer; out Wnd: TArray<double>);
var
 I: integer;
begin
 SetLength(Wnd, N);
 if N <= 1 then
 begin
  for I := 0 to N - 1 do Wnd[I] := 1.0;
  Exit;
 end;

 for I := 0 to N - 1 do
  Wnd[I] := 0.5 * (1.0 - Cos(2 * Pi * I / (N - 1)));
end;

function RawToHU(Raw: word; Slope, Intercept: double): double; inline;
begin
 Result := Raw * Slope + Intercept;
end;

function HUToRawWord(HU: double; Slope, Intercept: double): word;
var
 Raw: double;
 V: integer;
begin
 if Abs(Slope) < 1e-30 then
  Raw := HU
 else
  Raw := (HU - Intercept) / Slope;

 V := Round(Raw);
 if V < 0 then V := 0;
 if V > 65535 then V := 65535;
 Result := word(V);
end;

procedure Downsample2xWordAvg(const Src: TArray<word>; W, H: integer; out Dst: TArray<word>; out W2, H2: integer);
var
 X2, Y2: integer;
 X, Y: integer;
 Sum, Cnt: integer;
begin
 if Length(Src) <> W * H then
  raise Exception.Create('Downsample2xWordAvg: size mismatch');

 W2 := (W + 1) div 2;
 H2 := (H + 1) div 2;
 SetLength(Dst, W2 * H2);

 for Y2 := 0 to H2 - 1 do
  for X2 := 0 to W2 - 1 do
  begin
   Sum := 0;
   Cnt := 0;
   for Y := 2 * Y2 to 2 * Y2 + 1 do
    for X := 2 * X2 to 2 * X2 + 1 do
     if (X >= 0) and (X < W) and (Y >= 0) and (Y < H) then
     begin
      Inc(Sum, Src[Y * W + X]);
      Inc(Cnt);
     end;

   if Cnt = 0 then
    Dst[Y2 * W2 + X2] := 0
   else
    Dst[Y2 * W2 + X2] := word(Sum div Cnt);
  end;
end;

procedure ScaleROI2x(W2, H2: integer; RoiX, RoiY, RoiW, RoiH: integer; out RoiX2, RoiY2, RoiW2, RoiH2: integer);
begin
 if (RoiW > 0) and (RoiH > 0) then
 begin
  RoiX2 := EnsureRange(RoiX div 2, 0, W2 - 1);
  RoiY2 := EnsureRange(RoiY div 2, 0, H2 - 1);
  RoiW2 := EnsureRange(Max(1, RoiW div 2), 1, W2 - RoiX2);
  RoiH2 := EnsureRange(Max(1, RoiH div 2), 1, H2 - RoiY2);
 end
 else
 begin
  RoiX2 := 0;
  RoiY2 := 0;
  RoiW2 := 0;
  RoiH2 := 0;
 end;
end;

function OtsuThresholdHU(const Img: TArray<word>; W, H: integer; Slope, Intercept: double;
 RX, RY, RW, RH: integer): double;
const
 BINS = 4096;
var
 Hist: array[0..BINS - 1] of integer;
 X, Y: integer;
 Total: int64;
 SumAll, SumB: double;
 Wb, Wf: int64;
 Mb, Mf, VarBetween, MaxVar: double;
 ThrBin, B: integer;
 HU: double;

 function ClampBin(Value: integer): integer; inline;
 begin
  if Value < 0 then Exit(0);
  if Value > BINS - 1 then Exit(BINS - 1);
  Result := Value;
 end;

begin
 FillChar(Hist, SizeOf(Hist), 0);
 Total := 0;

 // Маппинг HU -> бин: диапазон [-1500..2500] HU ~ 4000 значений
 for Y := RY to RY + RH - 1 do
  for X := RX to RX + RW - 1 do
  begin
   HU := RawToHU(Img[Y * W + X], Slope, Intercept);
   B  := ClampBin(Round(HU + 1500));
   Inc(Hist[B]);
   Inc(Total);
  end;

 if Total = 0 then Exit(0);

 SumAll := 0.0;
 for B := 0 to BINS - 1 do
  SumAll := SumAll + B * Hist[B];

 SumB := 0.0;
 Wb := 0;
 MaxVar := -1.0;
 ThrBin := 0;

 for B := 0 to BINS - 1 do
 begin
  Wb := Wb + Hist[B];
  if Wb = 0 then Continue;

  Wf := Total - Wb;
  if Wf = 0 then Break;

  SumB := SumB + B * Hist[B];

  Mb := SumB / Wb;
  Mf := (SumAll - SumB) / Wf;

  VarBetween := Wb * Wf * Sqr(Mb - Mf);
  if VarBetween > MaxVar then
  begin
   MaxVar := VarBetween;
   ThrBin := B;
  end;
 end;

 Result := ThrBin - 1500; // HU
end;

procedure PreprocessToFeature(const Img: PDAWord; W, H: integer; Slope, Intercept: double;
 UseEdges: boolean; RoiX, RoiY, RoiW, RoiH: integer; out Feat: TArray<double>);
var
 HannX, HannY: TArray<double>;
 RX, RY, RW, RH: integer;
 ThrHU: double;
 X, Y: integer;
 I: integer;

 function HU(Xi, Yi: integer): double;
 begin
  Result := RawToHU(Img^[Yi * W + Xi], Slope, Intercept);
 end;

 function IsHead(Xi, Yi: integer; ThrHU: double): boolean; inline;
 begin
  Result := HU(Xi, Yi) > ThrHU;
 end;

var
 Gx, Gy, V: double;
begin
 SetLength(Feat, W * H);
 for I := 0 to High(Feat) do Feat[I] := 0.0;

 // ROI: если не задана, берём весь кадр
 if (RoiW > 0) and (RoiH > 0) then
 begin
  RX := EnsureRange(RoiX, 0, W - 1);
  RY := EnsureRange(RoiY, 0, H - 1);
  RW := EnsureRange(RoiW, 1, W - RX);
  RH := EnsureRange(RoiH, 1, H - RY);
 end
 else
 begin
  RX := 0;
  RY := 0;
  RW := W;
  RH := H;
 end;

 ThrHU := OtsuThresholdHU(Img^, W, H, Slope, Intercept, RX, RY, RW, RH);

 BuildHann1D(W, HannX);
 BuildHann1D(H, HannY);

 // Заполняем только внутри ROI, остальное нули
 for Y := Max(1, RY) to Min(H - 2, RY + RH - 2) do
  for X := Max(1, RX) to Min(W - 2, RX + RW - 2) do
  begin
   if not IsHead(X, Y, ThrHU) then
    Continue;

   if UseEdges then
   begin
    // Sobel по HU
    Gx :=
     (HU(X + 1, Y - 1) + 2 * HU(X + 1, Y) + HU(X + 1, Y + 1)) - (HU(X - 1, Y - 1) + 2 *
     HU(X - 1, Y) + HU(X - 1, Y + 1));

    Gy :=
     (HU(X - 1, Y + 1) + 2 * HU(X, Y + 1) + HU(X + 1, Y + 1)) - (HU(X - 1, Y - 1) + 2 *
     HU(X, Y - 1) + HU(X + 1, Y - 1));

    V := (Abs(Gx) + Abs(Gy)) * HannX[X] * HannY[Y];
   end
   else
    V := HU(X, Y) * HannX[X] * HannY[Y];

   Feat[Y * W + X] := V;
  end;
end;

procedure RotateFeatureBilinear(const Src: TArray<double>; W, H: integer; AngleRad: double; out Dst: TArray<double>);
var
 X, Y: integer;
 Cx, Cy: double;
 CosA, SinA: double;
 Xc, Yc, Xs, Ys: double;
 X0, Y0: integer;
 Fx, Fy: double;
 V00, V10, V01, V11: double;

 function GetF(Xi, Yi: integer): double;
 begin
  Result := Src[Yi * W + Xi];
 end;

begin
 SetLength(Dst, W * H);
 for X := 0 to High(Dst) do Dst[X] := 0.0;

 Cx := (W - 1) * 0.5;
 Cy := (H - 1) * 0.5;

 SinCos(AngleRad, SinA, CosA);

 {CosA := Cos(AngleRad);
 SinA := Sin(AngleRad);}

 for Y := 0 to H - 1 do
  for X := 0 to W - 1 do
  begin
   Xc := X - Cx;
   Yc := Y - Cy;

   Xs := CosA * Xc + SinA * Yc + Cx;
   Ys := -SinA * Xc + CosA * Yc + Cy;

   X0 := Floor(Xs);
   Y0 := Floor(Ys);

   if (X0 < 0) or (Y0 < 0) or (X0 >= W - 1) or (Y0 >= H - 1) then
    Continue;

   Fx := Xs - X0;
   Fy := Ys - Y0;

   V00 := GetF(X0, Y0);
   V10 := GetF(X0 + 1, Y0);
   V01 := GetF(X0, Y0 + 1);
   V11 := GetF(X0 + 1, Y0 + 1);

   Dst[Y * W + X] :=
    (1 - Fx) * (1 - Fy) * V00 + (Fx) * (1 - Fy) * V10 + (1 - Fx) * (Fy) * V01 + (Fx) * (Fy) * V11;
  end;
end;

procedure PadFeatureToComplex(const Feat: TArray<double>; W, H, PW, PH: integer; out C: TArray<TComplex>);
var
 Idx, X, Y: integer;
begin
 SetLength(C, PW * PH);
 for Idx := 0 to High(C) do
 begin
  C[Idx].Re := 0.0;
  C[Idx].Im := 0.0;
 end;

 for Y := 0 to H - 1 do
  for X := 0 to W - 1 do
  begin
   Idx := Y * PW + X;
   C[Idx].Re := Feat[Y * W + X];
  end;
end;

function ParabolicSubpixel(Fm1, F0, Fp1: double): double; inline;
var
 Den: double;
begin
 Den := (Fm1 - 2 * F0 + Fp1);
 if Abs(Den) < 1e-18 then Exit(0.0);
 Result := 0.5 * (Fm1 - Fp1) / Den;
end;

procedure PhaseCorrShiftPeak(const FRefFFT, FMovFFT: TArray<TComplex>; PW, PH: integer;
 MaxShiftPx: integer; UseSubpixel: boolean; out Dx, Dy, Peak: double);
var
 CPS, Corr: TArray<TComplex>;
 I, X, Y: integer;
 Mag, V0: double;
 MaxVal:  double;
 PeakX, PeakY: integer;
 SubX, SubY: double;
 Vxm1, Vxp1, Vym1, Vyp1: double;

 function Idx2(Xi, Yi, PW: integer): integer; inline;
 begin
  Result := Yi * PW + Xi;
 end;

 function WrapX(Xi, PW: integer): integer; inline;
 begin
  if Xi < 0 then Result := Xi + PW
  else
  if Xi >= PW then Result := Xi - PW
  else
   Result := Xi;
 end;

 function WrapY(Yi, PH: integer): integer; inline;
 begin
  if Yi < 0 then Result := Yi + PH
  else
  if Yi >= PH then Result := Yi - PH
  else
   Result := Yi;
 end;

 function WithinMaxShift(Xi, Yi, MaxShiftPx, PW, PH: integer): boolean; inline;
 var
  SX, SY: integer;
 begin
  if MaxShiftPx <= 0 then Exit(True);
  SX := Xi;
  SY := Yi;
  if SX > PW div 2 then SX := SX - PW;
  if SY > PH div 2 then SY := SY - PH;
  Result := (Abs(SX) <= MaxShiftPx) and (Abs(SY) <= MaxShiftPx);
 end;

begin
 SetLength(CPS, PW * PH);
 for I := 0 to PW * PH - 1 do
 begin
  CPS[I] := CMul(FRefFFT[I], CConj(FMovFFT[I]));
  Mag := Hypot(CPS[I].Re, CPS[I].Im);
  if Mag < EPS then
  begin
   CPS[I].Re := 0.0;
   CPS[I].Im := 0.0;
  end
  else
  begin
   CPS[I].Re := CPS[I].Re / Mag;
   CPS[I].Im := CPS[I].Im / Mag;
  end;
 end;

 Corr := CPS;
 FFT2D(Corr, PW, PH, True);

 // find peak (optionally restricted)
 MaxVal := -1e300;
 PeakX  := 0;
 PeakY  := 0;

 for Y := 0 to PH - 1 do
  for X := 0 to PW - 1 do
  begin
   if not WithinMaxShift(X, Y, MaxShiftPx, PW, PH) then Continue;
   V0 := Corr[Idx2(X, Y, PW)].Re;
   if V0 > MaxVal then
   begin
    MaxVal := V0;
    PeakX  := X;
    PeakY  := Y;
   end;
  end;

 // subpixel
 SubX := 0.0;
 SubY := 0.0;
 if UseSubpixel then
 begin
  V0 := Corr[Idx2(PeakX, PeakY, PW)].Re;

  Vxm1 := Corr[Idx2(WrapX(PeakX - 1, PW), PeakY, PW)].Re;
  Vxp1 := Corr[Idx2(WrapX(PeakX + 1, PW), PeakY, PW)].Re;
  SubX := ParabolicSubpixel(Vxm1, V0, Vxp1);

  Vym1 := Corr[Idx2(PeakX, WrapY(PeakY - 1, PH), PW)].Re;
  Vyp1 := Corr[Idx2(PeakX, WrapY(PeakY + 1, PH), PW)].Re;
  SubY := ParabolicSubpixel(Vym1, V0, Vyp1);
 end;

 // unwrap to signed shift
 if PeakX > PW div 2 then PeakX := PeakX - PW;
 if PeakY > PH div 2 then PeakY := PeakY - PH;

 Dx := PeakX + SubX;
 Dy := PeakY + SubY;
 Peak := MaxVal;
end;

procedure ApplyRigidBilinearWord(const Src: PDAWord; W, H: integer; AngleRad: double;
 Dx, Dy: double; FillValue: word; out Dst: TArray<word>);
var
 X, Y: integer;
 Cx, Cy: double;
 CosA, SinA: double;
 X1, Y1, Xc, Yc, Xs, Ys: double;
 X0, Y0: integer;
 Fx, Fy: double;
 I00, I10, I01, I11: integer;
 V: double;

 function GetPix(Xi, Yi: integer): integer;
 begin
  Result := Src^[Yi * W + Xi];
 end;

 function ClampWord(Value: integer): word; inline;
 begin
  if Value < 0 then Exit(0);
  if Value > 65535 then Exit(65535);
  Result := word(Value);
 end;

begin
 SetLength(Dst, W * H);

 Cx := (W - 1) * 0.5;
 Cy := (H - 1) * 0.5;

 SinCos(AngleRad, SinA, CosA);

 {CosA := Cos(AngleRad);
 SinA := Sin(AngleRad);}

 for Y := 0 to H - 1 do
  for X := 0 to W - 1 do
  begin
   X1 := X - Dx;
   Y1 := Y - Dy;

   Xc := X1 - Cx;
   Yc := Y1 - Cy;

   Xs := CosA * Xc + SinA * Yc + Cx;
   Ys := -SinA * Xc + CosA * Yc + Cy;

   X0 := Floor(Xs);
   Y0 := Floor(Ys);

   if (X0 < 0) or (Y0 < 0) or (X0 >= W - 1) or (Y0 >= H - 1) then
   begin
    Dst[Y * W + X] := FillValue;
    Continue;
   end;

   Fx := Xs - X0;
   Fy := Ys - Y0;

   I00 := GetPix(X0, Y0);
   I10 := GetPix(X0 + 1, Y0);
   I01 := GetPix(X0, Y0 + 1);
   I11 := GetPix(X0 + 1, Y0 + 1);

   V :=
    (1 - Fx) * (1 - Fy) * I00 + (Fx) * (1 - Fy) * I10 + (1 - Fx) * (Fy) * I01 + (Fx) * (Fy) * I11;

   Dst[Y * W + X] := ClampWord(Round(V));
  end;
end;

function EstimateRigidAngleCentered(const RefFFT: TArray<TComplex>; PW, PH: integer;
 const MovFeat: TArray<double>; W, H: integer; MaxShiftPx: integer; UseSubpixel: boolean;
 CenterAngleRad: double; SpanDeg: double; StepDeg: double): TRigid;
var
 Best: TRigid;
 CenterDeg, AngleDeg: double;
 AngleRad: double;
 RotFeat: TArray<double>;
 MovC, MovFFT: TArray<TComplex>;
 Dx, Dy, Peak: double;

 procedure EvalAngle(aRad: double);
 begin
  RotateFeatureBilinear(MovFeat, W, H, aRad, RotFeat);
  PadFeatureToComplex(RotFeat, W, H, PW, PH, MovC);
  MovFFT := MovC;
  FFT2D(MovFFT, PW, PH, False);

  PhaseCorrShiftPeak(RefFFT, MovFFT, PW, PH, MaxShiftPx, UseSubpixel, Dx, Dy, Peak);

  if Peak > Best.Peak then
  begin
   Best.AngleRad := aRad;
   Best.Dx := Dx;
   Best.Dy := Dy;
   Best.Peak := Peak;
  end;
 end;

begin
 Best.AngleRad := CenterAngleRad;
 Best.Dx := 0;
 Best.Dy := 0;
 Best.Peak := -1e300;

 if StepDeg <= 0 then StepDeg := 1.0;
 if SpanDeg < 0 then SpanDeg := 0;

 CenterDeg := CenterAngleRad * 180.0 / Pi;

 AngleDeg := CenterDeg - SpanDeg;
 while AngleDeg <= CenterDeg + SpanDeg + 1e-9 do
 begin
  AngleRad := AngleDeg * Pi / 180.0;
  EvalAngle(AngleRad);
  AngleDeg := AngleDeg + StepDeg;
 end;

 // tiny local refinement
 EvalAngle(Best.AngleRad - (StepDeg * 0.25) * Pi / 180.0);
 EvalAngle(Best.AngleRad + (StepDeg * 0.25) * Pi / 180.0);
 EvalAngle(Best.AngleRad - (StepDeg * 0.10) * Pi / 180.0);
 EvalAngle(Best.AngleRad + (StepDeg * 0.10) * Pi / 180.0);

 Result := Best;
end;

function EstimateRigidAngleGlobal(const RefFFT: TArray<TComplex>; PW, PH: integer;
 const MovFeat: TArray<double>; W, H: integer; MaxShiftPx: integer; UseSubpixel: boolean;
 MaxAngleDeg: double; StepDeg: double): TRigid;
begin
 Result := EstimateRigidAngleCentered(RefFFT, PW, PH, MovFeat, W, H, MaxShiftPx, UseSubpixel,
  0.0, Abs(MaxAngleDeg), StepDeg);
end;

procedure RegisterCoreRigidPyramidInPlace(const Frames: TArray<PDAWord>; W, H: integer;
 Slope, Intercept: double; FillValue: word; MaxShiftPx: integer; UseSubpixel, UseEdges: boolean;
 RoiX, RoiY, RoiW, RoiH: integer; MaxAngleDeg, AngleStepDeg: double);
var
 // full
 RefFeat0: TArray<double>;
 PW0, PH0: integer;
 RefC0, RefFFT0: TArray<TComplex>;
 RefImg1:  TArray<word>;
 // half
 W1, H1, PW1, PH1: integer;
 RoiX1, RoiY1, RoiW1, RoiH1: integer;
 RefFeat1: TArray<double>;
 RefC1, RefFFT1: TArray<TComplex>;
 MaxShift1: integer;
 FineSpanDeg, FineStepDeg: double;
begin
 if Length(Frames) = 0 then Exit;
 if Length(Frames[0]^) <> W * H then
  raise Exception.Create('RegisterCoreRigidPyramidInPlace: reference frame size mismatch');

 // level 0 ref FFT
 PW0 := NextPow2(W);
 PH0 := NextPow2(H);

 PreprocessToFeature(Frames[0], W, H, Slope, Intercept, UseEdges, RoiX, RoiY, RoiW, RoiH, RefFeat0);
 PadFeatureToComplex(RefFeat0, W, H, PW0, PH0, RefC0);
 RefFFT0 := RefC0;
 FFT2D(RefFFT0, PW0, PH0, False);

 // level 1 ref FFT
 Downsample2xWordAvg(Frames[0]^, W, H, RefImg1, W1, H1);
 PW1 := NextPow2(W1);
 PH1 := NextPow2(H1);

 ScaleROI2x(W1, H1, RoiX, RoiY, RoiW, RoiH, RoiX1, RoiY1, RoiW1, RoiH1);
 PreprocessToFeature(@RefImg1, W1, H1, Slope, Intercept, UseEdges, RoiX1, RoiY1, RoiW1, RoiH1, RefFeat1);
 PadFeatureToComplex(RefFeat1, W1, H1, PW1, PH1, RefC1);
 RefFFT1 := RefC1;
 FFT2D(RefFFT1, PW1, PH1, False);

 if MaxShiftPx > 0 then
  MaxShift1 := Max(1, MaxShiftPx div 2)
 else
  MaxShift1 := 0;

 FineSpanDeg := Min(Abs(MaxAngleDeg), Max(2.0, AngleStepDeg * 2.0));
 FineStepDeg := Max(0.2, AngleStepDeg * 0.25);

 TParallel.for(1, High(Frames),
  procedure(i: integer; LoopState: TParallel.TLoopState)
 var
  MovImg1: TArray<word>;
  MovFeat1, MovFeat0: TArray<double>;
  Rigid1, Rigid0: TRigid;
  Aligned: TArray<word>;
  W1, H1: integer;
 begin
  ParallelProgress;
  if Length(Frames[I]^) <> W * H then
   raise Exception.Create('RegisterCoreRigidPyramidInPlace: frame size mismatch');

  // coarse on half-res
  Downsample2xWordAvg(Frames[I]^, W, H, MovImg1, W1, H1);
  PreprocessToFeature(@MovImg1, W1, H1, Slope, Intercept, UseEdges, RoiX1, RoiY1, RoiW1, RoiH1, MovFeat1);

  Rigid1 := EstimateRigidAngleGlobal(RefFFT1, PW1, PH1, MovFeat1, W1, H1, MaxShift1, UseSubpixel,
  MaxAngleDeg, AngleStepDeg);

  // refine on full-res around coarse angle
  PreprocessToFeature(Frames[I], W, H, Slope, Intercept, UseEdges, RoiX, RoiY, RoiW, RoiH, MovFeat0);

  Rigid0 := EstimateRigidAngleCentered(RefFFT0, PW0, PH0, MovFeat0, W, H, MaxShiftPx,
  UseSubpixel, Rigid1.AngleRad, FineSpanDeg, FineStepDeg);

  ApplyRigidBilinearWord(
  Frames[I], W, H,
  Rigid0.AngleRad, Rigid0.Dx, Rigid0.Dy,
  FillValue, Aligned
  );

  Move(Aligned[0], Frames[I]^[0], Length(Aligned) * SizeOf(word));

  if ProgressForm.ProgressCancelled then
   LoopState.Break;
 end);
end;

Джейсоны, увы. Достаточно лютая математика, среди прочего. Область: медицина.

Ничего нового: «Вычислительная машина ценна ровно настолько, насколько ценен человек, который ее использует» (ц) Винер

1
23 ...

Информация

В рейтинге
3 643-й
Зарегистрирован
Активность