48 ProcessRubberSheet::ProcessRubberSheet(
int startSize,
int endSize) {
49 p_bandChangeFunct = NULL;
54 p_startQuadSize = startSize;
55 p_endQuadSize = endSize;
58 m_patchStartSample = 1;
62 m_patchSampleIncrement = 4;
63 m_patchLineIncrement = 4;
113 void ProcessRubberSheet::setPatchParameters(
int startSample,
int startLine,
114 int samples,
int lines,
115 int sampleIncrement,
int lineIncrement) {
116 m_patchStartSample = startSample;
117 m_patchStartLine = startLine;
118 m_patchSamples = samples;
119 m_patchLines = lines;
120 m_patchSampleIncrement = sampleIncrement;
121 m_patchLineIncrement = lineIncrement;
144 if(InputCubes.size() != 1) {
145 string m =
"You must specify exactly one input cube";
148 else if(OutputCubes.size() != 1) {
149 string m =
"You must specify exactly one output cube";
154 p_lineMap.resize(p_startQuadSize);
155 p_sampMap.resize(p_startQuadSize);
157 for(
unsigned int pos = 0; pos < p_lineMap.size(); pos++) {
158 p_lineMap[pos].resize(p_startQuadSize);
159 p_sampMap[pos].resize(p_startQuadSize);
163 TileManager otile(*OutputCubes[0], p_startQuadSize, p_startQuadSize);
167 InputCubes[0]->pixelType() ,
171 p_progress->SetMaximumSteps(otile.
Tiles());
172 p_progress->CheckStatus();
174 if(p_bandChangeFunct == NULL) {
180 int tilesPerBand = otile.
Tiles() / OutputCubes[0]->bandCount();
182 for(
int tile = 1; tile <= tilesPerBand; tile++) {
183 bool useLastTileMap =
false;
184 for(
int band = 1; band <= OutputCubes[0]->bandCount(); band++) {
187 if(p_startQuadSize <= 2) {
188 SlowGeom(otile, iportal, trans, interp);
191 QuadTree(otile, iportal, trans, interp, useLastTileMap);
194 useLastTileMap =
true;
196 OutputCubes[0]->write(otile);
197 p_progress->CheckStatus();
202 int lastOutputBand = -1;
204 for(otile.
begin(); !otile.
end(); otile++) {
206 if(lastOutputBand != otile.
Band()) {
207 lastOutputBand = otile.
Band();
209 p_bandChangeFunct(lastOutputBand);
212 if(p_startQuadSize <= 2) {
213 SlowGeom(otile, iportal, trans, interp);
216 QuadTree(otile, iportal, trans, interp,
false);
219 OutputCubes[0]->write(otile);
220 p_progress->CheckStatus();
237 void ProcessRubberSheet::BandChange(
void (*funct)(
const int band)) {
238 p_bandChangeFunct = funct;
243 double outputSamp, outputLine;
244 double inputSamp, inputLine;
245 int outputBand = otile.
Band();
247 for(
int i = 0; i < otile.
size(); i++) {
248 outputSamp = otile.
Sample(i);
249 outputLine = otile.
Line(i);
252 if(trans.
Xform(inputSamp, inputLine, outputSamp, outputLine)) {
253 if((inputSamp < 0.5) || (inputLine < 0.5) ||
254 (inputLine > InputCubes[0]->lineCount() + 0.5) ||
255 (inputSamp > InputCubes[0]->sampleCount() + 0.5)) {
260 iportal.
SetPosition(inputSamp, inputLine, outputBand);
261 InputCubes[0]->read(iportal);
262 otile[i] = interp.
Interpolate(inputSamp, inputLine,
272 void ProcessRubberSheet::QuadTree(TileManager &otile, Portal &iportal,
273 Transform &trans, Interpolator &interp,
274 bool useLastTileMap) {
276 vector<Quad *> quadTree;
278 if(!useLastTileMap) {
280 Quad *quad =
new Quad;
281 quad->sline = otile.Line();
282 quad->ssamp = otile.Sample();
284 quad->eline = otile.Line(otile.size() - 1);
285 quad->esamp = otile.Sample(otile.size() - 1);
286 quad->slineTile = otile.Line();
287 quad->ssampTile = otile.Sample();
289 quadTree.push_back(quad);
293 while(quadTree.size() > 0) {
294 ProcessQuad(quadTree, trans, p_lineMap, p_sampMap);
299 int outputBand = otile.Band();
300 for(
int i = 0, line = 0; line < p_startQuadSize; line++) {
301 for(
int samp = 0; samp < p_startQuadSize; samp++, i++) {
302 double inputLine = p_lineMap[line][samp];
303 double inputSamp = p_sampMap[line][samp];
304 if(inputLine != NULL8) {
305 iportal.SetPosition(inputSamp, inputLine, outputBand);
306 InputCubes[0]->read(iportal);
307 otile[i] = interp.Interpolate(inputSamp, inputLine,
308 iportal.DoubleBuffer());
330 bool ProcessRubberSheet::TestLine(
Transform &trans,
int ssamp,
int esamp,
331 int sline,
int eline,
int increment) {
332 for(
int line = sline; line <= eline; line += increment) {
333 for(
int sample = ssamp; sample <= esamp; sample += increment) {
337 if(trans.
Xform(sjunk, ljunk, sample, line)) {
347 void ProcessRubberSheet::ProcessQuad(std::vector<Quad *> &quadTree,
349 std::vector< std::vector<double> > &lineMap,
350 std::vector< std::vector<double> > &sampMap) {
351 Quad *quad = quadTree[0];
352 double oline[4], osamp[4];
353 double iline[4], isamp[4];
357 oline[0] = quad->sline;
358 osamp[0] = quad->ssamp;
359 if(!trans.
Xform(isamp[0], iline[0], osamp[0], oline[0])) {
364 oline[1] = quad->sline;
365 osamp[1] = quad->esamp;
366 if(!trans.
Xform(isamp[1], iline[1], osamp[1], oline[1])) {
371 oline[2] = quad->eline;
372 osamp[2] = quad->ssamp;
373 if(!trans.
Xform(isamp[2], iline[2], osamp[2], oline[2])) {
378 oline[3] = quad->eline;
379 osamp[3] = quad->esamp;
380 if(!trans.
Xform(isamp[3], iline[3], osamp[3], oline[3])) {
388 if((quad->eline - quad->sline) < p_endQuadSize) {
389 SlowQuad(quadTree, trans, lineMap, sampMap);
392 if(p_forceSamp !=
Null && p_forceLine !=
Null) {
393 if(p_forceSamp >= quad->ssamp && p_forceSamp <= quad->esamp &&
394 p_forceLine >= quad->sline && p_forceLine <= quad->eline) {
400 int centerSample = (quad->ssamp + quad->esamp) / 2;
401 int centerLine = (quad->sline + quad->eline) / 2;
418 if(TestLine(trans, quad->ssamp + 1, quad->esamp - 1,
419 quad->sline, quad->sline, 4) ||
421 TestLine(trans, quad->ssamp + 1, quad->esamp - 1,
422 quad->eline, quad->eline, 4) ||
424 TestLine(trans, quad->ssamp, quad->ssamp,
425 quad->sline + 1, quad->eline - 1, 4) ||
427 TestLine(trans, quad->esamp, quad->esamp,
428 quad->sline + 1, quad->eline - 1, 4) ||
430 TestLine(trans, centerSample, centerSample,
431 quad->sline + 1, quad->eline - 1, 4) ||
433 TestLine(trans, quad->ssamp + 1, quad->esamp - 1,
434 centerLine, centerLine, 4)) {
442 for(
int i = quad->sline; i <= quad->eline; i++) {
443 for(
int j = quad->ssamp; j <= quad->esamp; j++) {
444 lineMap[i-quad->slineTile][j-quad->ssampTile] = NULL8;
448 quadTree.erase(quadTree.begin());
472 if((quad->eline - quad->sline) < p_endQuadSize) {
473 SlowQuad(quadTree, trans, lineMap, sampMap);
486 for(
int i = 0; i < 4; i++) {
489 A[i][2] = oline[i] * osamp[i];
497 if((detA = Det4x4(A)) == 0.0) {
498 if((quad->eline - quad->sline) < p_endQuadSize) {
499 SlowQuad(quadTree, trans, lineMap, sampMap);
511 for(
int j = 0; j < 4; j++) {
512 memmove(B, A, 16 *
sizeof(
double));
514 for(
int i = 0; i < 4; i++) {
518 lineCoef[j] = Det4x4(B) / detA;
523 for(
int j = 0; j < 4; j++) {
524 memmove(B, A, 16 *
sizeof(
double));
526 for(
int i = 0; i < 4; i++) {
530 sampCoef[j] = Det4x4(B) / detA;
534 double quadMidLine = (quad->sline + quad->eline) / 2.0;
535 double quadMidSamp = (quad->ssamp + quad->esamp) / 2.0;
536 double midLine, midSamp;
538 if(!trans.
Xform(midSamp, midLine, quadMidSamp, quadMidLine)) {
539 if((quad->eline - quad->sline) < p_endQuadSize) {
540 SlowQuad(quadTree, trans, lineMap, sampMap);
548 double cmidLine = lineCoef[0] * quadMidLine +
549 lineCoef[1] * quadMidSamp +
550 lineCoef[2] * quadMidLine * quadMidSamp +
553 double cmidSamp = sampCoef[0] * quadMidLine +
554 sampCoef[1] * quadMidSamp +
555 sampCoef[2] * quadMidLine * quadMidSamp +
558 if((abs(cmidSamp - midSamp) > 0.5) || (abs(cmidLine - midLine) > 0.5)) {
559 if((quad->eline - quad->sline) < p_endQuadSize) {
560 SlowQuad(quadTree, trans, lineMap, sampMap);
570 double ulLine = lineCoef[0] * (double) quad->sline +
571 lineCoef[1] * (
double) quad->ssamp +
572 lineCoef[2] * (double) quad->sline * (
double) quad->ssamp +
575 double ulSamp = sampCoef[0] * (double) quad->sline +
576 sampCoef[1] * (
double) quad->ssamp +
577 sampCoef[2] * (double) quad->sline * (
double) quad->ssamp +
582 double lineChangeWrLine = lineCoef[0] + lineCoef[2] * (double) quad->ssamp;
583 double sampChangeWrLine = sampCoef[0] + sampCoef[2] * (
double) quad->ssamp;
585 for(
int ol = quad->sline; ol <= quad->eline; ol++) {
588 double lineChangeWrSamp = lineCoef[1] + lineCoef[2] * (double) ol;
589 double sampChangeWrSamp = sampCoef[1] + sampCoef[2] * (double) ol;
592 double cline = ulLine;
593 double csamp = ulSamp;
596 int startSamp = quad->ssamp - quad->ssampTile;
597 std::vector<double> &lineVect = lineMap[ol-quad->slineTile];
598 std::vector<double> &sampleVect = sampMap[ol-quad->slineTile];
601 for(
int os = quad->ssamp; os <= quad->esamp; os++) {
602 lineVect[startSamp] = cline;
603 sampleVect[startSamp] = csamp;
607 cline += lineChangeWrSamp;
608 csamp += sampChangeWrSamp;
612 ulLine += lineChangeWrLine;
613 ulSamp += sampChangeWrLine;
618 quadTree.erase(quadTree.begin());
622 void ProcessRubberSheet::SplitQuad(std::vector<Quad *> &quadTree) {
624 Quad *quad = quadTree[0];
625 int n = (quad->eline - quad->sline + 1) / 2;
630 q1->eline = quad->sline + n - 1;
631 q1->esamp = quad->ssamp + n - 1;
632 quadTree.push_back(q1);
637 q2->eline = quad->sline + n - 1;
638 q2->ssamp = quad->ssamp + n;
639 quadTree.push_back(q2);
644 q3->sline = quad->sline + n;
645 q3->esamp = quad->ssamp + n - 1;
646 quadTree.push_back(q3);
651 q4->sline = quad->sline + n;
652 q4->ssamp = quad->ssamp + n;
653 quadTree.push_back(q4);
657 quadTree.erase(quadTree.begin());
662 void ProcessRubberSheet::SlowQuad(std::vector<Quad *> &quadTree,
664 std::vector< std::vector<double> > &lineMap,
665 std::vector< std::vector<double> > &sampMap)
668 Quad *quad = quadTree[0];
672 for(
int oline = quad->sline; oline <= quad->eline; oline++) {
673 int lineIndex = oline - quad->slineTile;
674 for(
int osamp = quad->ssamp; osamp <= quad->esamp; osamp++) {
675 int sampIndex = osamp - quad->ssampTile;
676 lineMap[lineIndex][sampIndex] = NULL8;
677 if(trans.Xform(isamp, iline, (
double) osamp, (
double) oline)) {
680 (iline <= InputCubes[0]->lineCount() + 0.5) ||
681 (isamp <= InputCubes[0]->sampleCount() + 0.5)) {
682 lineMap[lineIndex][sampIndex] = iline;
683 sampMap[lineIndex][sampIndex] = isamp;
691 quadTree.erase(quadTree.begin());
695 double ProcessRubberSheet::Det4x4(
double m[4][4]) {
698 cofact [0][0] = m[1][1];
699 cofact [0][1] = m[1][2];
700 cofact [0][2] = m[1][3];
701 cofact [1][0] = m[2][1];
702 cofact [1][1] = m[2][2];
703 cofact [1][2] = m[2][3];
704 cofact [2][0] = m[3][1];
705 cofact [2][1] = m[3][2];
706 cofact [2][2] = m[3][3];
707 double det = m[0][0] * Det3x3(cofact);
709 cofact [0][0] = m[1][0];
710 cofact [0][1] = m[1][2];
711 cofact [0][2] = m[1][3];
712 cofact [1][0] = m[2][0];
713 cofact [1][1] = m[2][2];
714 cofact [1][2] = m[2][3];
715 cofact [2][0] = m[3][0];
716 cofact [2][1] = m[3][2];
717 cofact [2][2] = m[3][3];
718 det -= m[0][1] * Det3x3(cofact);
720 cofact [0][0] = m[1][0];
721 cofact [0][1] = m[1][1];
722 cofact [0][2] = m[1][3];
723 cofact [1][0] = m[2][0];
724 cofact [1][1] = m[2][1];
725 cofact [1][2] = m[2][3];
726 cofact [2][0] = m[3][0];
727 cofact [2][1] = m[3][1];
728 cofact [2][2] = m[3][3];
729 det += m[0][2] * Det3x3(cofact);
731 cofact [0][0] = m[1][0];
732 cofact [0][1] = m[1][1];
733 cofact [0][2] = m[1][2];
734 cofact [1][0] = m[2][0];
735 cofact [1][1] = m[2][1];
736 cofact [1][2] = m[2][2];
737 cofact [2][0] = m[3][0];
738 cofact [2][1] = m[3][1];
739 cofact [2][2] = m[3][2];
740 det -= m[0][3] * Det3x3(cofact);
746 double ProcessRubberSheet::Det3x3(
double m[3][3]) {
747 return m[0][0] * m[1][1] * m[2][2] -
748 m[0][0] * m[1][2] * m[2][1] -
749 m[0][1] * m[1][0] * m[2][2] +
750 m[0][1] * m[1][2] * m[2][0] +
751 m[0][2] * m[1][0] * m[2][1] -
752 m[0][2] * m[1][1] * m[2][0];
778 void ProcessRubberSheet::processPatchTransform(
Transform &trans,
781 if(InputCubes.size() != 1) {
782 string m =
"You must specify exactly one input cube";
785 else if(OutputCubes.size() != 1) {
786 string m =
"You must specify exactly one output cube";
792 InputCubes[0]->pixelType() ,
796 Brick obrick(*OutputCubes[0],1,1,1);
799 int patchesPerBand = 0;
800 for (
int line = m_patchStartLine; line <= InputCubes[0]->lineCount();
801 line += m_patchLineIncrement) {
802 for (
int samp = m_patchStartSample;
803 samp <= InputCubes[0]->sampleCount();
804 samp += m_patchSampleIncrement) {
809 int patchCount = InputCubes[0]->bandCount() * patchesPerBand;
810 p_progress->SetMaximumSteps(patchCount);
811 p_progress->CheckStatus();
833 for (
int band=1; band <= InputCubes[0]->bandCount(); band++) {
834 if(p_bandChangeFunct != NULL) p_bandChangeFunct(band);
835 iportal.SetPosition(1,1,band);
837 for (
int line = m_patchStartLine; line <= InputCubes[0]->lineCount();
838 line += m_patchLineIncrement) {
839 for (
int samp = m_patchStartSample;
840 samp <= InputCubes[0]->sampleCount();
841 samp += m_patchSampleIncrement, p_progress->CheckStatus()) {
842 transformPatch((
double)samp, (
double)(samp + m_patchSamples - 1),
843 (
double)line, (
double)(line + m_patchLines - 1),
844 obrick, iportal, trans, interp);
852 void ProcessRubberSheet::transformPatch (
double ssamp,
double esamp,
853 double sline,
double eline,
859 if (esamp > InputCubes[0]->sampleCount()) {
860 esamp = InputCubes[0]->sampleCount();
861 if (ssamp == esamp) ssamp = esamp - 1;
863 if (eline > InputCubes[0]->lineCount()) {
864 eline = InputCubes[0]->lineCount();
865 if (sline == eline) sline = eline - 1;
876 if (trans.
Xform(tsamp,tline,ssamp,sline)) {
877 isamps.push_back(ssamp);
878 ilines.push_back(sline);
879 osamps.push_back(tsamp);
880 olines.push_back(tline);
884 if (trans.
Xform(tsamp,tline,esamp,sline)) {
885 isamps.push_back(esamp);
886 ilines.push_back(sline);
887 osamps.push_back(tsamp);
888 olines.push_back(tline);
892 if (trans.
Xform(tsamp,tline,ssamp,eline)) {
893 isamps.push_back(ssamp);
894 ilines.push_back(eline);
895 osamps.push_back(tsamp);
896 olines.push_back(tline);
900 if (trans.
Xform(tsamp,tline,esamp,eline)) {
901 isamps.push_back(esamp);
902 ilines.push_back(eline);
903 osamps.push_back(tsamp);
904 olines.push_back(tline);
910 if (isamps.size() == 0)
return;
913 if (isamps.size() < 4) {
914 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
919 int osampMin = (int) osamps[0];
920 int osampMax = (int) (osamps[0] + 0.5);
921 int olineMin = (int) olines[0];
922 int olineMax = (int) (olines[0] + 0.5);
923 for (
int i = 1; i < osamps.size(); i++) {
924 if (osamps[i] < osampMin) osampMin = (int) osamps[i];
925 if (osamps[i] > osampMax) osampMax = (int) (osamps[i] + 0.5);
926 if (olines[i] < olineMin) olineMin = (int) olines[i];
927 if (olines[i] > olineMax) olineMax = (int) (olines[i] + 0.5);
932 if (osampMax < 1)
return;
933 if (olineMax < 1)
return;
934 if (osampMin > OutputCubes[0]->sampleCount())
return;
935 if (olineMin > OutputCubes[0]->lineCount())
return;
938 if (osampMin < 1) osampMin = 1;
939 if (olineMin < 1) olineMin = 1;
940 if (osampMax > OutputCubes[0]->sampleCount()) {
941 osampMax = OutputCubes[0]->sampleCount();
943 if (olineMax > OutputCubes[0]->lineCount()) {
944 olineMax = OutputCubes[0]->lineCount();
956 if (osampMax - osampMin + 1.0 > OutputCubes[0]->sampleCount() * 0.50) {
957 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
960 if (olineMax - olineMin + 1.0 > OutputCubes[0]->lineCount() * 0.50) {
961 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
966 BasisFunction isampFunc(
"Ax+By+C",3,3);
967 LeastSquares isampLSQ(isampFunc);
969 BasisFunction ilineFunc(
"Dx+Ey+F",3,3);
970 LeastSquares ilineLSQ(ilineFunc);
973 for (
int i=0; i < isamps.size(); i++) {
974 std::vector<double> vars;
975 vars.push_back(osamps[i]);
976 vars.push_back(olines[i]);
979 isampLSQ.AddKnown(vars,isamps[i]);
980 ilineLSQ.AddKnown(vars,ilines[i]);
983 isampLSQ.Solve(LeastSquares::QRD);
984 ilineLSQ.Solve(LeastSquares::QRD);
986 catch (IException &e) {
987 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
992 for (
int i=0; i<isamps.size(); i++) {
993 if (fabs(isampLSQ.Residual(i)) > 0.5) {
994 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
997 if (fabs(ilineLSQ.Residual(i)) > 0.5) {
998 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
1007 double csamp = (ssamp + esamp) / 2.0;
1008 double cline = (sline + eline) / 2.0;
1009 if (trans.
Xform(tsamp,tline,csamp,cline)) {
1010 std::vector<double> vars;
1011 vars.push_back(tsamp);
1012 vars.push_back(tline);
1013 vars.push_back(1.0);
1015 double isamp = isampFunc.Evaluate(vars);
1016 double iline = ilineFunc.Evaluate(vars);
1018 double err = (csamp - isamp) * (csamp - isamp) +
1019 (cline - iline) * (cline - iline);
1021 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
1026 splitPatch(ssamp, esamp, sline, eline, obrick, iportal, trans, interp);
1031 double A = isampFunc.Coefficient(0);
1032 double B = isampFunc.Coefficient(1);
1033 double C = isampFunc.Coefficient(2);
1035 double D = ilineFunc.Coefficient(0);
1036 double E = ilineFunc.Coefficient(1);
1037 double F = ilineFunc.Coefficient(2);
1041 for (
int oline = olineMin; oline <= olineMax; oline++) {
1042 double isamp = A * osampMin + B * oline + C;
1043 double iline = D * osampMin + E * oline + F;
1044 double isampChangeWRTosamp = A;
1045 double ilineChangeWRTosamp = D;
1046 for (
int osamp = osampMin; osamp <= osampMax; osamp++,
1047 isamp += isampChangeWRTosamp, iline += ilineChangeWRTosamp) {
1051 InputCubes[0]->read(iportal);
1058 OutputCubes[0]->write(obrick);
1066 void ProcessRubberSheet::splitPatch (
double ssamp,
double esamp,
1067 double sline,
double eline,
1068 Brick &obrick, Portal &iportal,
1069 Transform &trans, Interpolator &interp) {
1072 if ((esamp - ssamp < 0.1) && (eline - sline < 0.1))
return;
1075 double midSamp = (esamp + ssamp) / 2.0;
1076 double midLine = (eline + sline) / 2.0;
1078 transformPatch(ssamp, midSamp,
1080 obrick, iportal, trans, interp);
1081 transformPatch(midSamp, esamp,
1083 obrick, iportal, trans, interp);
1084 transformPatch(ssamp, midSamp,
1086 obrick, iportal, trans, interp);
1087 transformPatch(midSamp, esamp,
1089 obrick, iportal, trans, interp);