11 #include "IException.h"
12 #include "Interpolator.h"
14 #include "LineManager.h"
15 #include "PolygonTools.h"
17 #include "Projection.h"
18 #include "Statistics.h"
19 #include "TProjection.h"
26 #include <geos/geom/Point.h>
27 #include <tnt/tnt_array2d_utils.h>
44 Chip::Chip(
const Chip &other) {
63 m_clipPolygon = PolygonTools::MakeMultiPolygon(
82 Chip::Chip(
const int samples,
const int lines) {
89 if (m_clipPolygon != NULL)
delete m_clipPolygon;
100 void Chip::SetAllValues(
const double &d) {
101 for (
unsigned int i = 0; i < m_buf.size(); i++) {
102 fill(m_buf[i].begin(), m_buf[i].end(), d);
116 void Chip::Init(
const int samples,
const int lines) {
117 SetReadInterpolator(Interpolator::CubicConvolutionType);
118 SetSize(samples, lines);
120 m_clipPolygon = NULL;
134 void Chip::SetSize(
const int samples,
const int lines) {
135 if ( samples <= 0.0 || lines <= 0.0 ) {
137 msg +=
"Unable to set chip size to [";
139 msg +=
toString(lines) +
"]. Samples and lines must be greater than zero.";
140 throw IException(IException::User, msg, _FILEINFO_);
142 m_chipSamples = samples;
146 for (
int i = 0; i < lines; i++) {
147 m_buf[i].resize(samples);
150 m_tackSample = ((samples - 1) / 2) + 1;
151 m_tackLine = ((lines - 1) / 2) + 1;
162 bool Chip::IsInsideChip(
double sample,
double line) {
163 double minSamp = m_cubeTackSample - ((m_chipSamples - 1) / 2);
164 double maxSamp = m_cubeTackSample + ((m_chipSamples - 1) / 2);
166 double minLine = m_cubeTackLine - ((m_chipLines - 1) / 2);
167 double maxLine = m_cubeTackLine + ((m_chipLines - 1) / 2);
169 if ( sample < minSamp || sample > maxSamp )
return false;
170 if ( line < minLine || line > maxLine )
return false;
182 void Chip::TackCube(
const double cubeSample,
const double cubeLine) {
183 m_cubeTackSample = cubeSample;
184 m_cubeTackLine = cubeLine;
186 m_affine.Translate(m_cubeTackSample, m_cubeTackLine);
203 void Chip::Load(
Cube &cube,
const double rotation,
const double scale,
const int band) {
209 m_affine.Scale(scale);
210 m_affine.Rotate(rotation);
211 m_affine.Translate(m_cubeTackSample, m_cubeTackLine);
244 void Chip::Load(
Cube &cube,
const Affine &affine,
const bool &keepPoly,
const int band) {
247 SetTransform(affine);
248 SetChipPosition(TackSample(), TackLine());
252 delete m_clipPolygon;
300 void Chip::Load(
Cube &cube,
Chip &match,
Cube &matchChipCube,
const double scale,
const int band)
306 matchCam = matchChipCube.
camera();
313 QString msg =
"Can not geom chip. ";
314 msg +=
"Match chip cube [" + matchChipCube.
fileName();
315 msg +=
"] is not a camera or map projection";
317 IException fullError(IException::User, msg, _FILEINFO_);
336 QString msg =
"Can not geom chip. ";
337 msg +=
"Chip cube [" + cube.
fileName();
338 msg +=
"] is not a camera or map projection";
340 IException fullError(IException::User, msg, _FILEINFO_);
351 vector<double> x(4), y(4);
352 vector<double> xp(4), yp(4);
360 for (
int i = 0; i < (int) xp.size(); i++) {
364 int endSamp = Samples() - 1;
365 int endLine = Lines() - 1;
367 bool pointfound =
false;
368 while (!pointfound) {
371 if (startSamp < 1 || startSamp > Samples() - 1 ||
372 endSamp < 1 || endSamp > Samples() - 1 ||
373 startLine < 1 || startLine > Lines() - 1 ||
374 endLine < 1 || endLine > Lines() - 1) {
377 x.erase(x.begin() + i);
378 y.erase(y.begin() + i);
379 xp.erase(xp.begin() + i);
380 yp.erase(yp.begin() + i);
384 int chipSamp, chipLine;
386 chipSamp = startSamp;
392 chipLine = startLine;
398 int sampOffset = chipSamp - TackSample();
399 int lineOffset = chipLine - TackLine();
402 double matchChipSamp = match.
TackSample() + sampOffset;
403 double matchChipLine = match.
TackLine() + lineOffset;
408 if (matchCam != NULL) {
411 vector<int> newlocation = MovePoints(startSamp, startLine, endSamp, endLine);
412 startSamp = newlocation[0];
413 startLine = newlocation[1];
414 endSamp = newlocation[2];
415 endLine = newlocation[3];
423 if (!matchProj->
IsGood() ) {
424 vector<int> newlocation = MovePoints(startSamp, startLine, endSamp, endLine);
425 startSamp = newlocation[0];
426 startLine = newlocation[1];
427 endSamp = newlocation[2];
428 endLine = newlocation[3];
440 vector<int> newlocation = MovePoints(startSamp, startLine, endSamp, endLine);
441 startSamp = newlocation[0];
442 startLine = newlocation[1];
443 endSamp = newlocation[2];
444 endLine = newlocation[3];
453 vector<int> newlocation = MovePoints(startSamp, startLine, endSamp, endLine);
454 startSamp = newlocation[0];
455 startLine = newlocation[1];
456 endSamp = newlocation[2];
457 endLine = newlocation[3];
478 if (xp.size() == 3) {
482 if ( PointsColinear(xp[0], yp[0], xp[1], yp[1], xp[i], yp[i], tol) ) {
485 vector<int> newlocation = MovePoints(startSamp, startLine, endSamp, endLine);
486 startSamp = newlocation[0];
487 startLine = newlocation[1];
488 endSamp = newlocation[2];
489 endLine = newlocation[3];
497 QString msg =
"Cannot find enough points to perform Affine transformation. ";
498 msg +=
"Unable to load chip from [" + cube.
fileName();
499 msg +=
"] to match chip from [" + matchChipCube.
fileName() +
"].";
500 throw IException(IException::User, msg, _FILEINFO_);
504 m_affine.Solve(&x[0], &y[0], &xp[0], &yp[0], (
int)x.size());
507 m_affine.Scale(scale);
510 m_affine.Compute(0.0, 0.0);
511 double cubeSampleOffset = m_cubeTackSample - m_affine.xp();
512 double cubeLineOffset = m_cubeTackLine - m_affine.yp();
513 m_affine.Translate(cubeSampleOffset, cubeLineOffset);
546 bool Chip::PointsColinear(
const double x0,
const double y0,
const double x1,
const double y1,
547 const double x2,
const double y2,
const double tol) {
565 vector<double> v01, v12, v20;
566 v01.push_back(x1 - x0);
567 v01.push_back(y1 - y0);
568 v12.push_back(x2 - x1);
569 v12.push_back(y2 - y1);
570 v20.push_back(x0 - x2);
571 v20.push_back(y0 - y2);
575 double sinP0 = fabs(v01[0] * v20[1] - v01[1] * v20[0]) /
576 sqrt((pow(v01[0], 2) + pow(v01[1], 2)) * (pow(v12[0], 2) + pow(v12[1], 2)));
578 double sinP1 = fabs(v12[0] * v20[1] - v12[1] * v20[0]) /
579 sqrt((pow(v12[0], 2) + pow(v12[1], 2)) * (pow(v20[0], 2) + pow(v20[1], 2)));
581 double sinP2 = fabs(v20[0] * v01[1] - v20[1] * v01[0]) /
582 sqrt((pow(v20[0], 2) + pow(v20[1], 2)) * (pow(v01[0], 2) + pow(v01[1], 2)));
589 double minSinValue = min(sinP0, min(sinP1, sinP2));
590 if ( minSinValue < sin(tol *
PI / 180) ) {
615 vector<int> Chip::MovePoints(
const int startSamp,
const int startLine,
616 const int endSamp,
const int endLine) {
617 vector<int> newlocations(4);
618 int sinc = (endSamp - startSamp) / 4;
623 int linc = (endLine - startLine) / 3;
628 newlocations[0] = startSamp + sinc;
629 newlocations[1] = startLine + linc;
630 newlocations[2] = endSamp - sinc;
631 newlocations[3] = endLine - linc;
643 void Chip::SetChipPosition(
const double sample,
const double line) {
644 m_chipSample = sample;
646 m_affine.Compute(sample - TackSample(), line - TackLine());
647 m_cubeSample = m_affine.xp();
648 m_cubeLine = m_affine.yp();
660 void Chip::SetCubePosition(
const double sample,
const double line) {
661 m_cubeSample = sample;
663 m_affine.ComputeInverse(sample, line);
664 m_chipSample = m_affine.x() + TackSample();
665 m_chipLine = m_affine.y() + TackLine();
681 void Chip::SetValidRange(
const double minimum,
const double maximum) {
682 if (minimum >= maximum) {
683 QString msg =
"Unable to set valid chip range to [" +
toString(minimum);
684 msg +=
", " +
toString(maximum) +
"]. First parameter must be smaller than the second.";
685 throw IException(IException::Programmer, msg, _FILEINFO_);
688 m_validMinimum = minimum;
689 m_validMaximum = maximum;
701 bool Chip::IsValid(
double percentage) {
703 for (
int samp = 1; samp <= Samples(); samp++) {
704 for (
int line = 1; line <= Lines(); line++) {
705 if (IsValid(samp, line)) validCount++;
708 double validPercentage = 100.0 * (double) validCount /
709 (
double)(Samples() * Lines());
710 if (validPercentage < percentage)
return false;
727 Chip Chip::Extract(
int samples,
int lines,
int samp,
int line) {
728 if ( samples > Samples() || lines > Lines() ) {
729 QString msg =
"Cannot extract sub-chip of size [" +
toString(samples);
730 msg +=
", " +
toString(lines) +
"] from chip of size [" +
toString(Samples());
731 msg +=
", " +
toString(Lines()) +
"]";
732 throw IException(IException::Programmer, msg, _FILEINFO_);
735 Chip chipped(samples, lines);
736 for (
int oline = 1; oline <= lines; oline++) {
737 for (
int osamp = 1; osamp <= samples; osamp++) {
738 int thisSamp = samp + (osamp - chipped.
TackSample());
739 int thisLine = line + (oline - chipped.
TackLine());
740 if ((thisSamp < 1) ||
742 (thisSamp > Samples()) ||
743 (thisLine > Lines()) ) {
747 chipped.
SetValue(osamp, oline, GetValue(thisSamp, thisLine));
772 void Chip::Extract(
int samp,
int line,
Chip &chipped) {
773 int samples = chipped.
Samples();
774 int lines = chipped.
Lines();
779 for (
int oline = 1; oline <= lines; oline++) {
780 for (
int osamp = 1; osamp <= samples; osamp++) {
781 int thisSamp = samp + (osamp - chipped.
TackSample());
782 int thisLine = line + (oline - chipped.
TackLine());
783 if ((thisSamp < 1) ||
785 (thisSamp > Samples()) ||
786 (thisLine > Lines()) ) {
790 chipped.
SetValue(osamp, oline, GetValue(thisSamp, thisLine));
852 int samples = chipped.
Samples();
853 int lines = chipped.
Lines();
855 for (
int oline = 1; oline <= lines; oline++) {
856 int relativeLine = oline - chipped.
TackLine();
857 for (
int osamp = 1; osamp <= samples; osamp++) {
858 int relativeSamp = osamp - chipped.
TackSample();
859 affine.
Compute(relativeSamp, relativeLine);
860 double xp = affine.
xp() + TackSample();
861 double yp = affine.
yp() + TackLine();
863 for (
int i = 0; i < port.
size(); i++) {
864 int csamp = port.
Sample(i);
865 int cline = port.
Line(i);
868 (csamp > Samples()) ||
869 (cline > Lines()) ) {
873 port[i] = GetValue(csamp, cline);
912 stats->SetValidRange(m_validMinimum, m_validMaximum);
914 for (
int i = 0; i < m_chipSamples; i++) {
915 stats->
AddData(&m_buf[i][0], m_chipLines);
939 void Chip::Read(
Cube &cube,
const int band) {
945 for (
int line = 1; line <= Lines(); line++) {
946 for (
int samp = 1; samp <= Samples(); samp++) {
947 SetChipPosition((
double)samp, (
double)line);
948 if ((CubeSample() < 0.5) ||
949 (CubeLine() < 0.5) ||
953 m_buf[line-1][samp-1] = Isis::NULL8;
955 else if (m_clipPolygon == NULL) {
958 m_buf[line-1][samp-1] =
962 geos::geom::Point *pnt = globalFactory->createPoint(
963 geos::geom::Coordinate(CubeSample(), CubeLine()));
964 if (pnt->within(m_clipPolygon)) {
967 m_buf[line-1][samp-1] =
971 m_buf[line-1][samp-1] = Isis::NULL8;
985 void Chip::Write(
const QString &filename) {
990 for (
int i = 1; i <= Lines(); i++) {
992 for (
int j = 1; j <= Samples(); j++) {
993 line[j-1] = GetValue(j, i);
1008 void Chip::SetClipPolygon(
const geos::geom::MultiPolygon &clipPolygon) {
1009 if (m_clipPolygon != NULL)
delete m_clipPolygon;
1010 m_clipPolygon = PolygonTools::CopyMultiPolygon(clipPolygon);
1024 m_buf = other.
m_buf;
1040 if (m_clipPolygon) {
1041 delete m_clipPolygon;
1042 m_clipPolygon = NULL;
1046 m_clipPolygon = PolygonTools::MakeMultiPolygon(